From 277b5bf04d136f34a3f2856d8b2f5ff0ced309e1 Mon Sep 17 00:00:00 2001
From: Purnendu Chakraborty <purnendu.chakraborty@nasa.gov>
Date: Fri, 12 May 2023 19:46:15 -0400
Subject: [PATCH 01/13] Added GEOShs_GridComp and associated changes

---
 CMakeLists.txt                                |    5 +-
 GEOSagcm_GridComp/CMakeLists.txt              |   10 +-
 GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90 |  389 ++++++
 .../GEOShs_GridComp/CMakeLists.txt            |   13 +
 .../GEOShs_GridComp/GEOS_HsGridComp.F90       | 1129 +++++++++++++++++
 GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc   |    2 +
 .../GEOShs_GridComp/held_suarez_oacc.F90      |  238 ++++
 .../GEOSsuperdyn_GridComp/CMakeLists.txt      |    6 +
 8 files changed, 1790 insertions(+), 2 deletions(-)
 create mode 100644 GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
 create mode 100644 GEOSagcm_GridComp/GEOShs_GridComp/CMakeLists.txt
 create mode 100644 GEOSagcm_GridComp/GEOShs_GridComp/GEOS_HsGridComp.F90
 create mode 100755 GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc
 create mode 100644 GEOSagcm_GridComp/GEOShs_GridComp/held_suarez_oacc.F90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4c9699c3e..fecd215ba 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,8 +7,11 @@ set (alldirs
   GEOSogcm_GridComp
   )
 
+if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90)
 
