VERY IMPORTANT: Miscelaneous Update

[arango]

This update enhances several routines:

• The routine ad_main3d is rewritten to include the nesting logic. The 4D-Var nesting scheme are working and will be released sometime next year. The complicated nested 4D-Var algorithms were developed in collaboration with Andy Moore. It was a herculean and tricky effort involving lots of testing.

• For consistency, the time-stepping counter is updated at the bottom of main3d, tl_main3d, and ad_main3d. For example, in main3d.F we have:

  !
  !-----------------------------------------------------------------------
  !  Advance time index and time clock.
  !-----------------------------------------------------------------------
  !
              DO ig=1,GridsInLayer(nl)
                ng=GridNumber(ig,nl)
                iic(ng)=iic(ng)+1
                time(ng)=time(ng)+dt(ng)
                step_counter(ng)=step_counter(ng)-1
                CALL time_string (time(ng), time_code(ng))
              END DO
  

  It cleaned the awkward logic at the bottom of initial, tl_initial, and ad_initial. This change was required for a more straightforward nested 4D-Var design.

• Adds AD_OUTPUT_STATE option to process and write out the whole adjoint solution. Due to the exact discrete adjoint, the predictor/corrector time-stepping scheme with multiple time levels, pieces of the adjoint solution are in two-time levels and need to be added in the ad_*_sol arrays for output purposes. We need this capability in the stochastic algorithms and the ROMS-JEDI interface.

Thus, changes are made to ad_ini_fields, ad_main2d, ad_main3d, and ad_step3d_uv. The two pieces for ad_u_sol and ad_v_sol are computed in ad_step3d_uv (index nnew). The other piece is computed in ad_pre_step3d, and added in ad_out_fields:

!
!  Add a second piece of the 3D adjoint momentum solution into IO arrays.
!  The first is loaded in "ad_step3d_uv".
!
      DO j=JstrB,JendB
        IF (j.ge.JstrM) THEN
          DO k=1,N(ng)
            DO i=IstrB,IendB
              ad_v_sol(i,j,k)=ad_v_sol(i,j,k)+ad_v(i,j,k,nnew)
            END DO
          END DO
        END IF
        DO k=1,N(ng)
          DO i=IstrM,IendB
            ad_u_sol(i,j,k)=ad_u_sol(i,j,k)+ad_u(i,j,k,nnew)
          END DO
        END DO
      END DO

• The time index logic in ad_wrt_his was slightly modified to account for the change of location of iic(ng) counter update in ad_main3d:

  !
  !  Determine time index to write. The "nout" index is updated to the
  !  version of "ad_main3d" that updates the "iic" counter at the bottom.
  !  Therefore, we need to change the conditional "iic(ng).ne.ntend(ng)"
  !  to "iic(ng).gt.ntend(ng)" to get identical solutions.
  !
  # ifdef SOLVE3D
        kout=kstp(ng)
  # else
        kout=kstp(ng)
  # endif
  # if defined WEAK_CONSTRAINT
        kfout=2
  # endif
  # ifdef SOLVE3D
        IF (iic(ng).gt.ntend(ng)) THEN
          nout=nnew(ng)
  # ifdef AD_OUTPUT_STATE
          LwrtState3d(ng)=.FALSE.
  # endif
        ELSE
  # ifdef AD_OUTPUT_STATE
          LwrtState3d(ng)=.TRUE.
  # endif
          nout=nstp(ng)
        END IF
  # endif
  

• The jedi_roms.h driver was modified to facilitate the LinearModel Class for the ROMS-JEDI interface. OOPS will advance the ROMS kernels by smaller intervals, usually a single timestep. The strategy here is different from that used for coupling because ROMS delayed output. The delayed output last half timestep will affect the OOPS trajectory logic needed to save the NLM background fields required to linearize the TLM and ADM kernels. Therefore, we need to address some technical issues for the correct TL and AD state vector to exchange between ROMS and OOPS.

