Update to ROMS 2D Kernel and Implicit Vertical Advection

ROMS/Adjoint/ad_main3d.F

@@ -14,6 +14,14 @@ SUBROUTINE ad_main3d (RunInterval)
 !  backwards the adjoint model equations for all nested grids, if      !
 !  any, by the specified time interval (seconds), RunInterval.         !
 !                                                                      !
+# if defined STEP2D_FB_LF_AM3
+!  Numerical 2D time-stepping kernel: FB AB3-AM4                       !
+# elif defined STEP2D_FB_LF_AM3
+!  Numerical 2D time-stepping kernel: FB LF-AM3                        !
+# else
+!  Numerical 2D time-stepping kernel: LF-AM3 (Legacy scheme)           !
+# endif
+!                                                                      !
 !=======================================================================
 !
       USE mod_param
@@ -56,7 +64,7 @@ SUBROUTINE ad_main3d (RunInterval)
 # ifdef BVF_MIXING_NOT_YET
 !!    USE ad_bvf_mix_mod,          ONLY : ad_bvf_mix
 # endif
-      USE ad_diag_mod, ONLY : ad_diag
+      USE ad_diag_mod,             ONLY : ad_diag
 # if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
       USE ad_frc_adjust_mod,       ONLY : ad_frc_adjust
 # endif
@@ -119,7 +127,9 @@ SUBROUTINE ad_main3d (RunInterval)
 !!    USE ad_set_tides_mod,        ONLY : ad_set_tides
 # endif
       USE ad_set_vbc_mod,          ONLY : ad_set_vbc
+# if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
       USE ad_set_zeta_mod,         ONLY : ad_set_zeta
+# endif
       USE ad_step2d_mod,           ONLY : ad_step2d
 # ifndef TS_FIXED
       USE ad_step3d_t_mod,         ONLY : ad_step3d_t
@@ -135,15 +145,15 @@ SUBROUTINE ad_main3d (RunInterval)
 # ifdef TIDE_GENERATING_FORCES
       USE equilibrium_tide_mod,    ONLY : equilibrium_tide
 # endif
-# ifdef WEAK_CONSTRAINT
-      USE frc_weak_mod,            ONLY : frc_ADgather, frc_clear
-# endif
 # if defined ATM_COUPLING_NOT_YET && defined MCT_LIB
       USE mct_coupler_mod,         ONLY : ocn2atm_coupling
 # endif
 # if defined WAV_COUPLING_NOT_YET && defined MCT_LIB
       USE mct_coupler_mod,         ONLY : ocn2wav_coupling
 # endif
+# ifdef WEAK_CONSTRAINT
+      USE frc_weak_mod,            ONLY : frc_ADgather, frc_clear
+# endif
 # if (defined FOUR_DVAR    && !defined I4DVAR_ANA_SENSITIVITY) && \
       defined OBSERVATIONS
       USE obs_read_mod,            ONLY : obs_read
@@ -507,12 +517,7 @@ SUBROUTINE ad_main3d (RunInterval)
                   CALL ad_step3d_uv (ng, tile)
                 END DO
               END DO
-!
-!-----------------------------------------------------------------------
-!  Adjoint of recompute depths and thicknesses using the new time
-!  filtered free-surface. This call was moved from "ad_step2d" to here.
-!-----------------------------------------------------------------------
-!
+
 # ifdef NESTING
 !
 !  If nesting, determine vertical indices and vertical interpolation
@@ -523,12 +528,21 @@ SUBROUTINE ad_main3d (RunInterval)
                 CALL ad_nesting (ng, iADM, nzwgt)
               END DO
 # endif
+
+# if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
+!
+!-----------------------------------------------------------------------
+!  Adjoint of recompute depths and thicknesses using the new time
+!  filtered free-surface. This call was moved from "ad_step2d" to here.
+!-----------------------------------------------------------------------
+!
               DO ig=1,GridsInLayer(nl)
                 ng=GridNumber(ig,nl)
                 DO tile=last_tile(ng),first_tile(ng),-1
                   CALL ad_set_depth (ng, tile, iADM)
                 END DO
               END DO
+# endif
 
 # ifdef NESTING
 !
@@ -556,6 +570,127 @@ SUBROUTINE ad_main3d (RunInterval)
               END DO
 #  endif
 # endif