-if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_GcmGridComp.F90)
+  esma_add_subdirectories (GEOSagcm_GridComp)
+
+elseif (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_GcmGridComp.F90)
 
  ecbuild_declare_project()
   esma_add_library(${this}
diff --git a/GEOSagcm_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/CMakeLists.txt
index 5cbdbb1c7..734df590f 100644
--- a/GEOSagcm_GridComp/CMakeLists.txt
+++ b/GEOSagcm_GridComp/CMakeLists.txt
@@ -5,7 +5,15 @@ set (alldirs
   GEOSphysics_GridComp
   )
 
-if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_AgcmGridComp.F90)
+if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_AgcmSimpleGridComp.F90)
+
+  message(STATUS "Building Held-Suarez GridComp")
+  esma_add_library (${this}
+    SRCS GEOS_AgcmSimpleGridComp.F90
+    SUBCOMPONENTS GEOSsuperdyn_GridComp GEOShs_GridComp
+    DEPENDENCIES MAPL esmf)
+
+elseif (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_AgcmGridComp.F90)
 
   esma_add_library (${this}
     SRCS GEOS_AgcmGridComp.F90
diff --git a/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90 b/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
new file mode 100644
index 000000000..66c5099cd
--- /dev/null
+++ b/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
@@ -0,0 +1,389 @@
+! $Id: GEOS_AgcmSimpleGridComp.F90,v 1.1 2007/05/16 15:33:09 trayanov Exp $
+
+#include "MAPL_Generic.h"
+
+
+module GEOS_AgcmSimpleGridCompMod
+
+!=============================================================================
+!BOP
+! \renewcommand{\comp}{\tt GEOS\_AgcmSimpleGridCompMod}
+!
+! !MODULE: GEOS_AgcmSimpleGridCompMod
+!
+! !DESCRIPTION: This gridded component (GC) combines the the GC that
+!   implements the Finite-Volume (FV) dynamics, with a simple physics
+!   component that implements the Held-Suarez benchmark forcing for 
+!   testing dry dynamical cores. 
+!
+! !USES:
+
+  use ESMF
+  use MAPL_Mod
+  use GEOS_TopoGetMod
+
+  use GEOS_superdynGridCompMod,  only:  SDYN_SetServices => SetServices
+  use GEOS_hsGridCompMod,     only:  PHS_SetServices => SetServices
+
+  implicit none
+  private
+
+  integer :: SDYN
+  integer :: PHS
+
+! !PUBLIC MEMBER FUNCTIONS:
+
+  public SetServices
+
+ 
+!EOP
+
+contains
+
+!BOP
+
+! !IROUTINE: SetServices
+
+! !DESCRIPTION:  This is the only method of this component, since all of its
+!   registered methods can be defaulted to their MAPL generic versions.
+!   SetServices merely creates the children through MAPL and connects
+!   their Import-Export states. FV has some Imports that are
+!   not used in the Held-Suarez benchmark. Since these have proper defaults
+!   in FV, they can simply be terminated here. 
+
+
+! !INTERFACE:
+
+    subroutine SetServices ( GC, RC )
+
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(INOUT) :: GC  ! gridded component
+    integer, optional,   intent(  OUT) :: RC  ! return code
+
+
+!EOP
+!BOC
+!=============================================================================
+!
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)              :: IAm
+    integer                                 :: STATUS
+    character(len=ESMF_MAXSTR)              :: COMP_NAME
+
+! Locals
+
+
+!=============================================================================
+
+! Begin...
+
+! Get my name and set-up traceback handle
+! ---------------------------------------
+
+    Iam = 'SetServices'
+    call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // trim(Iam)
+
+    call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS )
+    VERIFY_(STATUS)
+
+    call MAPL_GridCompSetEntryPoint ( gc, ESMF_METHOD_RUN,  Run, __RC__)
+
+! dummy import for testing comcurrent ens
+!    call MAPL_AddImportSpec ( gc,                                  &
+!         SHORT_NAME = 'DTDT',                                      &
+!         LONG_NAME  = 'temperature increment',                     &
+!         UNITS      = 'K s-1',                                     &
+!         DIMS       = MAPL_DimsHorzVert,                           &
+!         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+!    VERIFY_(STATUS)
+
+! these internal spec are "fake" and here only to provide moisture to FV
+! this component will not touch them
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'Q',                                         &
+         LONG_NAME  = 'specific_humidity',                         &
+         UNITS      = 'kg kg-1',                                   &
+         FRIENDLYTO = 'DYNAMICS',    &
+         default    = 1.0e-6,                                      &
+!         RESTART    = MAPL_RestartRequired,                        &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                        &
+         SHORT_NAME = 'QLLS',                                            &
+         LONG_NAME  = 'mass_fraction_of_large_scale_cloud_liquid_water', &
+         UNITS      = 'kg kg-1',                                         &
+         FRIENDLYTO = 'DYNAMICS',                             &
+         DIMS       = MAPL_DimsHorzVert,                                 &
+         VLOCATION  = MAPL_VLocationCenter,                   RC=STATUS  )  
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                       &
+         SHORT_NAME = 'QLCN',                                           &
+         LONG_NAME  = 'mass_fraction_of_convective_cloud_liquid_water', &
+         UNITS      = 'kg kg-1',                                        &
+         FRIENDLYTO = 'DYNAMICS',                            &
+         DIMS       = MAPL_DimsHorzVert,                                &
+         VLOCATION  = MAPL_VLocationCenter,                  RC=STATUS  )  
+    VERIFY_(STATUS)                                                                          
+    call MAPL_AddInternalSpec(GC,                                     &
+         SHORT_NAME = 'QILS',                                         &
+         LONG_NAME  = 'mass_fraction_of_large_scale_cloud_ice_water', &
+         UNITS      = 'kg kg-1',                                      &
+         FRIENDLYTO = 'DYNAMICS',                          &
+         DIMS       = MAPL_DimsHorzVert,                              &
+         VLOCATION  = MAPL_VLocationCenter,                RC=STATUS  )  
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                    &
+         SHORT_NAME = 'QICN',                                        &
+         LONG_NAME  = 'mass_fraction_of_convective_cloud_ice_water', &
+         UNITS      = 'kg kg-1',                                     &
+         FRIENDLYTO = 'DYNAMICS',                         &
+         DIMS       = MAPL_DimsHorzVert,                             &
+         VLOCATION  = MAPL_VLocationCenter,               RC=STATUS  )  
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'CLLS',                                      &
+         LONG_NAME  = 'large_scale_cloud_area_fraction',           &
+         UNITS      = '1',                                         &
+         FRIENDLYTO = 'DYNAMICS',                                  &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
+    VERIFY_(STATUS)                                                                          
+
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'CLCN',                                      &
+         LONG_NAME  = 'convective_cloud_area_fraction',            &
+         UNITS      = '1',                                         &
+         FRIENDLYTO = 'DYNAMICS',                                  &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
+    VERIFY_(STATUS)                                                                          
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'QRAIN',                                     &
+         LONG_NAME  = 'mass_fraction_of_rain',                     & 
+         UNITS      = 'kg kg-1',                                   &
+         FRIENDLYTO = 'DYNAMICS',                       &
+         default    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'QSNOW',                                     &
+         LONG_NAME  = 'mass_fraction_of_snow',                     &
+         UNITS      = 'kg kg-1',                                   &
+         FRIENDLYTO = 'DYNAMICS',                       &
+         default    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(GC,                                  &
+         SHORT_NAME = 'QGRAUPEL',                                  &
+         LONG_NAME  = 'mass_fraction_of_graupel',                  &
+         UNITS      = 'kg kg-1',                                   &
+         FRIENDLYTO = 'DYNAMICS',                       &
+         default    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+    VERIFY_(STATUS)
+                                                                          
+
+! Register children with MAPL and go down their SS hierarchy
+! ----------------------------------------------------------
+
+    SDYN  = MAPL_AddChild(GC, NAME='SUPERDYNAMICS', SS=SDYN_SetServices, RC=STATUS)
+    VERIFY_(STATUS)
+    PHS  = MAPL_AddChild(GC, NAME='HSPHYSICS',  SS=PHS_SetServices, RC=STATUS)
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec ( GC   ,                          &
+         SHORT_NAME = 'T',                                    &
+         CHILD_ID   = SDYN,                                    &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec ( GC   ,                          &
+         SHORT_NAME = 'PS',                                   &
+         CHILD_ID   = SDYN,                                    &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+! Register connections between children
+! -------------------------------------
+
+    call MAPL_AddConnectivity ( GC,                                    &
+         SHORT_NAME  = (/ 'DUDT', 'DVDT', 'DTDT' /),                   &
+         SRC_ID      =  PHS,                                           &
+         DST_ID      =  SDYN,                                           &
+                                                            RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddConnectivity ( GC,                                    &
+         SRC_NAME = (/ 'U    ', 'V    ', 'T    ', 'PLE  ' /),          &
+         DST_NAME = (/ 'U    ', 'V    ', 'TEMP ', 'PLE  ' /),          &
+         SRC_ID   =  SDYN,                                              &
+         DST_ID   =  PHS,                                              &
+                                                            RC=STATUS  )
+    VERIFY_(STATUS)
+
+
+! SetServices clean-up on the way back up through the hierarchy
+!--------------------------------------------------------------
+
+    call MAPL_TerminateImport(GC, SHORT_NAME = (/'PHIS ','DPEDT'/), &
+         CHILD = SDYN, RC=STATUS)
+    VERIFY_(STATUS)
+
+    call MAPL_GenericSetServices( GC, RC=STATUS )
+    VERIFY_(STATUS)
+   
+    RETURN_(ESMF_SUCCESS)
+  end subroutine SetServices
+!EOC
+
+  subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC )
+
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_State),    intent(inout) :: IMPORT ! Import state
+    type(ESMF_State),    intent(inout) :: EXPORT ! Export state
+    type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
+    integer, optional,   intent(  out) :: RC     ! Error code
+
+! !DESCRIPTION: The Initialize method of this Gridded Component.
+
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)              :: IAm 
+    integer                                 :: STATUS
+    character(len=ESMF_MAXSTR)              :: COMP_NAME
+
+! local vars
+   type (MAPL_MetaComp),      pointer  :: MAPL
+   type (ESMF_State),         pointer  :: GIM(:)
+   type (ESMF_Field)                   :: FIELD
+   type (ESMF_FieldBundle)             :: BUNDLE
+   type (ESMF_Config)                  :: cf
+
+    Iam = "Initialize"
+    call ESMF_GridCompGet ( GC, name=COMP_NAME, config=cf, RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // Iam
+
+    call MAPL_GridCreate(GC, rc=status)
+    VERIFY_(STATUS)
+
+    call MAPL_GenericInitialize ( GC, IMPORT, EXPORT, CLOCK,  RC=STATUS)
+    VERIFY_(STATUS)
+
+! Get my MAPL_Generic state
+!--------------------------
+
+    call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS)
+    VERIFY_(STATUS)
+
+! Fill Childrens TOPO variables and Diagnostics
+!----------------------------------------------
+    call MAPL_Get ( MAPL, GIM=GIM, RC=STATUS )
+    VERIFY_(STATUS)
+
+! PHIS ...
+!---------
+    call ESMF_StateGet( GIM(SDYN), 'PHIS', FIELD, rc=STATUS )
+    VERIFY_(STATUS)
+    Call GEOS_TopoGet ( cf, MEAN=FIELD, rc=STATUS )
+    VERIFY_(STATUS)
+
+! TRADV ...
+!----------
+    call ESMF_StateGet(GIM(SDYN), 'TRADV', BUNDLE, RC=STATUS )
+    VERIFY_(STATUS)
+    call MAPL_GridCompGetFriendlies(GC, "DYNAMICS", BUNDLE, RC=STATUS )
+    VERIFY_(STATUS)
+
+    RETURN_(ESMF_SUCCESS)
+  end subroutine Initialize
+
+  subroutine Run( GC, IMPORT, EXPORT, CLOCK, RC )
+
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_State),    intent(inout) :: IMPORT ! Import state
+    type(ESMF_State),    intent(inout) :: EXPORT ! Export state
+    type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
+    integer, optional,   intent(  out) :: RC     ! Error code
+
+ 
+!EOP
+
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)        :: IAm 
+    integer                           :: STATUS
+    character(len=ESMF_MAXSTR)        :: COMP_NAME
+    
+! Local derived type aliases
+
+    type (MAPL_MetaComp    ), pointer :: MAPL
+    type (ESMF_GridComp),      pointer  :: GCS(:)
+    type (ESMF_State),         pointer  :: GIM(:)
+    type (ESMF_State),         pointer  :: GEX(:)
+    
+!=============================================================================
+
+! Begin... 
+
+! Get the target components name and set-up traceback handle.
+! -----------------------------------------------------------
+
+    Iam = "Run"
+    call ESMF_GridCompGet( GC, name=COMP_NAME, RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // trim(Iam)
+
+! Get my MAPL_MetaComp
+!---------------------
+
+    call MAPL_GetObjectFromGC( GC, MAPL, RC=STATUS)
+    VERIFY_(STATUS)
+
+! Start the TOTAL timer
+!----------------------
+
+    call MAPL_TimerOn(MAPL,"TOTAL")
+    call MAPL_TimerOn(MAPL,"RUN")
+
+! Get esmf internal state from generic state.
+!-------------------------------------------
+
+    call MAPL_Get(MAPL, GCS=GCS, GIM=GIM, GEX=GEX, __RC__ )
+
+    call ESMF_GridCompRun(GCS(SDYN), importState=GIM(SDYN), exportState=GEX(SDYN), clock=CLOCK, PHASE=1, userRC=STATUS)
+    VERIFY_(STATUS)
+
+    call MAPL_GenericRunCouplers( MAPL, CHILD=SDYN, CLOCK=clock, __RC__ )
+
+    call ESMF_GridCompRun(GCS(PHS), importState=GIM(PHS), exportState=GEX(PHS), clock=CLOCK, userRC=STATUS)
+    VERIFY_(STATUS)
+
+    call ESMF_GridCompRun(GCS(SDYN), importState=GIM(SDYN), exportState=GEX(SDYN), clock=CLOCK, PHASE=2, userRC=STATUS)
+    VERIFY_(STATUS)
+
+    call MAPL_TimerOff(MAPL,"RUN")
+    call MAPL_TimerOff(MAPL,"TOTAL")
+
+    RETURN_(ESMF_SUCCESS)
+  end subroutine Run
+end module GEOS_AgcmSimpleGridCompMod
+
diff --git a/GEOSagcm_GridComp/GEOShs_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOShs_GridComp/CMakeLists.txt
new file mode 100644
index 000000000..76e93eb3b
--- /dev/null
+++ b/GEOSagcm_GridComp/GEOShs_GridComp/CMakeLists.txt
@@ -0,0 +1,13 @@
+esma_set_this()
+
+set (srcs
+  held_suarez_oacc.F90
+  GEOS_HsGridComp.F90
+  )
+
+esma_add_library (${this}
+  SRCS ${srcs}
+  DEPENDENCIES GEOS_Shared MAPL)
+target_include_directories (${this} PUBLIC ${INC_ESMF} ${INC_NETCDF})
+
+target_compile_options(${this} PRIVATE -fopenacc -foffload=nvptx-none)
diff --git a/GEOSagcm_GridComp/GEOShs_GridComp/GEOS_HsGridComp.F90 b/GEOSagcm_GridComp/GEOShs_GridComp/GEOS_HsGridComp.F90
new file mode 100644
index 000000000..281300d5e
--- /dev/null
+++ b/GEOSagcm_GridComp/GEOShs_GridComp/GEOS_HsGridComp.F90
@@ -0,0 +1,1129 @@
+!  $Id: GEOS_HsGridComp.F90,v 1.27 2007/05/17 13:09:31 f4mjs Exp $
+
+#include "MAPL_Generic.h"
+
+  module GEOS_HSGridCompMod
+
+!BOP
+! \renewcommand{\comp}{{\tt GEOS\_HSGridCompMod}}
+!
+! !MODULE: GEOS_HSGridCompMod
+!
+! !DESCRIPTION:
+!
+!  \comp  implements the Held and Suarez (1994) benchmark
+!  forcing for global atmospheric dynamical cores.
+!  The Williamson, et al. (1998) stratospheric modifications
+!  to the benchmark are also implemented, and a mechanism is included to add
+!  a localized, gaussian heat source.
+!
+!
+!  {\it Scientific Basis:} 
+!
+! {\bf The Held-Suarez forcing:}
+!The wind tendencies are given by
+!$$
+!\frac{\partial {\bf v}}{\partial t} = \cdots - k_v(\sigma) \, {\bf v}
+!$$
+!where
+!$$
+! k_v(\sigma) = k_f f_1(\sigma)
+!$$
+!and
+!$$
+! f_1(\sigma) = \max \left[\,0, \frac{\sigma - \sigma_b}{1 - \sigma_b} \right]
+!$$
+!
+! The heating tendencies are given by:
+!$$
+!\frac{\partial T}{\partial t} = \cdots - k_T(\phi,\sigma) \left[\, T - T_{eq}(\phi,p) \,\right]
+!$$
+!$$
+!T_{eq} = \max \left[   T_{strat}, 
+!   \left(\frac{p}{p_\circ}\right)^\kappa 
+!   \left( T_0  - (\Delta T)_y\,\sin^2\!\phi 
+!          - (\Delta\theta)_z  \log\left(\frac{p}{p_\circ}\right) \cos^2\!\phi\, 
+!   \right)     \right] 
+!$$
+!where
+!$$
+!k_T = k_a + (k_s - k_a) \cos^4\!\phi \, f_1(\sigma).
+!$$
+!
+! The parameter values suggested by Held and Suarez (1994) are:
+! $$\sigma_b = .7,$$
+! $$k_f = 1 \hbox{ day}^{-1},$$
+! $$k_a = \frac{1}{40 \hbox{ days}},$$
+! $$k_s = \frac{1}{4 \hbox{ days}},$$
+! $$(\Delta T)_y = 60  \hbox{K },$$
+! $$(\Delta \theta)_z = 10 \hbox{K },$$
+! $$T_{strat} = 200 \hbox{K },$$
+! $$T_0 = 315 \hbox{K }.$$
+!
+! {\bf The Williamson Modification:}
+!  Williamson's modification in the stratosphere simply alters $T_{eq}$ at upper levels.
+!  For $p<p_D$, the constant HS value is replaced with the following distribution:
+!$$
+!T_{eq} = T_{strat} \,\left[
+! \min \left(1,\frac{p}{p_D}\right)^\frac{R_d \gamma_D}{g} +  
+! \min \left(1,\frac{p}{p_I}\right)^\frac{R_d \gamma_I}{g} - 1 \right], 
+!$$
+! where $T_{strat}$ is as in HS,
+!$$
+!  p_I = p_D - \frac{p_D-p(\sigma=0)}{2}( 1 + \tanh(\alpha\frac{|\phi|-\phi_0}{(\delta\phi)_0})),
+!$$
+!and defaults for the other parameters are
+!  $\gamma_I=-3.345 \times 10^{-3}$ K m$^{-1}$, 
+!                $\gamma_D=2.0 \times 10^{-3}$ K m$^{-1}$,
+!$p_D = 100 \hbox{ hPa} $,
+!$\alpha = 2.65 $,
+!$(\delta\phi)_0 = 15^\circ $,
+!$\phi_0 = 60^\circ $.
+!
+! {\bf The Local Heat Source:}
+! A local heat form can also be included. It has the form
+!$$
+!  Q(\lambda,\phi) = Q_{max} e^{ -\frac{1}{2} 
+!    \left( \frac{(\phi_h-\phi)^2}{(\delta\phi)_h^2}
+!   +       \frac{(\lambda_h-\lambda)^2}{(\delta\lambda)_h^2} \right) }
+!     h(p),
+!$$
+!where $\lambda$ is the longitude, $p_s$ is the surface pressure,
+!\begin{equation*}
+!     h(p) = \left\{ \begin{array}{ll}
+!                     0   &  \mbox{ $p < p_1$} \\
+
+!                     \sin(\pi  \frac{p_s-p}{p_s-p_1}) & \mbox{$ p_1 \le p < p_s$}, \\
+!                   \end{array}
+!            \right.
+!\end{equation*}
+!  $Q_{max}$ is in K day$^{-1}$, with a default of zero, and defaults for the
+!  other parameters are $\phi_h=\lambda_h=0$, 
+!  $(\delta\phi)_h=5^\circ$, $(\delta\lambda)_h=30^\circ$,
+!  $p_1=$0 hPa.
+!  
+!
+!
+!  {\it Code Implementation:} 
+!  
+!  
+!   All HS parameters are
+!  optionally adjustable from the configuration, with published values
+!  being the defaults.  Its parameters are also adjustable
+!  from the configuration and also default to published values. The Williamson
+!  modification can be disabled by choosing $p_D=0$. Finally, we also include
+!  a mechanism for adding a localized, gaussian heat source whose strength and shape can
+!  be controlled from the configuration. This is zero by default.
+!
+! {\it References:}
+!
+!  Held, I. M., and M. J. Suarez, 1994: A proposal for the intercomparison
+!  of the dynamical cores of atmospheric general circulation models. {\em Bulletin
+!  of the American Meteorological Society}, {\bf 75(10)}, 1825-1830. 
+!
+!  Williamson, D.L., J. G. Olson, and B. A. Boville, 1998: A comparison of 
+!  semi-Lagrangian and Eulerian tropical climate simulations. {\em Mon. Wea. Rev.}, 
+!  {\bf 126}, 1001-1012
+!
+
+! !USES:
+
+  use ESMF
+  use MAPL_Mod
+  use hs_oacc_mod
+  use iso_c_binding
+  use, intrinsic :: ieee_arithmetic
+  use ieee_exceptions, only: ieee_get_halting_mode, ieee_set_halting_mode, ieee_all
+
+  implicit none
+  private
+!
+! !PUBLIC MEMBER FUNCTIONS:
+
+  public  SetServices
+
+!EOP
+
+   contains
+
+!BOP
+
+! !IROUTINE: SetServices -- Sets ESMF services for this component
+
+! !DESCRIPTION:  This version uses the MAPL\_GenericSetServices, which sets
+!                the Initialize and Finalize services, as well as allocating
+!   our instance of the MAPL\_MetaComp and putting it in the 
+!   gridded component (GC). The MAPL\_MetaComp contains an ESMF state to use as
+!   an internal state. This routine describes the contents of this Internal state,
+!   as well as of the conventional Import and Export states by making call 
+!   to MAPL.
+!
+!
+
+! !INTERFACE:
+
+  subroutine SetServices ( GC, RC )
+
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(INOUT) :: GC  ! gridded component
+    integer, optional,   intent(  OUT) :: RC  ! return code
+
+!EOP
+
+!=============================================================================
+
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)        :: IAm
+    integer                           :: STATUS
+    character(len=ESMF_MAXSTR)        :: COMP_NAME
+
+!=============================================================================
+
+! Begin...
+
+! Get my name and set-up traceback handle
+! ---------------------------------------
+
+    Iam = 'SetServices'
+    call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // trim(Iam)
+
+! Set Data Services for the GC
+! ----------------------------
+
+!BOP
+
+! !IROUTINE: State Descriptions
+
+! !DESCRIPTION: The component uses all three states (Import, Export
+!  and Internal). There is no Private (non-ESMF) Internal state. All
+!  three are managed by MAPL. The Internal state contains only invariant,
+!  horizontally dependent quantities that are set in Initialize, therefore
+!  it never needs checkpointing.
+
+
+! !IMPORT STATE:
+
+    call MAPL_AddImportSpec(GC,                                    &
+         SHORT_NAME = 'U',                                         &
+         LONG_NAME  = 'eastward_wind',                             &
+         UNITS      = 'm s-1',                                     &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddImportSpec(GC,                                    &
+         SHORT_NAME = 'V',                                         &
+         LONG_NAME  = 'northward_wind',                            &
+         UNITS      = 'm s-1',                                     &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddImportSpec(GC,                                    &
+         SHORT_NAME = 'TEMP',                                      &
+         LONG_NAME  = 'air_temperature',                           &
+         UNITS      = 'K',                                         &
+         DEFAULT    = 300.0,                                       &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationCenter,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddImportSpec(GC,                                    &
+         SHORT_NAME = 'PLE',                                       &
+         LONG_NAME  = 'edge_pressure',                             &
+         UNITS      = 'Pa',                                        &
+         DIMS       = MAPL_DimsHorzVert,                           &
+         VLOCATION  = MAPL_VLocationEdge,                          &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+
+! !INTERNAL STATE:
+
+    call MAPL_AddInternalSpec(gc,                                  &
+         SHORT_NAME = 'SPHI2',                                     &
+         LONG_NAME  = 'sine_latitude_squared',                     &
+         UNITS      = '1',                                         &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+         RESTART    = MAPL_RestartSkip,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(gc,                                  &
+         SHORT_NAME = 'CPHI2',                                     &
+         LONG_NAME  = 'cosine_latitude_squared',                   &
+         UNITS      = '1',                                         &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+         RESTART    = MAPL_RestartSkip,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(gc,                                  &
+         SHORT_NAME = 'HFCN',                                      &
+         LONG_NAME  = 'horizontal_structure_of_localized_heating', &
+         UNITS      = '1',                                         &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+         RESTART    = MAPL_RestartSkip,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddInternalSpec(gc,                                  &
+         SHORT_NAME = 'P_I',                                       &
+         LONG_NAME  = 'Williamson_interface_pressure',             &
+         UNITS      = '1',                                         &
+         DEFAULT    = 0.0,                                         &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+         RESTART    = MAPL_RestartSkip,                        &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+
+! !EXPORT STATE:
+
+    call MAPL_AddExportSpec(GC,                                    &
+         SHORT_NAME = 'DTDT',                                      &
+         LONG_NAME  = 'pressure_weighted_air_temperature_tendency',&
+         UNITS      = 'Pa K s-1',                                  &
+         DIMS       =  MAPL_DimsHorzVert,                          &
+         VLOCATION  =  MAPL_VLocationCenter,                       &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+    
+    call MAPL_AddExportSpec(GC,                                    &
+         SHORT_NAME = 'DUDT',                                      &
+         LONG_NAME  = 'eastward_wind_tendency',                    &
+         UNITS      = 'm s-2',                                     &
+         DIMS       =  MAPL_DimsHorzVert,                          &
+         VLOCATION  =  MAPL_VLocationCenter,                       &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    & 
+         SHORT_NAME = 'DVDT',                                      &
+         LONG_NAME  = 'northward_wind_tendency',                   &
+         UNITS      = 'm s-2',                                     &
+         DIMS       =  MAPL_DimsHorzVert,                          &
+         VLOCATION  =  MAPL_VLocationCenter,                       &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    & 
+         SHORT_NAME = 'T_EQ',                                      &
+         LONG_NAME  = 'equilibrium_temperature',                   &
+         UNITS      = 'K',                                         &
+         DIMS       =  MAPL_DimsHorzVert,                          &
+         VLOCATION  =  MAPL_VLocationCenter,                       &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    & 
+         SHORT_NAME = 'THEQ',                                      &
+         LONG_NAME  = 'equilibrium_potential_temperature',         &
+         UNITS      = 'K',                                         &
+         DIMS       =  MAPL_DimsHorzVert,                          &
+         VLOCATION  =  MAPL_VLocationCenter,                       &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    &
+         SHORT_NAME = 'TAUX',                                      &
+         LONG_NAME  = 'eastward_surface_stress',                   &
+         UNITS      = 'N m-2',                                     &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    &
+         SHORT_NAME = 'TAUY',                                      &
+         LONG_NAME  = 'northward_surface_stress',                  &
+         UNITS      = 'N m-2',                                     &
+         DIMS       = MAPL_DimsHorzOnly,                           &
+         VLOCATION  = MAPL_VLocationNone,                          &
+                                                        RC=STATUS  )
+     VERIFY_(STATUS)
+
+    call MAPL_AddExportSpec(GC,                                    & 
+         SHORT_NAME = 'DISS',                                      &
+         LONG_NAME  = 'frictional_dissipation',                    & 
+         UNITS      = 'W m-2',                                     &
+         DIMS       =  MAPL_DimsHorzOnly,                          &
+         VLOCATION  =  MAPL_VLocationNone,                         &
+                                                        RC=STATUS  )
+    VERIFY_(STATUS)
+!EOP
+
+! Set the Initialize and Run entry points
+! --------------------------------------------------
+
+    call MAPL_GridCompSetEntryPoint ( gc, ESMF_METHOD_INITIALIZE, Initialize, rc=status)
+    VERIFY_(STATUS)
+    call MAPL_GridCompSetEntryPoint ( gc, ESMF_METHOD_RUN,  Run,        rc=status)
+    VERIFY_(STATUS)
+
+! Set the Profiling timers
+! ------------------------
+
+    call MAPL_TimerAdd(GC,    name="RUN"   ,RC=STATUS); VERIFY_(STATUS)
+    call MAPL_TimerAdd(GC,    name="INIT"  ,RC=STATUS); VERIFY_(STATUS)
+
+! SetServices clean-up on the way back up through the hierarchy
+!--------------------------------------------------------------
+
+    call MAPL_GenericSetServices    ( GC, RC=STATUS)
+    VERIFY_(STATUS)
+
+! All Done
+!---------
+
+    RETURN_(ESMF_SUCCESS)
+  end subroutine SetServices
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!BOP
+! !IROUTINE: Initialize
+!
+! !DESCRIPTION: Here we initialize the internal state, which in HS contains
+!   two-dimensional invariant arrays used in the forcing calculations. They
+!   are done here simply for economy. Initialize calls MAPL\_GenericInitialize.  
+!   If the import state needs to be restarted, MAPL will do it 
+!   if a restart file is provided. 
+!
+! !INTERFACE:
+
+  subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC )
+
+    use openacc
+    
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_State),    intent(inout) :: IMPORT ! Import state
+    type(ESMF_State),    intent(inout) :: EXPORT ! Export state
+    type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
+    integer, optional,   intent(  out) :: RC     ! Error code
+
+ 
+!EOP
+
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)    :: IAm 
+    integer                       :: STATUS
+    character(len=ESMF_MAXSTR)    :: COMP_NAME
+
+! Pointers to geography info in the MAPL MetaComp
+
+    real,                 pointer :: LATS (:,:)
+    real,                 pointer :: LONS (:,:)
+
+! Pointers to internals
+
+    real,                 pointer :: SPHI2(:,:)
+    real,                 pointer :: CPHI2(:,:)
+    real,                 pointer :: HFCN (:,:)
+    real,                 pointer :: P_I  (:,:)
+
+! Local variables
+
+    real                          :: DX, DY, X0, Y0
+    real                          :: AFAC, PHI0, QMAX
+    logical                       :: halting_mode(5)
+
+! Local derived type aliases
+
+    type (MAPL_MetaComp), pointer :: MAPL
+    type (ESMF_State   )          :: INTERNAL 
+
+!=============================================================================
+
+! Begin... 
+
+! Get the target components name and set-up traceback handle.
+! -----------------------------------------------------------
+
+    Iam = "Initialize"
+    call ESMF_GridCompGet       ( GC, name=COMP_NAME,        RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // Iam
+
+! Get my instance of the MAPL MetaComp
+!-----------------------------------
+
+    call MAPL_GetObjectFromGC   ( GC, MAPL,                   RC=STATUS)
+    VERIFY_(STATUS)
+
+! Call Generic Initialize 
+! -----------------------
+
+    call MAPL_GenericInitialize ( GC, IMPORT, EXPORT, CLOCK,  RC=STATUS)
+    VERIFY_(STATUS)
+
+! Start total timer
+!------------------
+
+    call MAPL_TimerOn(MAPL,"TOTAL")
+    call MAPL_TimerOn(MAPL,"INIT" )
+
+! Get coordinate information from MAPL_MetaComp
+!--------------------------------------------
+
+    call MAPL_Get(MAPL,                &
+       LATS          = LATS,           & ! These are in radians
+       LONS          = LONS,           & ! These are in radians
+       INTERNAL_ESMF_STATE=INTERNAL,   &
+                             RC=STATUS )
+    VERIFY_(STATUS)
+
+! Get pointers to internal variables
+!-----------------------------------
+
+    call MAPL_GetPointer(INTERNAL, SPHI2, 'SPHI2', RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, CPHI2, 'CPHI2', RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, HFCN , 'HFCN' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, P_I  , 'P_I'  , RC=STATUS); VERIFY_(STATUS)
+
+! Initialize geometric factors
+!-----------------------------
+
+    SPHI2 = sin(LATS)**2
+    CPHI2 = cos(LATS)**2
+
+! Precompute Local heating distribution
+!--------------------------------------
+
+
+!BOR
+
+! !RESOURCE_ITEM: K day-1 :: Local heating 
+    call MAPL_GetResource(MAPL,QMAX,'QMAX:',default=.000    , RC=STATUS ) 
+    VERIFY_(STATUS)
+
+    if(QMAX /= 0.0) then
+
+! Get local heating parameters from the configuration
+!----------------------------------------------------
+
+! !RESOURCE_ITEM: degrees :: Central longitude of local heating, $\lambda_h$
+       call MAPL_GetResource(MAPL,X0  ,'X0:'  ,default=0.0     , RC=STATUS)
+       VERIFY_(STATUS)                                               
+! !RESOURCE_ITEM: degrees :: Central latitude of local heating, $\phi_h$
+       call MAPL_GetResource(MAPL,Y0  ,'Y0:'  ,default=0.0     , RC=STATUS)
+       VERIFY_(STATUS)                                               
+! !RESOURCE_ITEM: degrees :: Longitudinal width of local heating, $(\delta\lambda)_h$
+       call MAPL_GetResource(MAPL,DX  ,'DX:'  ,default=30.     , RC=STATUS)
+       VERIFY_(STATUS)                                               
+! !RESOURCE_ITEM: degrees :: Latitudinal width of local heating, $(\delta\phi)_h$
+       call MAPL_GetResource(MAPL,DY  ,'DY:'  ,default=5.0     , RC=STATUS)
+       VERIFY_(STATUS)                             
+
+       ! Horizontal structure of optional stationary heat source
+       !--------------------------------------------------------
+
+       X0 = X0 * (MAPL_PI/180.0)
+       Y0 = Y0 * (MAPL_PI/180.0)
+       DX = DX * (MAPL_PI/180.0)
+       DY = DY * (MAPL_PI/180.0)
+
+       HFCN =          ((LONS     -X0)/DX)**2
+       HFCN = min(HFCN,((LONS+360.-X0)/DX)**2)
+       HFCN = min(HFCN,((LONS-360.-X0)/DX)**2)
+
+       HFCN = HFCN +   ((LATS     -Y0)/DY)**2
+
+       where(HFCN < 10.)
+          HFCN = exp(-0.5*HFCN)
+       elsewhere
+          HFCN = 0.0
+       end where
+    else
+       HFCN = 0.0
+    endif
+                        
+! Get Williamson parameters from the configuration
+!-------------------------------------------------
+                                                
+! !RESOURCE_ITEM: 1 :: Williamson's $\alpha$ parameter
+    call MAPL_GetResource(MAPL,AFAC,'AFAC:',default=2.65/15., RC=STATUS)
+    VERIFY_(STATUS)
+! !RESOURCE_ITEM: degrees :: Williamson's $\phi_o$ parameter
+    call MAPL_GetResource(MAPL,PHI0,'PHI0:',default=60.     , RC=STATUS)
+    VERIFY_(STATUS)                                      
+!EOR
+                                           
+! Horizontal structure of Willianson parameterization
+!----------------------------------------------------
+
+    P_I = ( 1.0 + tanh( AFAC*( abs(LATS) - PHI0 ) ) )
+
+! All Done
+!---------
+
+    call MAPL_TimerOff(MAPL,"INIT" )
+    call MAPL_TimerOff(MAPL,"TOTAL")
+
+    !OACC GCC10 Implementation can trigger floating-point exception
+    !per https://gcc.gnu.org/wiki/OpenACC/Implementation%20Status
+    !we saw those on PTX load so we deactivate the ieee mode around
+    !the init
+    call ieee_get_halting_mode(ieee_all, halting_mode)
+    call ieee_set_halting_mode(ieee_all, .false.)
+    call acc_set_device_type(acc_device_default)
+    call ieee_set_halting_mode(ieee_all, halting_mode)
+
+    RETURN_(ESMF_SUCCESS)
+  end subroutine INITIALIZE
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+!BOP
+
+! !IROUTINE: Run
+
+! !DESCRIPTION: The Run method of the HS Gridded Component. It computes
+!  tendencies for temperature and winds as described above,
+!  as well as a number of
+!  diagnostics. These are all optional Exports. The input winds 
+!  and/or temperatures may be Friendly. All Friendlies are updated
+!  using a time step
+!  obtained from 'RUN\_DT:' in the configuration. Both U and V must be
+!  Friendly for the wind to be updated. It is an error for one to be Friendly
+!  and not the other.  
+!
+! !INTERFACE:
+
+  subroutine Run( GC, IMPORT, EXPORT, CLOCK, RC )
+    include 'mpif.h'
+
+    INTERFACE
+        subroutine c_call_hs_oacc(CPHI2, DISS, DTDT, DUDT, DVDT, &
+                          	HFCN, P_I, PLE, SPHI2, TAUX, TAUY, T, &
+                          	THEQ, T_EQ, U, V, &
+                          	DAYLEN, DELH, DELV1, DT, FRICQ, FriendlyTemp, &
+                          	FriendlyWind, GAM_D, GAM_I, IM, JM, LM, P_1, P_D, QMAX, &
+                          	SIG1, TAUA, TAUF, TAUS, TSTRT, T0, compType, rank) bind(C, name='c_call_hs_oacc')
+        
+            use, intrinsic :: iso_c_binding
+
+            integer(C_INT),  value :: FRICQ, IM, JM, LM, compType, rank
+            ! The logical variables are declared with (C_INT) since (C_BOOL) seems
+            ! to be 1 byte, while Fortran (with gfortran) uses 4-byte logicals
+            logical(C_INT), value :: FriendlyTemp, FriendlyWind
+            real(C_FLOAT),   value :: DAYLEN, DELH, DELV1, DT, GAM_D, GAM_I, P_D, P_1
+            real(C_FLOAT),   value :: QMAX, SIG1, T0, TAUA, TAUF, TAUS, TSTRT
+
+            type(C_PTR), value :: CPHI2, HFCN, P_I, SPHI2
+            type(C_PTR), value :: DISS, TAUX, TAUY
+            type(C_PTR), value :: PLE
+            type(C_PTR), value :: DTDT, DUDT, DVDT, T, T_EQ, THEQ, U, V
+        end subroutine
+    end interface
+
+! !ARGUMENTS:
+
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_State),    intent(inout) :: IMPORT ! Import state
+    type(ESMF_State),    intent(inout) :: EXPORT ! Export state
+    type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
+    integer, optional,   intent(  out) :: RC     ! Error code
+
+ 
+!EOP
+
+! ErrLog Variables
+
+    character(len=ESMF_MAXSTR)        :: IAm 
+    integer                           :: STATUS
+    character(len=ESMF_MAXSTR)        :: COMP_NAME
+    
+! Local derived type aliases
+
+    type (MAPL_MetaComp    ), pointer :: MAPL
+    type (ESMF_State       )          :: INTERNAL
+    type (ESMF_Grid        )          :: GRID
+    type (ESMF_Logical     )          :: Friendly
+
+    type(ESMF_VM) :: vm
+
+    integer :: comm, rank, ierr, length
+    character * (MPI_MAX_PROCESSOR_NAME) nameNode
+
+! Pointers to inputs
+
+    real, pointer, dimension(:,:,:)   :: PLE, T, U, V
+
+! Pointers to internals
+
+    real, pointer, dimension(:,:  )   :: HFCN, SPHI2, CPHI2, P_I
+
+! Pointers to outputs
+
+    real, pointer, dimension(:,:,:)   :: DTDT, DUDT, DVDT, DTDT_p, DUDT_p, DVDT_p
+    real, pointer, dimension(:,:,:)   :: THEQ, T_EQ, THEQ_p, T_EQ_p
+    real, pointer, dimension(:,:  )   :: TAUX, TAUY, DISS, TAUX_p, TAUY_p, DISS_p
+
+! Scratch arrays and working pointers
+
+    real, allocatable, dimension(:,:) :: PII
+    real, allocatable, dimension(:,:) :: DP
+    real, allocatable, dimension(:,:) :: PL
+    real, allocatable, dimension(:,:) :: UU
+    real, allocatable, dimension(:,:) :: VV
+    real, allocatable, dimension(:,:) :: VR
+    real, allocatable, dimension(:,:) :: TE
+    real, allocatable, dimension(:,:) :: F1
+    real, allocatable, dimension(:,:) :: RR
+    real, allocatable, dimension(:,:) :: DS
+    real, allocatable, dimension(:,:) :: DM
+    real, allocatable, dimension(:,:) :: PK
+
+    real,     pointer, dimension(:,:) :: PS
+    real,     pointer, dimension(:,:) :: PT
+
+    integer :: L
+    integer :: IM, JM, LM
+    integer :: FRICQ
+
+    logical :: FriendlyTemp
+    logical :: FriendlyWind
+
+    real    :: DT
+    real    :: KA  
+    real    :: KS  
+    real    :: KF
+
+! 8 Held-Suarez parameters
+
+    real    :: SIG1
+    real    :: TAUA
+    real    :: TAUS
+    real    :: TAUF
+    real    :: DELV1
+    real    :: DELH
+    real    :: T0
+    real    :: TSTRT
+
+! 3 Localized heating parameters
+
+    real    :: P_1
+    real    :: QMAX
+
+! 3 Williamson parameters
+
+    real    :: P_D
+    real    :: GAM_I
+    real    :: GAM_D
+    
+    real, parameter :: DAYLEN=86400.
+
+    integer :: run_gtfv3 = 0
+!=============================================================================
+
+! Begin... 
+
+! Get the target components name and set-up traceback handle.
+! -----------------------------------------------------------
+
+    Iam = "Run"
+    call ESMF_GridCompGet( GC, name=COMP_NAME, RC=STATUS )
+    VERIFY_(STATUS)
+    Iam = trim(COMP_NAME) // trim(Iam)
+
+! Get my MAPL_MetaComp
+!---------------------
+
+    call MAPL_GetObjectFromGC( GC, MAPL, RC=STATUS)
+    VERIFY_(STATUS)
+
+! Start the TOTAL timer
+!----------------------
+
+    call MAPL_TimerOn(MAPL,"TOTAL")
+    call MAPL_TimerOn(MAPL,"RUN")
+
+! Get esmf internal state from generic state.
+!-------------------------------------------
+
+    call MAPL_Get(MAPL,                  &
+         IM = IM, JM=JM, LM=LM,          &
+         INTERNAL_ESMF_STATE=INTERNAL,   &
+                               RC=STATUS )
+    VERIFY_(STATUS)
+
+! Pointers TO INTERNALS
+!----------------------
+
+    call MAPL_GetPointer(INTERNAL, SPHI2, 'SPHI2', RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, CPHI2, 'CPHI2', RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, HFCN , 'HFCN' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(INTERNAL, P_I  , 'P_I'  , RC=STATUS); VERIFY_(STATUS)
+
+! Pointers TO IMPORTS
+!-------------------- 
+
+    call MAPL_GetPointer(IMPORT  , PLE  , 'PLE'  , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(IMPORT  , T    , 'TEMP' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(IMPORT  , U    , 'U'    , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(IMPORT  , V    , 'V'    , RC=STATUS); VERIFY_(STATUS)
+
+! Pointers TO EXPORTS
+!-------------------- 
+
+    call MAPL_GetPointer(EXPORT  , DUDT , 'DUDT' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , DVDT , 'DVDT' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , DTDT , 'DTDT' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , DISS , 'DISS' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , THEQ , 'THEQ' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , T_EQ , 'T_EQ' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , TAUX , 'TAUX' , RC=STATUS); VERIFY_(STATUS)
+    call MAPL_GetPointer(EXPORT  , TAUY , 'TAUY' , RC=STATUS); VERIFY_(STATUS)
+
+! Get parameters from the configuration
+!--------------------------------------
+
+! 8 Held-Suarez parameters
+
+!BOR
+
+! !RESOURCE_ITEM: days :: H-S $k_a^{-1}$ parameter
+    call MAPL_GetResource(MAPL,TAUA ,'TAUA:' ,default=40. ,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: days :: H-S $k_s^{-1}$ parameter
+    call MAPL_GetResource(MAPL,TAUS ,'TAUS:' ,default=4.0 ,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: days :: H-S $k_f^{-1}$ parameter
+    call MAPL_GetResource(MAPL,TAUF ,'TAUF:' ,default=1.0 ,     RC=STATUS )
+    VERIFY_(STATUS)
+! !RESOURCE_ITEM: days :: H-S $\sigma_b$ parameter
+    call MAPL_GetResource(MAPL,SIG1 ,'SIG1:' ,default=0.7 ,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: K :: H-S $T_o$ parameter
+    call MAPL_GetResource(MAPL,T0   ,'T0:'   ,default=315.,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: K :: H-S $(\Delta \theta)_z$ parameter
+    call MAPL_GetResource(MAPL,DELV1,'DELV1:',default=10. ,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: K :: H-S $(\Delta T)_y$ parameter
+    call MAPL_GetResource(MAPL,DELH ,'DELH:' ,default=60. ,     RC=STATUS )
+    VERIFY_(STATUS)                  
+! !RESOURCE_ITEM: K :: H-S $T_{strat}$ parameter
+    call MAPL_GetResource(MAPL,TSTRT,'TSTRT:',default=200.,     RC=STATUS )
+    VERIFY_(STATUS)
+! !RESOURCE_ITEM: 1 OR 0 :: Controls addition of dissipation heat.
+    call MAPL_GetResource(MAPL,FRICQ,'DISSQ:',default=1,        RC=STATUS )
+    VERIFY_(STATUS)
+
+! 3 Localized heating parameters
+
+! !RESOURCE_ITEM: Pa :: Local heating $p_1$ parameter
+    call MAPL_GetResource(mapl,P_1  ,'P1:'   ,default=0.0 ,     RC=STATUS )
+    VERIFY_(STATUS)
+    call MAPL_GetResource(mapl,QMAX ,'QMAX:' ,default=0.0 ,     RC=STATUS )
+    VERIFY_(STATUS)                                        
+
+! 3 Williamson parameters
+
+! !RESOURCE_ITEM: K m-1 :: Williamson's $\gamma_I$ parameter
+    call MAPL_GetResource(mapl,GAM_I,'GAM_I:',default=-3.345e-3,RC=STATUS)
+    VERIFY_(STATUS)                                        
+! !RESOURCE_ITEM: K m-1 :: Williamson's $\gamma_D$ parameter
+    call MAPL_GetResource(mapl,GAM_D,'GAM_D:',default=.002     ,RC=STATUS)
+    VERIFY_(STATUS)                                        
+! !RESOURCE_ITEM: Pa :: Williamson's $p_D$ parameter
+    call MAPL_GetResource(mapl,P_D  ,'P_D:'  ,default=1.E4     ,RC=STATUS)
+    VERIFY_(STATUS)
+
+    call MAPL_GetResource(MAPL, run_gtfv3, 'RUN_GTFV3:', default=0, RC=STATUS)
+    VERIFY_(STATUS)
+  
+!EOR
+
+! Check for Friendliness
+!-----------------------
+
+    Friendly = MAPL_VerifyFriendly(IMPORT,'V',trim(COMP_NAME),RC=STATUS)
+    VERIFY_(STATUS)
+
+    FriendlyWind = Friendly == ESMF_TRUE
+
+    Friendly = MAPL_VerifyFriendly(IMPORT,'U',trim(COMP_NAME),RC=STATUS)
+    VERIFY_(STATUS)
+
+    ASSERT_((Friendly == ESMF_TRUE) .eqv. FriendlyWind)
+
+    Friendly = MAPL_VerifyFriendly(IMPORT,'TEMP',trim(COMP_NAME),RC=STATUS)
+    VERIFY_(STATUS)
+
+    FriendlyTemp = Friendly == ESMF_TRUE
+
+! If we have a friendly, we need a time step
+!-------------------------------------------
+    
+    if(FriendlyTemp .or. FriendlyWind) then
+       call MAPL_GetResource(MAPL, DT,  'RUN_DT', RC=STATUS)
+       VERIFY_(STATUS)
+    else
+       DT = huge(0.0)
+    end if
+
+! Allocate 10 2D scratch arrays
+!------------------------------
+
+    allocate(DP (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(PL (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(UU (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(VV (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(VR (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(TE (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(DS (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(PII(IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(F1 (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(RR (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(DM (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+    allocate(PK (IM,JM),STAT=STATUS); VERIFY_(STATUS)
+
+! Get MPI Rank
+!-------------
+
+  call ESMF_VMGetCurrent(vm, RC=status)
+  call ESMF_VMGet(vm, mpiCommunicator=comm, RC=status)
+  call MPI_Comm_rank(comm, rank, ierr)
+
+  if (run_gtfv3 == 0) then
+    ! Begin calculations.
+    !-------------------
+
+        PS  => PLE(:,:,LM)
+        PT  => PLE(:,:, 0)
+
+        PII =  P_D - (P_D - PT)*0.5*P_I
+
+    ! Initialize vertically integrated diagnostics
+    !---------------------------------------------
+
+        if(associated(DISS)) DISS = 0.0
+        if(associated(TAUX)) TAUX = 0.0
+        if(associated(TAUY)) TAUY = 0.0
+
+    ! Loop invariants
+    !----------------
+
+        KA   = 1.0/(DAYLEN*TAUA)
+        KS   = 1.0/(DAYLEN*TAUS)
+        KF   = 1.0/(DAYLEN*TAUF)
+
+        LEVELS: do L = 1,LM
+
+        DP  = (PLE(:,:,L)-PLE(:,:,L-1))
+        PL  = (PLE(:,:,L)+PLE(:,:,L-1))*0.5
+        DM  = DP / MAPL_GRAV
+        PK  = (PL/MAPL_P00)**MAPL_KAPPA 
+
+    ! H&S equilibrium temperature
+    !----------------------------
+
+        TE  = PK*( T0 - DELH*SPHI2 - DELV1*CPHI2*log( PL/MAPL_P00 ) )
+        TE  = max( TE, TSTRT )
+
+    ! Williamson Stratospheric modifications to equilibrium temperature
+    ! -----------------------------------------------------------------
+
+        where( PL < P_D )
+            TE  = TSTRT*( min(1.0,PL/P_D)**(MAPL_RGAS*GAM_D/MAPL_GRAV)     &
+                        + min(1.0,PL/PII)**(MAPL_RGAS*GAM_I/MAPL_GRAV) - 1 )
+        end where
+
+    !  Exports of equilibrium T and Theta
+    !------------------------------------
+
+        if(associated(T_EQ)) T_EQ(:,:,L) = TE
+        if(associated(THEQ)) THEQ(:,:,L) = TE/PK
+
+    ! Vertical structure of timescales in H&S.
+    !---------------------------------------------
+
+        F1  = max(0.0, ( (PL/PS)-SIG1 )/( 1.0-SIG1 ) )
+
+    ! Atmospheric heating from H&S
+    !-----------------------------
+
+        RR = (KA + (KS-KA)*F1*CPHI2**2) * (TE-T(:,:,L))
+
+        if(associated(DTDT)) DTDT(:,:,L) =            DP*RR
+        if(FriendlyTemp    ) T   (:,:,L) = T(:,:,L) + DT*RR
+
+    ! Wind tendencies
+    !----------------
+
+        UU  = -U(:,:,L)*(F1*KF)
+        VV  = -V(:,:,L)*(F1*KF)
+
+        if(associated(DUDT)) DUDT(:,:,L) = UU
+        if(associated(DVDT)) DVDT(:,:,L) = VV
+
+        if(FriendlyWind) then
+            U(:,:,L) = U(:,:,L) + DT*UU
+            V(:,:,L) = V(:,:,L) + DT*VV
+        end if
+
+    !  Frictional heating from H&S drag
+    !----------------------------------
+
+        DS = U(:,:,L)*UU + V(:,:,L)*VV
+
+        if(associated(DISS)) DISS           = DISS        - DS*DM
+        if(FRICQ /= 0) then
+            if(associated(DTDT)) DTDT(:,:,L) = DTDT(:,:,L) - DS*(DP/MAPL_CP  )
+            if(FriendlYTemp    ) T   (:,:,L) = T   (:,:,L) - DS*(DT/MAPL_CP  )
+        end if
+
+    !  Surface stresses from vertically integrated H&S surface drag
+    !--------------------------------------------------------------
+
+        if(associated(TAUX)) TAUX = TAUX - UU*DM
+        if(associated(TAUY)) TAUY = TAUY - VV*DM
+
+    ! Localized heat source, if any
+    !------------------------------
+
+        if((associated(DTDT).or.FriendlyTemp) .and. QMAX/=0.0) then
+            where(PL > P_1)
+                VR = HFCN*(QMAX/DAYLEN)*sin( MAPL_PI*(PS-PL)/(PS-P_1) )
+            elsewhere
+                VR = 0.
+            end where
+
+            if(associated(DTDT)) DTDT(:,:,L) = DTDT(:,:,L) + DP*VR 
+            if(FriendlyTemp    ) T   (:,:,L) = T   (:,:,L) + DT*VR
+        end if
+
+        enddo LEVELS
+
+        ! if(associated(DISS) .eqv. .FALSE.) then
+        !     write(*,*) 'DISS is not associated'
+        ! else
+        !     DISS_p => DISS
+        ! endif
+
+        ! if(associated(TAUX) .eqv. .FALSE.) then
+        !     write(*,*) 'TAUX is not associated'
+        ! else
+        !     TAUX_p => TAUX
+        ! endif
+
+        ! if(associated(TAUY) .eqv. .FALSE.) then
+        !     write(*,*) 'TAUY is not associated'
+        ! else
+        !     TAUY_p => TAUY
+        ! endif
+
+        ! if(associated(DTDT) .eqv. .FALSE.) then
+        !     write(*,*) 'DTDT is not associated'
+        ! else
+        !     DTDT_p => DTDT
+        ! endif
+
+        ! if(associated(DUDT) .eqv. .FALSE.) then
+        !     write(*,*) 'DUDT is not associated'
+        ! else
+        !     DUDT_p => DUDT
+        ! endif
+
+        ! if(associated(DVDT) .eqv. .FALSE.) then
+        !     write(*,*) 'DVDT is not associated'
+        ! else
+        !     DVDT_p => DVDT
+        ! endif
+
+        ! if(associated(THEQ) .eqv. .FALSE.) then
+        !     write(*,*) 'THEQ is not associated'
+        ! else
+        !     THEQ_p => THEQ
+        ! endif
+
+        ! if(associated(T_EQ) .eqv. .FALSE.) then
+        !     write(*,*) 'T_EQ is not associated'
+        ! else
+        !     T_EQ_p => T_EQ
+        ! endif
+
+        ! if(associated(CPHI2) .eqv. .FALSE.) write(*,*) 'CPHI2 is not associated'
+        ! if(associated(HFCN) .eqv. .FALSE.) write(*,*) 'HFCN is not associated'
+        ! if(associated(P_I) .eqv. .FALSE.) write(*,*) 'P_I is not associated'
+        ! if(associated(PLE) .eqv. .FALSE.) write(*,*) 'PLE is not associated'
+        ! if(associated(SPHI2) .eqv. .FALSE.) write(*,*) 'SPHI2 is not associated'
+        ! if(associated(U) .eqv. .FALSE.) write(*,*) 'U is not associated'
+        ! if(associated(V) .eqv. .FALSE.) write(*,*) 'V is not associated'
+        ! if(associated(T) .eqv. .FALSE.) write(*,*) 'T is not associated'
+
+        ! write(*,*) 'Calling c_call_hs_oacc'
+        ! call c_call_hs_oacc(C_LOC(CPHI2), C_LOC(DISS_p), C_LOC(DTDT_p), C_LOC(DUDT_p), C_LOC(DVDT_p), &
+        !                     C_LOC(HFCN), C_LOC(P_I), C_LOC(PLE), C_LOC(SPHI2), C_LOC(TAUX_p), C_LOC(TAUY_p), C_LOC(T), &
+        !                     C_LOC(THEQ_p), C_LOC(T_EQ), C_LOC(U), C_LOC(V), &
+        !                     DAYLEN, DELH, DELV1, DT, FRICQ, FriendlyTemp, &
+        !                     FriendlyWind, GAM_D, GAM_I, IM, JM, LM, P_1, P_D, QMAX, &
+        !                     SIG1, TAUA, TAUF, TAUS, TSTRT, T0, 0)
+    else
+        call MPI_GET_PROCESSOR_NAME(nameNode, length, ierr)
+        !write(*,*) 'Rank = ', rank, 'on node', nameNode
+        DISS_P => null()
+        TAUX_P => null()
+        TAUY_P => null()
+        THEQ_P => null()
+        call held_suarez_oacc(CPHI2, DISS_P, DTDT, DUDT, DVDT, &
+                            HFCN, P_I, PLE, SPHI2, TAUX_P, TAUY_P, T, &
+                            THEQ_P, T_EQ, U, V, &
+                            DAYLEN, DELH, DELV1, DT, FRICQ, FriendlyTemp, &
+                            FriendlyWind, GAM_D, GAM_I, IM, JM, LM, P_1, P_D, QMAX, &
+                            SIG1, TAUA, TAUF, TAUS, TSTRT, T0, 2, rank)
+    end if
+! Free 10 scratch arrays
+!-----------------------
+
+    deallocate(PK ) 
+    deallocate(DM ) 
+    deallocate(RR ) 
+    deallocate(F1 ) 
+    deallocate(PII) 
+    deallocate(DS ) 
+    deallocate(TE ) 
+    deallocate(VR ) 
+    deallocate(VV ) 
+    deallocate(UU )
+    deallocate(PL ) 
+    deallocate(DP ) 
+
+! Close timers
+!-------------
+
+    call MAPL_TimerOff(MAPL,"RUN")
+    call MAPL_TimerOff(MAPL,"TOTAL")
+
+! All Done
+!---------
+
+    RETURN_(ESMF_SUCCESS)
+  end subroutine RUN
+
+end module GEOS_HSGridCompMod
+
+
+
diff --git a/GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc b/GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc
new file mode 100755
index 000000000..758c67ce7
--- /dev/null
+++ b/GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc
@@ -0,0 +1,2 @@
+P_D: 0.0
+
diff --git a/GEOSagcm_GridComp/GEOShs_GridComp/held_suarez_oacc.F90 b/GEOSagcm_GridComp/GEOShs_GridComp/held_suarez_oacc.F90
new file mode 100644
index 000000000..641ac9e0b
--- /dev/null
+++ b/GEOSagcm_GridComp/GEOShs_GridComp/held_suarez_oacc.F90
@@ -0,0 +1,238 @@
+module hs_oacc_mod
+
+  use MAPL_ConstantsMod, only: MAPL_CP, MAPL_GRAV, MAPL_KAPPA, MAPL_P00, MAPL_PI, MAPL_RGAS
+
+contains
+
+  subroutine held_suarez_oacc( &
+       CPHI2, DISS, DTDT, DUDT, DVDT, &
+       HFCN, P_I, PLE, SPHI2, TAUX, TAUY, T, &
+       THEQ, T_EQ, U, V, &
+       DAYLEN, DELH, DELV1, DT, FRICQ, FriendlyTemp, &
+       FriendlyWind, GAM_D, GAM_I, IM, JM, LM, P_1, P_D, QMAX, &
+       SIG1, TAUA, TAUF, TAUS, TSTRT, T0, compType, rank)
+
+    use openacc
+
+    implicit none
+
+    integer, intent(in) :: FRICQ, IM, JM, LM, compType, rank
+    logical, intent(in) :: FriendlyTemp, FriendlyWind
+
+    real, dimension(:,:),   pointer,  intent(in)    :: CPHI2, HFCN, P_I, SPHI2
+    real, dimension(:,:),   pointer,  intent(inout) :: DISS, TAUX, TAUY
+    real, dimension(:,:,:), pointer,  intent(in)	:: PLE
+    real, dimension(:,:,:), pointer,  intent(inout) :: DTDT, DUDT, DVDT, T, T_EQ, THEQ, U, V
+
+    real, intent(in) :: DAYLEN, DELH, DELV1, DT, GAM_D, GAM_I, P_D, P_1
+    real, intent(in) :: QMAX, SIG1, T0, TAUA, TAUF, TAUS, TSTRT
+
+    real, dimension(IM,JM) :: DP, PL, UU, VV, VR, TE, DS, PII, F1, RR, DM, PK
+
+    real :: DP_s, PL_s, UU_s, VV_s, VR_s, TE_s, DS_s, PII_s, F1_s, RR_s, DM_s, PK_s
+
+    real, pointer, dimension(:,:) :: PS
+    real, pointer, dimension(:,:) :: PT
+
+    real :: t1, t2
+
+    integer ::  L, I, J
+    real    :: KA, KF, KS
+
+    integer :: ngpus
+
+    real :: CP, GRAV, KAPPA, P00, PI, RGAS
+
+    CP = MAPL_CP
+    GRAV = MAPL_GRAV
+    KAPPA = MAPL_KAPPA
+    P00 = MAPL_P00
+    PI = MAPL_PI
+    RGAS = MAPL_RGAS
+
+    ! Initialize
+    DP_s = 0.0
+    PL_s = 0.0
+    UU_s = 0.0
+    VV_s = 0.0
+    VR_s = 0.0
+    TE_s = 0.0
+    DS_s = 0.0
+    PII_s = 0.0
+    F1_s= 0.0
+    RR_s = 0.0
+    DM_s = 0.0
+    PK_s = 0.0
+    t1 = 0.0
+    t2 = 0.0
+
+    DP = 0.0
+    PL = 0.0
+    UU = 0.0
+    VV = 0.0
+    VR = 0.0
+    TE = 0.0
+    DS = 0.0
+    PII = 0.0
+    F1= 0.0
+    RR = 0.0
+    DM = 0.0
+    PK = 0.0
+
+    ngpus = acc_get_num_devices(acc_device_default)
+
+    write(*,*) 'HS: rank =', rank, ', comp type: ', compType, ', Number of GPUS:', ngpus
+    ! call acc_set_device_num(mod(rank,4),acc_device_nvidia)
+
+    !write(*,*) 'From HS routine, rank =', rank, ": Number of GPUS = ", ngpus, ': This process is using', acc_get_device_num(acc_device_nvidia)
+
+    ! If running with OpenMP, use these pointer assignments for PS and PT
+    PS  => PLE(:,:,LM)
+    PT  => PLE(:,:, 0)
+
+    ! If running with OpenACC, use these pointer assignments for PS and PT
+    !PS => PLE(:,:,LM+1)
+    !PT => PLE(:,:,1)
+
+    !if(associated(DISS)) DISS = 0.0
+    !if(associated(TAUX)) TAUX = 0.0
+    !if(associated(TAUY)) TAUY = 0.0
+
+    KA   = 1.0/(DAYLEN*TAUA)
+    KS   = 1.0/(DAYLEN*TAUS)
+    KF   = 1.0/(DAYLEN*TAUF)
+
+    PII =  P_D - (P_D - PT)*0.5*P_I
+
+    !$acc data copyin(CPHI2, HFCN, P_I, SPHI2, PLE) &
+    !$acc      copy(DISS, TAUX, TAUY) &
+    !$acc      copy(DTDT, DUDT, DVDT, T, T_EQ, THEQ, U, V) &
+    !$acc      copyin(PS, PT)
+
+    ! write(*,*) 'Running HS Code with outer loop OpenMP/OpenACC Parallelism'
+    call cpu_time(t1)
+
+    !$acc parallel loop collapse(3) present(CPHI2, HFCN, P_I, SPHI2, PLE) &
+    !$acc                           present(DISS, TAUX, TAUY) &
+    !$acc                           present(DTDT, DUDT, DVDT, T, T_EQ, THEQ, U, V) &
+    !$acc                           present(PS, PT) &
+    !$acc                           private(DP_s, PL_s, UU_s, VV_s, VR_s, TE_s, DS_s, PII_s, F1_s, RR_s, DM_s, PK_s)
+    do L = 1,LM
+       do J = 1, JM
+          do I = 1, IM
+
+             ! If running with OpenMP, use these assignments for DP and PL
+             DP_s  = (PLE(I,J,L)-PLE(I,J,L-1))
+             PL_s  = (PLE(I,J,L)+PLE(I,J,L-1))*0.5
+
+             ! If running with OpenACC, use these assignments for DP and PL
+             !DP_s  = (PLE(I,J,L+1)-PLE(I,J,L))
+             !PL_s  = (PLE(I,J,L+1)+PLE(I,J,L))*0.5
+             DM_s  = DP_s / GRAV
+             PK_s  = (PL_s/P00)**KAPPA
+
+             ! H&S equilibrium temperature
+             !----------------------------
+
+             TE_s  = PK_s*( T0 - DELH*SPHI2(I,J) - DELV1*CPHI2(I,J)*log( PL_s/P00 ) )
+             TE_s  = max( TE_s, TSTRT )
+
+             ! Williamson Stratospheric modifications to equilibrium temperature
+             ! -----------------------------------------------------------------
+
+             if( PL_s < P_D ) then
+                TE_s  = TSTRT*( min(1.0,PL_s/P_D)**(RGAS*GAM_D/GRAV)     &
+                     + min(1.0,PL_s/PII_s)**(RGAS*GAM_I/GRAV) - 1 )
+             endif
+
+             !  Exports of equilibrium T and Theta
+             !------------------------------------
+
+             if(associated(T_EQ)) T_EQ(I,J,L) = TE_s
+             ! T_EQ(I,J,L) = TE_s
+             if(associated(THEQ)) THEQ(I,J,L) = TE_s/PK_s
+             ! THEQ(I,J,L) = TE_s/PK_s
+
+             ! Vertical structure of timescales in H&S.
+             !---------------------------------------------
+
+             F1_s  = max(0.0, ( (PL_s/PS(I,J))-SIG1 )/( 1.0-SIG1 ) )
+
+             ! Atmospheric heating from H&S
+             !-----------------------------
+
+             RR_s = (KA + (KS-KA)*F1_s*CPHI2(I,J)**2) * (TE_s-T(I,J,L))
+
+             if(associated(DTDT)) DTDT(I,J,L) = DP_s*RR_s
+             ! DTDT(I,J,L) = DP_s*RR_s
+             if(FriendlyTemp    ) T   (I,J,L) = T(I,J,L) + DT*RR_s
+
+             ! Wind tendencies
+             !----------------
+
+             UU_s  = -U(I,J,L)*(F1_s*KF)
+             VV_s  = -V(I,J,L)*(F1_s*KF)
+
+             if(associated(DUDT)) DUDT(I,J,L) = UU_s
+             ! DUDT(I,J,L) = UU_s
+             if(associated(DVDT)) DVDT(I,J,L) = VV_s
+             ! DVDT(I,J,L) = VV_s
+
+             if(FriendlyWind) then
+                U(I,J,L) = U(I,J,L) + DT*UU_s
+                V(I,J,L) = V(I,J,L) + DT*VV_s
+             end if
+
+             !  Frictional heating from H&S drag
+             !----------------------------------
+
+             DS_s = U(I,J,L)*UU_s + V(I,J,L)*VV_s
+
+             if(associated(DISS)) DISS(I,J) = DISS(I,J) - DS_s*DM_s
+             ! DISS(I,J) = DISS(I,J) - DS_s*DM_s
+             if(FRICQ /= 0) then
+                if(associated(DTDT)) DTDT(I,J,L) = DTDT(I,J,L) - DS_s*(DP_s/CP  )
+                ! DTDT(I,J,L) = DTDT(I,J,L) - DS_s*(DP_s/CP  )
+                if(FriendlYTemp    ) T   (I,J,L) = T   (I,J,L) - DS_s*(DT/CP  )
+             end if
+
+             !  Surface stresses from vertically integrated H&S surface drag
+             !--------------------------------------------------------------
+
+             if(associated(TAUX)) TAUX(I,J) = TAUX(I,J) - UU_s*DM_s
+             ! TAUX(I,J) = TAUX(I,J) - UU_s*DM_s
+             if(associated(TAUY)) TAUY(I,J) = TAUY(I,J) - VV_s*DM_s
+             ! TAUY(I,J) = TAUY(I,J) - VV_s*DM_s
+
+             ! Localized heat source, if any
+             !------------------------------
+
+             if((associated(DTDT).or.FriendlyTemp) .and. QMAX/=0.0) then
+                ! if(FriendlyTemp .and. QMAX/=0.0) then
+                if(PL_s > P_1) then
+                   VR_s = HFCN(I,J)*(QMAX/DAYLEN)*sin( PI*(PS(I,J)-PL_s)/(PS(I,J)-P_1) )
+                else
+                   VR_s = 0.
+                endif
+
+                if(associated(DTDT)) DTDT(I,J,L) = DTDT(I,J,L) + DP_s*VR_s
+                ! DTDT(I,J,L) = DTDT(I,J,L) + DP_s*VR_s
+                if(FriendlyTemp    ) T   (I,J,L) = T   (I,J,L) + DT*VR_s
+             end if
+          enddo
+       enddo
+    enddo
+    !$acc end parallel loop
+    call cpu_time(t2)
+    !!$acc update host (DTDT, DUDT, DVDT, T_EQ, T, U, V)
+
+    !$acc end data
+
+    ! write(*,*) 'sum(T_EQ) = ', sum(T_EQ)
+
+    write(*,*) 'Runtime = ', t2 - t1
+    ! write(*,*) '***'
+
+  end subroutine held_suarez_oacc
+
+end module hs_oacc_mod
diff --git a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
index 0fcc762fa..013e95b2d 100644
--- a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
+++ b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
@@ -9,6 +9,12 @@ set (alldirs
 
 if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_SuperdynGridComp.F90)
 
+  if (NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOSdatmodyn_GridComp)
+
+    set (rename_GEOSdatmodyn_GridComp GEOS_DatmodynGridComp)
+
+  endif ()
+
   esma_add_library (${this}
     SRCS GEOS_SuperdynGridComp.F90
     SUBCOMPONENTS ${alldirs}

From 424deaae67fc12da2b4baf2caff8531c54e0b01f Mon Sep 17 00:00:00 2001
From: Purnendu Chakraborty <purnendu.chakraborty@nasa.gov>
Date: Wed, 17 May 2023 14:25:55 -0400
Subject: [PATCH 02/13] Created alarm ReplayShutOff which is required by
 FVdycoreCubed GridComp

---
 GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90 | 80 ++++++++++++-------
 1 file changed, 50 insertions(+), 30 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90 b/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
index 66c5099cd..f79e52f57 100644
--- a/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
+++ b/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90
@@ -13,8 +13,8 @@ module GEOS_AgcmSimpleGridCompMod
 !
 ! !DESCRIPTION: This gridded component (GC) combines the the GC that
 !   implements the Finite-Volume (FV) dynamics, with a simple physics
-!   component that implements the Held-Suarez benchmark forcing for 
-!   testing dry dynamical cores. 
+!   component that implements the Held-Suarez benchmark forcing for
+!   testing dry dynamical cores.
 !
 ! !USES:
 
@@ -35,7 +35,7 @@ module GEOS_AgcmSimpleGridCompMod
 
   public SetServices
 
- 
+
 !EOP
 
 contains
@@ -49,7 +49,7 @@ module GEOS_AgcmSimpleGridCompMod
 !   SetServices merely creates the children through MAPL and connects
 !   their Import-Export states. FV has some Imports that are
 !   not used in the Held-Suarez benchmark. Since these have proper defaults
-!   in FV, they can simply be terminated here. 
+!   in FV, they can simply be terminated here.
 
 
 ! !INTERFACE:
@@ -111,7 +111,7 @@ subroutine SetServices ( GC, RC )
          default    = 1.0e-6,                                      &
 !         RESTART    = MAPL_RestartRequired,                        &
          DIMS       = MAPL_DimsHorzVert,                           &
-         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
     VERIFY_(STATUS)
 
     call MAPL_AddInternalSpec(GC,                                        &
@@ -120,7 +120,7 @@ subroutine SetServices ( GC, RC )
          UNITS      = 'kg kg-1',                                         &
          FRIENDLYTO = 'DYNAMICS',                             &
          DIMS       = MAPL_DimsHorzVert,                                 &
-         VLOCATION  = MAPL_VLocationCenter,                   RC=STATUS  )  
+         VLOCATION  = MAPL_VLocationCenter,                   RC=STATUS  )
     VERIFY_(STATUS)
 
     call MAPL_AddInternalSpec(GC,                                       &
@@ -129,15 +129,15 @@ subroutine SetServices ( GC, RC )
          UNITS      = 'kg kg-1',                                        &
          FRIENDLYTO = 'DYNAMICS',                            &
          DIMS       = MAPL_DimsHorzVert,                                &
-         VLOCATION  = MAPL_VLocationCenter,                  RC=STATUS  )  
-    VERIFY_(STATUS)                                                                          
+         VLOCATION  = MAPL_VLocationCenter,                  RC=STATUS  )
+    VERIFY_(STATUS)
     call MAPL_AddInternalSpec(GC,                                     &
          SHORT_NAME = 'QILS',                                         &
          LONG_NAME  = 'mass_fraction_of_large_scale_cloud_ice_water', &
          UNITS      = 'kg kg-1',                                      &
          FRIENDLYTO = 'DYNAMICS',                          &
          DIMS       = MAPL_DimsHorzVert,                              &
-         VLOCATION  = MAPL_VLocationCenter,                RC=STATUS  )  
+         VLOCATION  = MAPL_VLocationCenter,                RC=STATUS  )
     VERIFY_(STATUS)
 
     call MAPL_AddInternalSpec(GC,                                    &
@@ -146,7 +146,7 @@ subroutine SetServices ( GC, RC )
          UNITS      = 'kg kg-1',                                     &
          FRIENDLYTO = 'DYNAMICS',                         &
          DIMS       = MAPL_DimsHorzVert,                             &
-         VLOCATION  = MAPL_VLocationCenter,               RC=STATUS  )  
+         VLOCATION  = MAPL_VLocationCenter,               RC=STATUS  )
     VERIFY_(STATUS)
 
     call MAPL_AddInternalSpec(GC,                                  &
@@ -155,8 +155,8 @@ subroutine SetServices ( GC, RC )
          UNITS      = '1',                                         &
          FRIENDLYTO = 'DYNAMICS',                                  &
          DIMS       = MAPL_DimsHorzVert,                           &
-         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
-    VERIFY_(STATUS)                                                                          
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+    VERIFY_(STATUS)
 
     call MAPL_AddInternalSpec(GC,                                  &
          SHORT_NAME = 'CLCN',                                      &
@@ -164,11 +164,11 @@ subroutine SetServices ( GC, RC )
          UNITS      = '1',                                         &
          FRIENDLYTO = 'DYNAMICS',                                  &
          DIMS       = MAPL_DimsHorzVert,                           &
-         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )  
-    VERIFY_(STATUS)                                                                          
+         VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
+    VERIFY_(STATUS)
     call MAPL_AddInternalSpec(GC,                                  &
          SHORT_NAME = 'QRAIN',                                     &
-         LONG_NAME  = 'mass_fraction_of_rain',                     & 
+         LONG_NAME  = 'mass_fraction_of_rain',                     &
          UNITS      = 'kg kg-1',                                   &
          FRIENDLYTO = 'DYNAMICS',                       &
          default    = 0.0,                                         &
@@ -195,7 +195,7 @@ subroutine SetServices ( GC, RC )
          DIMS       = MAPL_DimsHorzVert,                           &
          VLOCATION  = MAPL_VLocationCenter,             RC=STATUS  )
     VERIFY_(STATUS)
-                                                                          
+
 
 ! Register children with MAPL and go down their SS hierarchy
 ! ----------------------------------------------------------
@@ -245,8 +245,9 @@ subroutine SetServices ( GC, RC )
 
     call MAPL_GenericSetServices( GC, RC=STATUS )
     VERIFY_(STATUS)
-   
+
     RETURN_(ESMF_SUCCESS)
+
   end subroutine SetServices
 !EOC
 
@@ -254,7 +255,7 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC )
 
 ! !ARGUMENTS:
 
-    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component
     type(ESMF_State),    intent(inout) :: IMPORT ! Import state
     type(ESMF_State),    intent(inout) :: EXPORT ! Export state
     type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
@@ -264,16 +265,19 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC )
 
 ! ErrLog Variables
 
-    character(len=ESMF_MAXSTR)              :: IAm 
+    character(len=ESMF_MAXSTR)              :: IAm
     integer                                 :: STATUS
     character(len=ESMF_MAXSTR)              :: COMP_NAME
 
 ! local vars
-   type (MAPL_MetaComp),      pointer  :: MAPL
-   type (ESMF_State),         pointer  :: GIM(:)
-   type (ESMF_Field)                   :: FIELD
-   type (ESMF_FieldBundle)             :: BUNDLE
-   type (ESMF_Config)                  :: cf
+    type (MAPL_MetaComp),      pointer  :: MAPL
+    type (ESMF_State),         pointer  :: GIM(:)
+    type (ESMF_Field)                   :: FIELD
+    type (ESMF_FieldBundle)             :: BUNDLE
+    type (ESMF_Config)                  :: cf
+    type (ESMF_Alarm)                   :: replay_shutoff_alarm
+    type(ESMF_TimeInterval)             :: shutoff
+    integer                             :: rplshut
 
     Iam = "Initialize"
     call ESMF_GridCompGet ( GC, name=COMP_NAME, config=cf, RC=STATUS )
@@ -311,38 +315,52 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC )
     call MAPL_GridCompGetFriendlies(GC, "DYNAMICS", BUNDLE, RC=STATUS )
     VERIFY_(STATUS)
 
+    ! Initialize alarms
+    call MAPL_GetResource(MAPL, rplshut, Label="REPLAY_SHUTOFF:", default=-3600, rc=status)
+    VERIFY_(status)
+    call ESMF_TimeIntervalSet(shutoff, S=abs(rplshut), rc=status)
+    VERIFY_(status)
+    replay_shutoff_alarm = ESMF_AlarmCreate( &
+         name="ReplayShutOff", &
+         clock=clock, &
+         ringInterval=shutoff, &
+         sticky=.true., &
+         rc=status)
+    VERIFY_(status)
+
     RETURN_(ESMF_SUCCESS)
+
   end subroutine Initialize
 
   subroutine Run( GC, IMPORT, EXPORT, CLOCK, RC )
 
 ! !ARGUMENTS:
 
-    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component 
+    type(ESMF_GridComp), intent(inout) :: GC     ! Gridded component
     type(ESMF_State),    intent(inout) :: IMPORT ! Import state
     type(ESMF_State),    intent(inout) :: EXPORT ! Export state
     type(ESMF_Clock),    intent(inout) :: CLOCK  ! The clock
     integer, optional,   intent(  out) :: RC     ! Error code
 
- 
+
 !EOP
 
 ! ErrLog Variables
 
-    character(len=ESMF_MAXSTR)        :: IAm 
+    character(len=ESMF_MAXSTR)        :: IAm
     integer                           :: STATUS
     character(len=ESMF_MAXSTR)        :: COMP_NAME
-    
+
 ! Local derived type aliases
 
     type (MAPL_MetaComp    ), pointer :: MAPL
     type (ESMF_GridComp),      pointer  :: GCS(:)
     type (ESMF_State),         pointer  :: GIM(:)
     type (ESMF_State),         pointer  :: GEX(:)
-    
+
 !=============================================================================
 
-! Begin... 
+! Begin...
 
 ! Get the target components name and set-up traceback handle.
 ! -----------------------------------------------------------
@@ -384,6 +402,8 @@ subroutine Run( GC, IMPORT, EXPORT, CLOCK, RC )
     call MAPL_TimerOff(MAPL,"TOTAL")
 
     RETURN_(ESMF_SUCCESS)
+
   end subroutine Run
+
 end module GEOS_AgcmSimpleGridCompMod
 

From b797bacb73f28a40352943656ce6778a2a85f57c Mon Sep 17 00:00:00 2001
From: Purnendu Chakraborty <purnendu.chakraborty@nasa.gov>
Date: Wed, 17 May 2023 14:27:04 -0400
Subject: [PATCH 03/13] Removed executable mode from GEOShs.rc

---
 GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 mode change 100755 => 100644 GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc

diff --git a/GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc b/GEOSagcm_GridComp/GEOShs_GridComp/GEOShs.rc
old mode 100755
new mode 100644

From 2565e5a108f8b16fa3d1c9273f2fcfc882439b85 Mon Sep 17 00:00:00 2001
From: Weiyuan Jiang <weiyuan.jiang@nasa.gov>
Date: Thu, 18 May 2023 13:49:28 -0400
Subject: [PATCH 04/13] allocate ity

---
 .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90      | 1 +
 1 file changed, 1 insertion(+)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index 477e17d08..ae2f7442c 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -424,6 +424,7 @@ subroutine allocate_catch(this,rc)
      allocate( this%      cdcr1(ntiles) )
      allocate( this%      cdcr2(ntiles) )
      allocate( this%      poros(ntiles) )
+     allocate( this%      ity(ntiles) )
 
      allocate( this%         tc(ntiles,4) )
      allocate( this%         qc(ntiles,4) )

From ff94c474df784dadbb72f3dc5002782037413545 Mon Sep 17 00:00:00 2001
From: Purnendu Chakraborty <purnendu.chakraborty@nasa.gov>
Date: Thu, 18 May 2023 17:15:48 -0400
Subject: [PATCH 05/13] Rolled back change to GEOSgcm_GridComp/CMakeLists.txt

A better solution for GEOShs is to not include GEOS_GcmGridComp.F90 at all. In that case the current CMakeLists.txt works just fine
---
 CMakeLists.txt | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index fecd215ba..4c9699c3e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,11 +7,8 @@ set (alldirs
   GEOSogcm_GridComp
   )
 
-if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOSagcm_GridComp/GEOS_AgcmSimpleGridComp.F90)
 
-  esma_add_subdirectories (GEOSagcm_GridComp)
-
-elseif (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_GcmGridComp.F90)
+if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_GcmGridComp.F90)
 
  ecbuild_declare_project()
   esma_add_library(${this}

From e01be3a4cb1dd0835951527b2030fbe7396845ff Mon Sep 17 00:00:00 2001
From: Purnendu Chakraborty <purnendu.chakraborty@nasa.gov>
Date: Thu, 18 May 2023 18:49:39 -0400
Subject: [PATCH 06/13] Better handling for when GEOSdatmodyn GridComp does not
 exist

---
 .../GEOSsuperdyn_GridComp/CMakeLists.txt             | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
index 013e95b2d..1a45ae14a 100644
--- a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
+++ b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/CMakeLists.txt
@@ -7,13 +7,13 @@ set (alldirs
   ARIESg3_GridComp
   )
 
-if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_SuperdynGridComp.F90)
+if (NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOSdatmodyn_GridComp)
 
-  if (NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOSdatmodyn_GridComp)
+  set (rename_GEOSdatmodyn_GridComp GEOS_DatmodynGridComp)
 
-    set (rename_GEOSdatmodyn_GridComp GEOS_DatmodynGridComp)
+endif ()
 
-  endif ()
+if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/GEOS_SuperdynGridComp.F90)
 
   esma_add_library (${this}
     SRCS GEOS_SuperdynGridComp.F90
@@ -25,7 +25,3 @@ else ()
   esma_add_subdirectories (${alldirs})
 
 endif ()
-
-
-
-

From fc8f57c58794db801faec7e27ec39fc0b06ddd14 Mon Sep 17 00:00:00 2001
From: Weiyuan Jiang <weiyuan.jiang@nasa.gov>
Date: Mon, 22 May 2023 09:46:46 -0400
Subject: [PATCH 07/13] ignore restarts' ity

---
 .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90   | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index ae2f7442c..9fd213fc4 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -148,7 +148,6 @@ function CatchmentRst_create(filename, time, rc) result (catch)
        if (myid ==0) then
           call catch%allocate_catch()
           call catch%read_shared_nc4(formatter, __RC__)
-          call MAPL_VarRead(formatter,"OLD_ITY",catch%ity, __RC__)
        endif
        call formatter%close()
      else
@@ -229,7 +228,7 @@ subroutine read_GEOSldas_rst_bin(this, filename, rc)
       read(unit) this%     atau
       read(unit) this%     btau
       read(unit) this%      ity
-      read(unit) this%       tc
+      read(unit) !this%       tc
       read(unit) this%       qc
       read(unit) this%    capac
       read(unit) this%   catdef
@@ -424,7 +423,6 @@ subroutine allocate_catch(this,rc)
      allocate( this%      cdcr1(ntiles) )
      allocate( this%      cdcr2(ntiles) )
      allocate( this%      poros(ntiles) )
-     allocate( this%      ity(ntiles) )
 
      allocate( this%         tc(ntiles,4) )
      allocate( this%         qc(ntiles,4) )

From 3da9f93f0f1ce253ca75d1cee92254e85fe3a8ed Mon Sep 17 00:00:00 2001
From: Weiyuan Jiang <52509753+weiyuan-jiang@users.noreply.github.com>
Date: Mon, 22 May 2023 10:04:59 -0400
Subject: [PATCH 08/13] Update CatchmentRst.F90

comment out the wrong variable
---
 .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90   | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index 9fd213fc4..9fc6f3745 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -227,8 +227,8 @@ subroutine read_GEOSldas_rst_bin(this, filename, rc)
       read(unit) this%     tsb2
       read(unit) this%     atau
       read(unit) this%     btau
-      read(unit) this%      ity
-      read(unit) !this%       tc
+      read(unit) !this%      ity
+      read(unit) this%       tc
       read(unit) this%       qc
       read(unit) this%    capac
       read(unit) this%   catdef

From 247623a3ada00c27209a4623ed62bd522d2a009f Mon Sep 17 00:00:00 2001
From: Matthew Thompson <matthew.thompson@nasa.gov>
Date: Thu, 25 May 2023 13:30:23 -0400
Subject: [PATCH 09/13] Use -O2 optimization in moist with GNU Release

---
 .../GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt    | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt
index a1f0435e1..af3e15f1f 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/CMakeLists.txt
@@ -20,6 +20,10 @@ if (CMAKE_Fortran_COMPILER_ID MATCHES Intel AND CMAKE_BUILD_TYPE MATCHES Aggress
    set (CMAKE_Fortran_FLAGS_AGGRESSIVE "${GEOS_Fortran_FLAGS_VECT}")
 endif ()
 
+if (CMAKE_Fortran_COMPILER_ID MATCHES GNU AND CMAKE_BUILD_TYPE MATCHES Release)
+   string (REPLACE "${FOPT3}" "${FOPT2}" CMAKE_Fortran_FLAGS_RELEASE ${CMAKE_Fortran_FLAGS_RELEASE})
+endif ()
+
 # Note For unknown reasons, BACM_1M_Interface takes 20 minutes to compile at O3
 #      and 10 minutes at O2. But only 7 seconds with O1. So we compile at O1
 if (CMAKE_Fortran_COMPILER_ID MATCHES Intel AND CMAKE_BUILD_TYPE MATCHES Release)

From cc7800e75ec119a87f9bc922fd3ebd46708c6110 Mon Sep 17 00:00:00 2001
From: Matthew Thompson <matthew.thompson@nasa.gov>
Date: Wed, 31 May 2023 12:22:40 -0400
Subject: [PATCH 10/13] Update CI Baselibs

---
 .circleci/config.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index cc09c5c54..260f043f6 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -1,7 +1,7 @@
 version: 2.1
 
 # Anchors to prevent forgetting to update a version
-baselibs_version: &baselibs_version v7.7.0
+baselibs_version: &baselibs_version v7.13.0
 bcs_version: &bcs_version v11.00.0
 
 orbs:

From 533953283ca49ea0db0795fd4abcefc949af057d Mon Sep 17 00:00:00 2001
From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com>
Date: Fri, 2 Jun 2023 08:59:21 -0400
Subject: [PATCH 11/13] added comment to clarify recent commit in
 CatchmentRst.F90

---
 .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90     | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index 9fc6f3745..b4da5ce66 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -227,7 +227,7 @@ subroutine read_GEOSldas_rst_bin(this, filename, rc)
       read(unit) this%     tsb2
       read(unit) this%     atau
       read(unit) this%     btau
-      read(unit) !this%      ity
+      read(unit)                  ! skip over "ity"
       read(unit) this%       tc
       read(unit) this%       qc
       read(unit) this%    capac

From 3eee6b00b9b64a6383ffeb0a5f351219d152dec6 Mon Sep 17 00:00:00 2001
From: Weiyuan Jiang <weiyuan.jiang@nasa.gov>
Date: Fri, 2 Jun 2023 13:42:44 -0400
Subject: [PATCH 12/13] skip more varaiables that are directly from BC

---
 .../Utils/mk_restarts/CatchmentRst.F90        | 58 +++++++++----------
 1 file changed, 29 insertions(+), 29 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index b4da5ce66..b50f8603b 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -198,36 +198,36 @@ subroutine read_GEOSldas_rst_bin(this, filename, rc)
       integer, optional, intent(out) :: rc
       integer :: unit
       open(newunit=unit, file=filename, form='unformatted', action='read')
-      read(unit) this%      bf1
-      read(unit) this%      bf2
-      read(unit) this%      bf3
+      read(unit) ! skip  bf1
+      read(unit) ! skip  bf2
+      read(unit) ! skip  bf3
       read(unit) this%   vgwmax
-      read(unit) this%    cdcr1
-      read(unit) this%    cdcr2
-      read(unit) this%     psis
-      read(unit) this%      bee
-      read(unit) this%    poros
-      read(unit) this%    wpwet
-      read(unit) this%     cond
-      read(unit) this%      gnu
-      read(unit) this%     ars1
-      read(unit) this%     ars2
-      read(unit) this%     ars3
-      read(unit) this%     ara1
-      read(unit) this%     ara2
-      read(unit) this%     ara3
-      read(unit) this%     ara4
-      read(unit) this%     arw1
-      read(unit) this%     arw2
-      read(unit) this%     arw3
-      read(unit) this%     arw4
-      read(unit) this%     tsa1
-      read(unit) this%     tsa2
-      read(unit) this%     tsb1
-      read(unit) this%     tsb2
-      read(unit) this%     atau
-      read(unit) this%     btau
-      read(unit)                  ! skip over "ity"
+      read(unit) this%   cdcr1
+      read(unit) this%   cdcr2
+      read(unit) ! skip  psis
+      read(unit) ! skip  bee
+      read(unit) this%   poros
+      read(unit) ! skip  wpwet
+      read(unit) ! skip    cond
+      read(unit) ! skip     gnu
+      read(unit) ! skip    ars1
+      read(unit) ! skip    ars2
+      read(unit) ! skip    ars3
+      read(unit) ! skip    ara1
+      read(unit) ! skip    ara2
+      read(unit) ! skip    ara3
+      read(unit) ! skip    ara4
+      read(unit) ! skip    arw1
+      read(unit) ! skip    arw2
+      read(unit) ! skip    arw3
+      read(unit) ! skip    arw4
+      read(unit) ! skip    tsa1
+      read(unit) ! skip    tsa2
+      read(unit) ! skip    tsb1
+      read(unit) ! skip    tsb2
+      read(unit) ! skip    atau
+      read(unit) ! skip    btau
+      read(unit) ! skip    ity
       read(unit) this%       tc
       read(unit) this%       qc
       read(unit) this%    capac

From fbc9b80aa6d3eba95e17d0dc361c298e6eb61cc6 Mon Sep 17 00:00:00 2001
From: Rolf Reichle <rolf.reichle@nasa.gov>
Date: Mon, 5 Jun 2023 09:29:51 -0400
Subject: [PATCH 13/13] white-space changes to improve readability
 (mk_restarts/CatchmentRst.F90)

---
 .../Utils/mk_restarts/CatchmentRst.F90        | 58 +++++++++----------
 1 file changed, 29 insertions(+), 29 deletions(-)

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
index b50f8603b..e8ef24138 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90
@@ -198,36 +198,36 @@ subroutine read_GEOSldas_rst_bin(this, filename, rc)
       integer, optional, intent(out) :: rc
       integer :: unit
       open(newunit=unit, file=filename, form='unformatted', action='read')
-      read(unit) ! skip  bf1
-      read(unit) ! skip  bf2
-      read(unit) ! skip  bf3
+      read(unit)                  ! skip     bf1
+      read(unit)                  ! skip     bf2
+      read(unit)                  ! skip     bf3
       read(unit) this%   vgwmax
-      read(unit) this%   cdcr1
-      read(unit) this%   cdcr2
-      read(unit) ! skip  psis
-      read(unit) ! skip  bee
-      read(unit) this%   poros
-      read(unit) ! skip  wpwet
-      read(unit) ! skip    cond
-      read(unit) ! skip     gnu
-      read(unit) ! skip    ars1
-      read(unit) ! skip    ars2
-      read(unit) ! skip    ars3
-      read(unit) ! skip    ara1
-      read(unit) ! skip    ara2
-      read(unit) ! skip    ara3
-      read(unit) ! skip    ara4
-      read(unit) ! skip    arw1
-      read(unit) ! skip    arw2
-      read(unit) ! skip    arw3
-      read(unit) ! skip    arw4
-      read(unit) ! skip    tsa1
-      read(unit) ! skip    tsa2
-      read(unit) ! skip    tsb1
-      read(unit) ! skip    tsb2
-      read(unit) ! skip    atau
-      read(unit) ! skip    btau
-      read(unit) ! skip    ity
+      read(unit) this%    cdcr1
+      read(unit) this%    cdcr2
+      read(unit)                  ! skip    psis
+      read(unit)                  ! skip     bee
+      read(unit) this%    poros
+      read(unit)                  ! skip   wpwet
+      read(unit)                  ! skip    cond
+      read(unit)                  ! skip     gnu
+      read(unit)                  ! skip    ars1
+      read(unit)                  ! skip    ars2
+      read(unit)                  ! skip    ars3
+      read(unit)                  ! skip    ara1
+      read(unit)                  ! skip    ara2
+      read(unit)                  ! skip    ara3
+      read(unit)                  ! skip    ara4
+      read(unit)                  ! skip    arw1
+      read(unit)                  ! skip    arw2
+      read(unit)                  ! skip    arw3
+      read(unit)                  ! skip    arw4
+      read(unit)                  ! skip    tsa1
+      read(unit)                  ! skip    tsa2
+      read(unit)                  ! skip    tsb1
+      read(unit)                  ! skip    tsb2
+      read(unit)                  ! skip    atau
+      read(unit)                  ! skip    btau
+      read(unit)                  ! skip     ity
       read(unit) this%       tc
       read(unit) this%       qc
       read(unit) this%    capac