• Enhanced function netcdf_inq_var in mod_netcdf.F to include more information about file global dimensions: names, ID, and size.

• Updated the lateral boundary conditions for the vertical mixing coefficients in the turbulent parameterizations schemes. For example, in gls_corstep.F, we have:

       DO k=0,N(ng)
         IF (DOMAIN(ng)%Western_Edge(tile)) THEN
           DO j=Jstr,Jend
             DO itrc=1,NAT
               Akt(Istr-1,j,k,itrc)=Akt(Istr,j,k,itrc)
             END DO
             Akv(Istr-1,j,k)=Akv(Istr,j,k)
           END DO
         END IF
         IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
           DO j=Jstr,Jend
             DO itrc=1,NAT
               Akt(Iend-1,j,k,itrc)=Akt(Iend,j,k,itrc)
             END DO
             Akv(Iend-1,j,k)=Akv(Iend,j,k)
           END DO
  
         END IF
         IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
           DO i=Istr,Iend
             DO itrc=1,NAT
               Akt(i,Jstr-1,k,itrc)=Akt(i,Jstr,k,itrc)
             END DO
             Akv(i,Jstr-1,k)=Akv(i,Jstr,k)
           END DO
         END IF
         IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
           DO i=Istr,Iend
             DO itrc=1,NAT
               Akt(i,Jend+1,k,itrc)=Akt(i,Jend,k,itrc)
             END DO
             Akv(i,Jend+1,k)=Akv(i,Jend,k)
           END DO
         END IF
         IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
           DO itrc=1,NAT
             Akt(Istr-1,Jstr-1,k,itrc)=0.5_r8*                           &
      &                                (Akt(Istr  ,Jstr-1,k,itrc)+       &
      &                                 Akt(Istr-1,Jstr  ,k,itrc))
           END DO
           Akv(Istr-1,Jstr-1,k)=0.5_r8*                                  &
      &                         (Akv(Istr  ,Jstr-1,k)+                   &
      &                          Akv(Istr-1,Jstr  ,k))
         END IF
         IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
           DO itrc=1,NAT
             Akt(Iend+1,Jstr-1,k,itrc)=0.5_r8*                           &
      &                                (Akt(Iend  ,Jstr-1,k,itrc)+       &
      &                                 Akt(Iend+1,Jstr  ,k,itrc))
           END DO
           Akv(Iend+1,Jstr-1,k)=0.5_r8*                                  &
      &                         (Akv(Iend  ,Jstr-1,k)+                   &
      &                          Akv(Iend+1,Jstr  ,k))
         END IF
         IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
           DO itrc=1,NAT
             Akt(Istr-1,Jend+1,k,itrc)=0.5_r8*                           &
      &                                (Akt(Istr  ,Jend+1,k,itrc)+       &
      &                                 Akt(Istr-1,Jend  ,k,itrc))
           END DO
           Akv(Istr-1,Jend+1,k)=0.5_r8*                                  &
      &                         (Akv(Istr  ,Jend+1,k)+                   &
      &                          Akv(Istr-1,Jend  ,k))
         END IF
         IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
           DO itrc=1,NAT
             Akt(Iend+1,Jend+1,k,itrc)=0.5_r8*                           &
      &                                (Akt(Iend  ,Jend+1,k,itrc)+       &
      &                                 Akt(Iend+1,Jend  ,k,itrc))
           END DO
           Akv(Iend+1,Jend+1,k)=0.5_r8*                                  &
      &                         (Akv(Iend  ,Jend+1,k)+                   &
      &                          Akv(Iend+1,Jend  ,k))
         END IF
       END DO
  
  

  The parameterization schemes compute the vertical mixing at the interior points, and we have boundary conditions for only periodic applications. This issue was discovered in the ROMS-JEDI interface. This change introduces minimal roundoff changes to the nonlinear model solution, which affects the TLM, RPM, ADM solutions.

• Corrected bug in def_avg_pio that gives error when synchronizing a PIO file to disk. Many thanks to Dave Robertson for bringing this issue to my attention.