Custom Query (964 matches)

Corrected critical bug in some of the tracer advection options:

• Corrected a bug in the HSIMT advection. Previously, I fixed the issue that three ghost points are needed for this scheme. At that time, I missed making a similar change in the computation of the inverse thickness (oHz). Therefore, in step3d_t.F we need to have

        logical :: LapplySrc, Lhsimt, Lmpdata
        ...
        
        Lhsimt =ANY(Hadvection(:,ng)%HSIMT).and.                          &
       &        ANY(Vadvection(:,ng)%HSIMT)
        Lmpdata=ANY(Hadvection(:,ng)%MPDATA).and.                         &
       &        ANY(Vadvection(:,ng)%MPDATA)
        ...
        
  !
  !  Compute reciprocal thickness, 1/Hz.
  !
        IF (Lmpdata.or.Lhsimt) THEN
          DO k=1,N(ng)
            DO j=Jstrm2,Jendp2
              DO i=Istrm2,Iendp2
                oHz(i,j,k)=1.0_r8/Hz(i,j,k)
              END DO
            END DO
          END DO
        ELSE
          DO k=1,N(ng)
            DO j=Jstr,Jend
              DO i=Istr,Iend
                oHz(i,j,k)=1.0_r8/Hz(i,j,k)
              END DO
            END DO
          END DO
        END IF
  

  I have been hunting on and off for a bug in the HSIMT scheme for awhile. Many thanks to Pierre St-Laurent for ​reporting this issue. It is always helpful when another set of eyes look into a code.

• Corrected a weird bug in the tangent linear advection switches. In inp_par.F we need to have instead:

  #if defined TANGENT || defined TL_IOMS
  !
  !  Set tracer advection scheme switches for the tangent linear models
  !  (TLM and RPM) to the same values as the adjoint model.
  !
        DO ng=1,Ngrids
          DO i=1,NT(ng)
            tl_Hadvection(i,ng)%AKIMA4    = ad_Hadvection(i,ng)%AKIMA4
            tl_Hadvection(i,ng)%CENTERED2 = ad_Hadvection(i,ng)%CENTERED2
            tl_Hadvection(i,ng)%CENTERED4 = ad_Hadvection(i,ng)%CENTERED4
            tl_Hadvection(i,ng)%HSIMT     = ad_Hadvection(i,ng)%HSIMT
            tl_Hadvection(i,ng)%MPDATA    = ad_Hadvection(i,ng)%MPDATA
            tl_Hadvection(i,ng)%SPLINES   = ad_Hadvection(i,ng)%SPLINES
            tl_Hadvection(i,ng)%SPLIT_U3  = ad_Hadvection(i,ng)%SPLIT_U3
            tl_Hadvection(i,ng)%UPSTREAM3 = ad_Hadvection(i,ng)%UPSTREAM3
  !
            tl_Vadvection(i,ng)%AKIMA4    = ad_Vadvection(i,ng)%AKIMA4
            tl_Vadvection(i,ng)%CENTERED2 = ad_Vadvection(i,ng)%CENTERED2
            tl_Vadvection(i,ng)%CENTERED4 = ad_Vadvection(i,ng)%CENTERED4
            tl_Vadvection(i,ng)%HSIMT     = ad_Vadvection(i,ng)%HSIMT
            tl_Vadvection(i,ng)%MPDATA    = ad_Vadvection(i,ng)%MPDATA
            tl_Vadvection(i,ng)%SPLINES   = ad_Vadvection(i,ng)%SPLINES
            tl_Vadvection(i,ng)%SPLIT_U3  = ad_Vadvection(i,ng)%SPLIT_U3
            tl_Vadvection(i,ng)%UPSTREAM3 = ad_Vadvection(i,ng)%UPSTREAM3
          END DO
        END DO
  #endif
  

  Before, the assigment was to use the nonlinear values instead. As a consequence, the 4D-Var cost function was not decreasing correctly and the minimization failed because the broken symmetry between tangent linear (TLM) and adjoint (ADM) model operators. We need to have the same advection scheme in both TLM and ADM. Recall that in roms.in, the user specify just the values for the adjoint model.

