From d87af81f4244b5cda3d8da13b10933db075c7583 Mon Sep 17 00:00:00 2001
From: Colin Lee <colin.lee@dal.ca>
Date: Fri, 4 Dec 2020 12:04:17 -0800
Subject: [PATCH] Forward and reverse models compile and run

---
 AdvCore_GridCompMod.F90 | 34 ++++++++++++++++++++++++++++++++++
 CubeToLatLon.F90        |  8 ++++----
 2 files changed, 38 insertions(+), 4 deletions(-)

diff --git a/AdvCore_GridCompMod.F90 b/AdvCore_GridCompMod.F90
index 1dff298d1..5d528974c 100755
--- a/AdvCore_GridCompMod.F90
+++ b/AdvCore_GridCompMod.F90
@@ -78,6 +78,10 @@ module AdvCore_GridCompMod
       logical     :: FV3_DynCoreIsRunning=.false.
       integer     :: AdvCore_Advection=1
       logical     :: chk_mass=.false.
+#ifdef ADJOINT
+      logical                    :: isAdjoint=.false.
+      character(len=ESMF_MAXSTR) :: modelPhase
+#endif
 
       integer,  parameter :: ntiles_per_pe = 1
 
@@ -355,6 +359,14 @@ subroutine SetServices(GC, rc)
       VERIFY_(STATUS)
       DT = ndt
 
+#ifdef ADJOINT
+      call MAPL_GetResource( MAPL, modelPhase, 'MODEL_PHASE:', default='FORWARD', RC=STATUS )
+      _VERIFY(STATUS)
+      isAdjoint = .false.
+      if (trim(ModelPhase) == 'ADJOINT') &
+           isAdjoint = .true.
+      if (isAdjoint) dt = -dt
+#endif
       ! Start up FV if AdvCore is running without FV3_DynCoreIsRunning
       !--------------------------------------------------
       if (.NOT. FV3_DynCoreIsRunning) then
@@ -556,6 +568,13 @@ subroutine Run(GC, IMPORT, EXPORT, CLOCK, RC)
       character(len=ESMF_MAXSTR), allocatable :: biggerlist(:)
       integer, parameter                  :: XLIST_MAX = 60
 
+#ifdef ADJOINT
+!     reverse time debug info
+      integer, parameter                :: DI = 3, DJ = 4, DL = 5
+      ! Debug variables
+      INTEGER, parameter             :: I_DBG = 6, J_DBG = 5, L_DBG=1
+#endif
+
 ! Get my name and set-up traceback handle
 ! ---------------------------------------
 
@@ -830,11 +849,22 @@ subroutine Run(GC, IMPORT, EXPORT, CLOCK, RC)
             endif
 
          endif
+#ifdef ADJOINT
+         if (.not. isAdjoint) &
+#endif
          firstRun=.false.
 
          ! Run FV3 advection
          !------------------
+#ifdef ADJOINT
+         if (AdvCore_Advection>0 .and. .not. firstRun) then
+         IF (MAPL_Am_I_Root()) THEN
+            WRITE(*,546) dt
+546         FORMAT(' calling offline_tracer_advection with timestep = ', f8.3)
+         ENDIF
+#else
          if (AdvCore_Advection>0) then
+#endif
          ! GCHP: use dry instead of moist pressure
          !call offline_tracer_advection(TRACERS, PLE0, PLE1, MFX, MFY, CX, CY, &
          call offline_tracer_advection(TRACERS, DryPLE0, DryPLE1, MFX, MFY, CX, CY, &
@@ -843,6 +873,10 @@ subroutine Run(GC, IMPORT, EXPORT, CLOCK, RC)
                                        FV_Atm(1)%domain, AK, BK, PTOP, FV_Atm(1)%npx, FV_Atm(1)%npy, FV_Atm(1)%npz,   &
                                        NQ, dt)
          endif
+#ifdef ADJOINT
+         if (isAdjoint) &
+              firstRun = .false.
+#endif
 
          ! Update tracer mass conservation
          !-------------------------------------------------------------------------
diff --git a/CubeToLatLon.F90 b/CubeToLatLon.F90
index 913ad951f..dcdf64851 100644
--- a/CubeToLatLon.F90
+++ b/CubeToLatLon.F90
@@ -1142,7 +1142,7 @@ subroutine C2LInterp(cubsph, latlon, index, weight, misval, subset, transpose)
                 ic=index(1,i,j)
                 jc=index(2,i,j)
 
-                ADJOINT: if(.not.transpose) then
+                ADJOINT_LOOP: if(.not.transpose) then
                    UNDEF: if( .not. present(misval)) then
                       latlon(ii,j) = weight(1,i,j)*cubsph(ic  ,jc  ,tile)  &
                                    + weight(2,i,j)*cubsph(ic  ,jc+1,tile)  &
@@ -1186,7 +1186,7 @@ subroutine C2LInterp(cubsph, latlon, index, weight, misval, subset, transpose)
                    cubsph(ic+1,jc+1,tile)=cubsph(ic+1,jc+1,tile)+weight(3,i,j)*latlon(ii,j)
                    cubsph(ic+1,jc  ,tile)=cubsph(ic+1,jc  ,tile)+weight(4,i,j)*latlon(ii,j)
 
-                end if ADJOINT
+                end if ADJOINT_LOOP
              endif HAVE_POINT
 
           enddo JLOOP
@@ -1239,7 +1239,7 @@ subroutine L2CInterp(latlon, cubsph, id1, id2, jdc, weight, misval, transpose)
              i1 = id1(i,j)
              i2 = id2(i,j)
 
-             ADJOINT: if(.not.transpose) then    
+             ADJOINT_LOOP: if(.not.transpose) then    
 
                 UNDEF: if(.not. present(misval)) then
                    cubsph(i,jx,k) = weight(1,i,j)*latlon(i1,j1  )   &
@@ -1282,7 +1282,7 @@ subroutine L2CInterp(latlon, cubsph, id1, id2, jdc, weight, misval, transpose)
                 latlon(i2,j1+1) = latlon(i2,j1+1) + weight(3,i,j)*cubsph(i,jx,k)
                 latlon(i1,j1+1) = latlon(i1,j1+1) + weight(4,i,j)*cubsph(i,jx,k)
 
-             end if ADJOINT
+             end if ADJOINT_LOOP
 
           enddo FACE_X
        enddo FACE_Y