Custom Query (986 matches)

Fixed couple of issues in the weak constraint 4D-Var algorithm:

• Removed some legacy code in routines forcing.F and tl_forcing.F when processing the impulse forcing for the weak constraint 4D-Var formulation. The legacy code that included the forcing for ubar and vbar in 3D solutions were removed from forcing.F when SOLVE3D is activated. Similarly, the forcing for tl_ubar and tl_vbar was removed from tl_forcing.F.

• Added logic to def_impulse.F, get_data.F, set_data.F, tl_get_data.F, tl_set_data.F, and checkvars.F to avoid processing the weak constraint impulse forcing for ubar (tl_ubar) and vbar (tl_vbar) when SOLVE3D is activated.

Many thanks to Andy Moore for his help in tracking this bug.

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.