+
+# ifdef STEP2D_FB_AB3_AM4
+!
+!-----------------------------------------------------------------------
+!  Solve the vertically integrated primitive equations for the
+!  free-surface and barotropic momentum components using a generalized
+!  Forward-Backward, 3rd-order Adams-Bashforth / 4th-order Adams-Moulton
+!  (FB AB3-AM4) time stepping scheme (Shchepetkin and McWilliams, 2009).
+!-----------------------------------------------------------------------
+!
+              LOOP_2D : DO my_iif=MAXVAL(nfast),1,-1
+
+#  ifdef NESTING
+!
+!  If composite or mosaic grids, process additional points in the
+!  contact zone between connected grids for the state variables
+!  associated with the 2D engine Corrector step section.
+!  If refinement, check mass flux conservation between coarse and
+!  fine grids during debugging.
+!
+                DO ig=1,GridsInLayer(nl)
+                  ng=GridNumber(ig,nl)
+                  IF (ANY(CompositeGrid(:,ng))) THEN
+                    CALL ad_nesting (ng, iADM, n2dCS)
+                  END IF
+                END DO
+#  endif
+
+                DO ig=1,GridsInLayer(nl)
+                  ng=GridNumber(ig,nl)
+                  IF (my_iif.le.nfast(ng)) THEN
+                    iif(ng)=my_iif
+!^                  kstp(ng)=knew(ng)
+!^                  knew(ng)=kstp(ng)+1
+!^                  IF (knew(ng).gt.4) knew(ng)=1
+!^
+                    knew(ng)=kstp(ng)
+                    kstp(ng)=knew(ng)-1
+                    IF (kstp(ng).lt.1) kstp(ng)=4
+
+                    IF (MOD(knew(ng),2).eq.0) THEN            ! zig-zag
+                      DO tile=first_tile(ng),last_tile(ng),+1 ! processing
+                        CALL ad_step2d (ng, tile)             ! sequence
+                      END DO
+                    ELSE
+                      DO tile=last_tile(ng),first_tile(ng),-1
+                        CALL ad_step2d (ng, tile)
+                      END DO
+                    END IF
+                  END IF
+                END DO
+
+              END DO LOOP_2D
+# else
+
+#  ifdef STEP2D_FB_LF_AM3
+!
+!-----------------------------------------------------------------------
+!  Solve the vertically integrated primitive equations for the
+!  free-surface and barotropic momentum components using a predictor-
+!  corrector LeapFrog / 3rd-order Adams-Moulton with a Forward-Backward
+!  feeback (FB LF-AM3) stepping scheme (Shchepetkin and McWilliams,
+!  2009).
+!-----------------------------------------------------------------------
+!
+              LOOP_2D : DO my_iif=MAXVAL(nfast),1,-1
+
+#   ifdef NESTING
+!
+!  If composite or mosaic grids, process additional points in the
+!  contact zone between connected grids for the state variables
+!  associated with the 2D engine Corrector step section.
+!  If refinement, check mass flux conservation between coarse and
+!  fine grids during debugging.
+!
+                DO ig=1,GridsInLayer(nl)
+                  ng=GridNumber(ig,nl)
+                  IF (ANY(CompositeGrid(:,ng))) THEN
+                    CALL ad_nesting (ng, iADM, n2dCS)
+                  END IF
+                END DO
+#   endif
+!
+!  Corrector AM3 substep with FB-feedback.
+!
+                DO ig=1,GridsInLayer(nl)
+                  ng=GridNumber(ig,nl)
+                  IF (my_iif.le.nfast(ng)) THEN
+                    DO tile=first_tile(ng),last_tile(ng),+1
+                      CALL ad_step2d (ng, tile)
+                    END DO
+
+!^                  knew(ng)=3-kstp(ng)
+!^                  next_kstp(ng)=knew(ng)
+!^
+                    knew(ng)=next_kstp(ng)
+                    kstp(ng)=3-knew(ng)
+                  END IF
+                END DO
+!
+!  Predictor LF substep with FB-feedback.
+!
+                DO ig=1,GridsInLayer(nl)
+                  ng=GridNumber(ig,nl)
+                  IF (my_iif.le.nfast(ng)) THEN
+                    iif(ng)=my_iif
+                    DO tile=last_tile(ng),first_tile(ng),-1
+                      CALL ad_step2d (ng, tile)
+                    END DO
+
+!^                  kstp(ng)=next_kstp(ng)
+!^                  knew(ng)=3
+!^
+!^ HGA: Need to compute the adjoint of the indices
+
+                  END IF
+                END DO
+
+              END DO LOOP_2D
+
+#  else
 !
 !-----------------------------------------------------------------------
 !  Solve adjoint vertically integrated primitive equations for the
@@ -564,7 +699,7 @@ SUBROUTINE ad_main3d (RunInterval)
 !
               LOOP_2D : DO my_iif=MAXVAL(nfast)+1,1,-1
 
-# ifdef NESTING
+#   ifdef NESTING
 !
 !  If composite or mosaic grids, process additional points in the
 !  contact zone between connected grids for the state variables
@@ -581,7 +716,7 @@ SUBROUTINE ad_main3d (RunInterval)
                     CALL ad_nesting (ng, iADM, nmflx)
                   END IF
                 END DO
-# endif
+#   endif
 !
 !  Corrector step - Apply 2D adjoint time-step corrector scheme.  Notice
 !  ==============    that there is not need for a corrector step during
@@ -614,7 +749,7 @@ SUBROUTINE ad_main3d (RunInterval)
                   END IF
                 END DO
 