Many thanks to Julia Levin for reporting this problem and Andy Moore for his help in debugging.

This update includes a couple of changes:

• Fixed the tracer advection logic in inp_par.F and tadv.F for the shallow-water (SOLVE3D undefined) applications not needing tracer advection.

  !
  !  Set switch for three ghost-points in the halo region.
  !
  #ifdef SOLVE3D
        ThreeGhostPoints=ANY(Hadvection(:,:)%MPDATA).or.                  &
       &                 ANY(Hadvection(:,:)%HSIMT)
  #endif
  #ifdef UV_VIS4
        ThreeGhostPoints=.TRUE.
  #endif
  

Many thanks to Guangyu Xu for ​reporting this issue.

• Updated inp_decode.F function load_s1d to have the following interface:

       INTERFACE load_s1d
         MODULE PROCEDURE load_s1d1       ! 1D structrure, S(:)
         MODULE PROCEDURE load_s1d2       ! 2D structrure, S(Ie,:)
       END INTERFACE load_s1d
  

  so we have specific funtions for 1D (load_s1d1) and 2D (load_s1d2) TYPE_IO structures. For example, in read_asspar.F we have:

             CASE ('STDnameI')
               label='STD - initial conditions standard deviation'
               Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid,       &
      &                      Ngrids, Nfiles, 4, 1, STD)
             CASE ('STDnameM')
               label='STD - model error standard deviation'
               Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid,       &
      &                      Ngrids, Nfiles, 4, 2, STD)
             CASE ('STDnameB')
               label='STD - boundary conditions standard deviation'
               Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid,       &
      &                      Ngrids, Nfiles, 4, 3, STD)
             CASE ('STDnameF')
               label='STD - surface forcing standard deviation'
               Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid,       &
      &                      Ngrids, Nfiles, 4, 4, STD)
  

  to process each element of the STD(1:4,Ngrids) standard deviation structure used in 4D-Var.

I have several research repositories for several ongoing new algorithms that will be available in the future. I have so many exceptional flags to update the codes that sometimes one fails, and undesired changes appear in public released version of ROMS and erase features in the research versions. To avoid doing that in the future, I am consolidating the research 4DVAR repository with the public version of the code.

For the last two months, we have been developing the saddle-point 4D-Var. By far, one of the most technically tricky algorithms that Andy and I have put our hands on. As formulated, the current available 4D-Var algorithms are sequential, and its inner loops cannot be parallelized in time to accelerate the analysis. The conjugate gradient solver used to minimize the cost function J, depends sequentially on the previous iteration.

Fisher and Gürol (2017) proposed a new weak constraint 4D-Var algorithm based on the saddle-point formulation that can be parallelized in time. The time trajectory can be split into the solution of Nsaddle intervals that can run concurrently to reduce the 4D-Var computational cost.

To achieve these technically challenging computations, we need to make changes to several ROMS routines to:

• To split the MPI communicator to carry out disjointed and concurrent runs of the nonlinear, tangent linear, and adjoint ROMS kernels. Therefore, several of the routines of distribute.F have an optional argument to pass the desired MPI communicators that are defined in mod_parallel.F.

• Add optional arguments to output NetCDF routines nf_fwrite2d.F, nf_fwrite2d_bry.F, nf_fwrite3d.F, nf_fwrite3d_bry.F, and nf_fwrite4d.F to get the minimun and maximun values ot the written fields.

• Include several new generic routines to manipulate ROMS state vector: ADfromTL.F, def_state.F, state_join.F, state_read.F, and wrt_state.F.

• No change is made to ROMS nonlinear, tangent linear, and adjoint time-stepping kernels.

• Added several routines from the LAPACK library that are needed for the Generalized Minimal Residual (GMRES) minimization solver.

• The saddle-point algorithm with be released in the future when it is fully tested and we have the proper documentation and relevant publication.