-# ifdef NESTING
+#   ifdef NESTING
 
 !
 !  If composite or mosaic grids, process additional points in the
@@ -632,7 +767,7 @@ SUBROUTINE ad_main3d (RunInterval)
                     CALL ad_nesting (ng, iADM, nmflx)
                   END IF
                 END DO
-# endif
+#   endif
 !
 !  Predictor step - Advance adjoint barotropic equations using 2D
 !  ==============   time-step predictor scheme.  No actual time-
@@ -670,7 +805,8 @@ SUBROUTINE ad_main3d (RunInterval)
                 END DO
 
               END DO LOOP_2D
-
+#  endif
+# endif
             END IF ADVANCE
 
 # if (defined FOUR_DVAR    && !defined I4DVAR_ANA_SENSITIVITY) && \
@@ -805,6 +941,9 @@ SUBROUTINE ad_main3d (RunInterval)
               END IF
             END DO
 # endif
+
+# if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3) || \
+       defined DIAGNOSTICS
 !
 !-----------------------------------------------------------------------
 !  Set adjoint free-surface to it time-averaged value.
@@ -813,12 +952,15 @@ SUBROUTINE ad_main3d (RunInterval)
             DO ig=1,GridsInLayer(nl)
               ng=GridNumber(ig,nl)
               DO tile=first_tile(ng),last_tile(ng),+1
-# ifdef DIAGNOSTICS
+#  ifdef DIAGNOSTICS
 !!              CALL set_diags (ng, tile)
-# endif
+#  endif
+#  if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
                 CALL ad_set_zeta (ng, tile)
+#  endif
               END DO
             END DO
+# endif
 !
 !-----------------------------------------------------------------------
 !  Compute adjoint vertical mixing coefficients for momentum and

ROMS/Adjoint/ad_set_data.F

@@ -1182,7 +1182,11 @@ SUBROUTINE ad_set_data_tile (ng, tile,                            &
 !
 !  Set forward free-surface.
 !
+#  if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
+      DO k=1,4
+#  else
       DO k=1,3
+#  endif
         CALL set_2dfldr_tile (ng, tile, iADM, idFsur,                   &
      &                        LBi, UBi, LBj, UBj,                       &
      &                        OCEAN(ng)%zetaG,                          &
@@ -1201,7 +1205,11 @@ SUBROUTINE ad_set_data_tile (ng, tile,                            &
 !
 !  Set forward 2D momentum.
 !
+#  if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
+      DO k=1,4
+#  else
       DO k=1,3
+#  endif
         CALL set_2dfldr_tile (ng, tile, iADM, idUbar,                   &
      &                        LBi, UBi, LBj, UBj,                       &
      &                        OCEAN(ng)%ubarG,                          &

ROMS/Adjoint/ad_set_zeta.F

@@ -1,7 +1,8 @@
 #include "cppdefs.h"
       MODULE ad_set_zeta_mod
 
-#if defined ADJOINT && defined SOLVE3D
+#if defined ADJOINT           && defined SOLVE3D           && \
+  !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
 !
 !git $Id$
 !================================================== Hernan G. Arango ===
@@ -85,7 +86,7 @@ SUBROUTINE ad_set_zeta_tile (ng, tile,                            &
       real(r8), intent(out) :: ad_zeta_sol(LBi:,LBj:)
 # else
       real(r8), intent(inout) :: ad_Zt_avg1(LBi:UBi,LBj:UBj)
-      real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,3)
+      real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
 
       real(r8), intent(out) :: ad_zeta_sol(LBi:UBi,LBj:UBj)
 # endif

ROMS/Adjoint/ad_step2d.F

@@ -8,16 +8,27 @@
 !    See License_ROMS.md                                               !
 !=======================================================================
 !                                                                      !
-!  This subroutine performs a fast (predictor or corrector) time-step  !
-!  for the free-surface and 2D momentum adjoint equations.             !
+!  This module timesteps the Adjoint Model (ADM) vertically-integrated !
+!  primitive (2D shallow-water) equations for the free-surface and 2D  !
+!  momentum. In 3D applications, the ROMS numerical kernel is split    !
+!  between baroclinic and barotropic dynamics. The barotropic engine   !
+!  uses a smaller timestep in this routine to resolve fast gravity     !
+!  wave processes.                                                     !
 # ifdef SOLVE3D
+!                                                                      !
 !  It also calculates the time filtering variables over all fast-time  !
 !  steps to damp high frequency signals in 3D applications.            !
 # endif
 !                                                                      !
 !=======================================================================
 !
-# include "ad_step2d_LF_AM3.h"
+# if defined STEP2D_FB_AB3_AM4
+#  include "ad_step2d_FB.h"
+# elif defined STEP2D_FB_LF_AM3
+#  include "ad_step2d_FB_LF_AM3.h"
+# else
+#  include "ad_step2d_LF_AM3.h"
+# endif
 #else
       MODULE ad_step2d_mod
       END MODULE ad_step2d_mod