OBCS_SEAICE_NEUMANN flag, more store directives, and cleanup

doc/phys_pkgs/obcs.rst

@@ -548,10 +548,10 @@ OBCS\_APPLY\_*:
 :filelink:`OBCS\_SPONGE <pkg/obcs/obcs_sponge.F>`:
 ##################################################
 
-The sponge layer code (turned on with CPP option :varlink:`ALLOW_OBCS_SPONGE` and run-time parameter
-:varlink:`useOBCSsponge`) adds a relaxation term to the right-hand-side of
-the momentum and tracer equations. The variables are relaxed towards
-the boundary values with a relaxation time scale that increases
+The sponge layer code (turned on with CPP option :varlink:`ALLOW_OBCS_SPONGE`
+and run-time parameter :varlink:`useOBCSsponge`) adds a relaxation term to the
+right-hand-side of the momentum and tracer equations. The variables are relaxed
+towards the boundary values with a relaxation time scale that increases
 linearly with distance from the boundary
 
 .. math::
@@ -563,27 +563,42 @@ linearly with distance from the boundary
    {[(1-l)\tau_{b}+l\tau_{i}]}
 
 where :math:`\chi` is the model variable (U/V/T/S) in the interior,
-:math:`\chi_{BC}` the boundary value, :math:`L` the thickness of the
-sponge layer (runtime parameter :varlink:`spongeThickness` in number
-of grid points), :math:`\delta{L}\in[0,L]`
-(:math:`\frac{\delta{L}}{L}=l\in[0,1]`) the distance from the boundary
-(also in grid points), and :math:`\tau_{b}` (runtime parameters
-:varlink:`Urelaxobcsbound` and :varlink:`Vrelaxobcsbound`) and
+:math:`\chi_{BC}` the boundary value, :math:`L` the thickness of the sponge
+layer (runtime parameter :varlink:`spongeThickness` in number of grid points),
+:math:`\delta{L}\in[0,L]` (:math:`\frac{\delta{L}}{L}=l\in[0,1]`) the distance
+from the boundary (also in grid points), and :math:`\tau_{b}` (runtime
+parameters :varlink:`Urelaxobcsbound` and :varlink:`Vrelaxobcsbound`) and
 :math:`\tau_{i}` (runtime parameters :varlink:`Urelaxobcsinner` and
-:varlink:`Vrelaxobcsinner`) the relaxation time scales on the boundary
-and at the interior termination of the sponge layer. The parameters
-:varlink:`Urelaxobcsbound` and :varlink:`Urelaxobcsinner` set the relaxation time scales for
-the Eastern and Western boundaries, :varlink:`Vrelaxobcsbound` and :varlink:`Vrelaxobcsinner`
-for the Northern and Southern boundaries.
+:varlink:`Vrelaxobcsinner`) the relaxation time scales on the boundary and at
+the interior termination of the sponge layer. The parameters
+:varlink:`Urelaxobcsbound` and :varlink:`Urelaxobcsinner` set the relaxation
+time scales for the Eastern and Western boundaries, :varlink:`Vrelaxobcsbound`
+and :varlink:`Vrelaxobcsinner` for the Northern and Southern boundaries.
 
 
 OB's with nonlinear free surface
 ################################
 
-
 OB's with sea ice
 #################
 
+Simple Dirichlet boundary conditions for sea ice parameters can be specified in
+anology to the ocean variables via filenames ``OB[N/S/E/W][a/h/sl/sn/u/v]File``
+(sea ice concentration, cell averaged sea ice thickness, salinity, cell
+averaged snow thickness, ice drift components). With CPP-flag
+:varlink:`ALLOW_OBCS_SEAICE_SPONGE` and runtime flags
+:varlink:`useSeaiceSponge`, :varlink:`seaiceSpongeThickness`, and
+``[A/H/SL/SN]relaxobcs[inner/bound]`` are available in analogy to the sponge
+parameters for the ocean variables.
+
+Neumann boundary conditions :math:`\frac{\partial\phi}{\partial{n}}=0` for all
+sea ice variables can be applied with runtime flag
+:varlink:`SEAICEuseNeumannBC`, which overrides the input files for the
+Dirichlet values.
+
+Defining CPP-flag :varlink:`OBCS_SEAICE_SMOOTH_EDGE` allows to smooth the
+tracer sea-ice variables near the edges.
+
 
 .. _ssub_phys_pkg_obcs_flowchart:
 
@@ -595,21 +610,59 @@ Flow chart
 
 
     C     !CALLING SEQUENCE:
-    c ...
+    C    [...]
+    C    | |-MAIN_DO_LOOP    :: Open-AD case: Main timestepping loop routine
+    C    | \                    otherwise: just call FORWARD_STEP
+    C    | |
+    C/\  | |-FORWARD_STEP        :: Step forward a time-step ( AT LAST !!! )
+    C    [...]
+    C/\  | | |-DO_OCEANIC_PHYS   :: Oceanic (& seaice) physics computation
+    C/\  | | | |
+    C/\  | | | |-OBCS_CALC       :: Open boundary. package (see pkg/obcs).
+    C/\  | | | |
+    C    [...]
+    C/\  | | | |-SEAICE_MODEL          :: pkg/seaice
+    C/\  | | | | |-SEAICE_DYNSOLVER    :: pkg/seaice
+    C/\  | | | | | |-OBCS_APPLY_UVICE  :: apply uIce/vIce boudnary conditions
+    C/\  | | | | |-OBCS_ADJUST_UVICE   :: (Only for OBCS_UVICE_OLD)
+    C/\  | | | | |-SEAICE_GROWTH
+    C/\  | | | | |-SEAICE_APPLY_SEAICE :: add OBCS for scalar variables
+    C    [...]
+    C/\  | | |-THERMODYNAMICS         :: theta, salt + tracer equations driver.
+    C/\  | | | |                         (synchronous time-stepping case)
+    C    [...]
+    C/\  | | | |-TEMP_INTEGRATE       :: Step forward Prognostic Eq for Temperature.
+    C/\  | | | |
+    C/\  | | | |-SALT_INTEGRATE       :: Step forward Prognostic Eq for Salinity.
+    C/\  | | | |                         same sequence of calls as in TEMP_INTEGRATE
+    C/\  | | | |
+    C/\  | | | |-PTRACERS_INTEGRATE   :: Integrate other tracer(s) (see pkg/ptracers).
+    C/\  | | | | |                     same sequence of calls as in TEMP_INTEGRATE
+    C/\  | | | | |-OBCS_APPLY_PTRACER :: Open boundary package for pTracers
+    C/\  | | | |
+    C/\  | | | |-OBCS_APPLY_TS        :: Open boundary package (see pkg/obcs ).
+    C/\  | | |
+    C    [...]
+    C/\  | | |
+    C/\  | | |-DYNAMICS       :: Momentum equations driver.
+    C/\  | | | |
+    C    [...]
+    C/\  | | | |-OBCS_APPLY_UV    :: Apply Open bndary Conditions to provisional U,V
+    C    [...]
+    C/\  | | |-MOMENTUM_CORRECTION_STEP :: Finalise momentum stepping
+    C    [...]
+    C/\  | | | |-OBCS_APPLY_UV       :: Open boundary package (see pkg/obcs).
+
 
 
 .. _ssub_phys_pkg_obcs_diagnostics:
 
 OBCS diagnostics
 ++++++++++++++++
 
-Diagnostics output is available via the diagnostics package (see :numref:`sub_outp_pkg_diagnostics`). Available output fields are summarized below:
-
-::
-
-    ------------------------------------------------------
-     <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
-    ------------------------------------------------------
+Diagnostics output is available via the diagnostics package (see
+:numref:`sub_outp_pkg_diagnostics`). Currently there are no OBCS-specific
+diagnostics available.
 
 
 .. _ssub_phys_pkg_obcs_experiments:

model/src/do_oceanic_phys.F

@@ -32,6 +32,9 @@
 # ifdef ALLOW_EXF
 #  include "EXF_OPTIONS.h"
 # endif
+#ifdef ALLOW_OBCS
+# include "OBCS_OPTIONS.h"
+#endif
 #endif /* ALLOW_AUTODIFF */
 
 CBOP
@@ -188,6 +191,10 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 # ifdef ALLOW_SALT_PLUME
 #  include "SALT_PLUME.h"
 # endif
+# ifdef ALLOW_OBCS
+#  include "OBCS_PARAMS.h"
+#  include "OBCS_FIELDS.h"
+# endif
 #endif /* ALLOW_AUTODIFF */
 
 #ifdef ALLOW_TAPENADE
@@ -238,16 +245,10 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 #endif /* ALLOW_AUTODIFF */
 #ifdef ALLOW_AUTODIFF_TAMC
 C     tkey :: tape key (tile dependent)
-C     kkey :: tape key (level and tile dependent)
-      INTEGER kkey, tkey
+      INTEGER tkey
 #endif
 CEOP
 
-#ifdef ALLOW_AUTODIFF_TAMC
-C--   dummy statement to end declaration part
-      tkey = 1
-#endif /* ALLOW_AUTODIFF_TAMC */
-
       kSrf = 1
       IF ( usingPCoords ) kSrf = Nr
 
@@ -287,13 +288,50 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 # ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
 CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE etaN  = comlev1, key=ikey_dynamics, kind=isbyte
+#  ifdef ALLOW_OBCS_STEVENS
+CADJ STORE uVel  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE vVel  = comlev1, key=ikey_dynamics, kind=isbyte
+#   ifdef ALLOW_OBCS_EAST
+CADJ STORE OBEtStevens = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE OBEsStevens = comlev1, key=ikey_dynamics, kind=isbyte
+#   endif
+#   ifdef ALLOW_OBCS_WEST
+CADJ STORE OBWtStevens = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE OBWsStevens = comlev1, key=ikey_dynamics, kind=isbyte
+#   endif
+#   ifdef ALLOW_OBCS_NORTH
+CADJ STORE OBNtStevens = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE OBNsStevens = comlev1, key=ikey_dynamics, kind=isbyte
+#   endif
+#   ifdef ALLOW_OBCS_SOUTH
+CADJ STORE OBStStevens = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE OBSsStevens = comlev1, key=ikey_dynamics, kind=isbyte
+#   endif
+#  endif /* ALLOW_OBCS_STEVENS */
 # endif
 # ifdef ALLOW_DEBUG
        IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
 # endif
        CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
      I                 uVel, vVel, wVel, theta, salt, myThid )
       ENDIF
+# if ( defined ALLOW_AUTODIFF_TAMC && defined ALLOW_OBCS_BALANCE )
+C     This needs to be done ***after*** the if-block to avoid calling
+C     S/R OBCS_CALC in the AD code.
+#  ifdef ALLOW_OBCS_NORTH
+CADJ STORE OBNv        = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+#  ifdef ALLOW_OBCS_SOUTH
+CADJ STORE OBSv        = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+#  ifdef ALLOW_OBCS_EAST
+CADJ STORE OBEu        = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+#  ifdef ALLOW_OBCS_WEST
+CADJ STORE OBWu        = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+# endif /* ALLOW_AUTODIFF_TAMC and ALLOW_OBCS_BALANCE */
 #endif  /* ALLOW_OBCS */
 
 #ifdef ALLOW_OCN_COMPON_INTERF
@@ -364,16 +402,16 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 # ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE qnet  = comlev1, key=ikey_dynamics, kind=isbyte
 CADJ STORE qsw   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
 CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
-#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
+#  if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
 CADJ STORE evap  = comlev1, key=ikey_dynamics, kind=isbyte
-#endif
-# endif /* ALLOW_AUTODIFF_TAMC */
-#ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE phiHydLow= comlev1, key=ikey_dynamics, byte=isbyte
-#endif
+#  endif
+CADJ STORE etan      = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE theta     = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE salt      = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel,vvel = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE phiHydLow = comlev1, key=ikey_dynamics, byte=isbyte
+# endif
       IF ( useSEAICE ) THEN
 # ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
@@ -421,7 +459,9 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 # endif
 # ifdef ALLOW_AUTODIFF
       ELSEIF ( SEAICEadjMODE .EQ. -1 ) THEN
+#  ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
         CALL SEAICE_FAKE( myTime, myIter, myThid )
 # endif /* ALLOW_AUTODIFF */
       ENDIF
@@ -505,7 +545,9 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 C---  Determines forcing terms based on external fields
 C     relaxation terms, etc.
 #ifdef ALLOW_AUTODIFF
+# ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
+# endif
 #else  /* ALLOW_AUTODIFF */
 C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
 C     and all vertical mixing schemes, but keep OBCS_CALC
@@ -517,6 +559,11 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
         CALL EXTERNAL_FORCING_SURF(
      I             iMin, iMax, jMin, jMax,
      I             myTime, myIter, myThid )
+#ifdef ALLOW_AUTODIFF_TAMC
+C     Avoid calling S/R EXTERNAL_FORCING_SURF in AD routine.
+CADJ STORE EmPmR             = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel, vvel        = comlev1, key=ikey_dynamics, kind=isbyte
+#endif
 
 #ifdef  ALLOW_OBCS
       IF (useOBCS) THEN
@@ -653,10 +700,6 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 CADJ STORE uvel (:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
 CADJ STORE vvel (:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
 # endif
-# ifdef ALLOW_SALT_PLUME
-CADJ STORE saltplumedepth(:,:,bi,bj)= comlev1_bibj,key=tkey,kind=isbyte
-CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
-# endif
 #endif /* ALLOW_AUTODIFF_TAMC */
 
 C--   Always compute density (stored in common block) here; even when it is not
@@ -679,9 +722,10 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 #endif /* ALLOW_DOWN_SLOPE */
 #ifdef ALLOW_BBL
          IF ( useBBL ) THEN
-C     pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
-C     To reduce computation and storage requirement, these densities are stored in the
-C     dry grid boxes of rhoInSitu.  See BBL_CALC_RHO for details.
+C     pkg/bbl requires in-situ bbl density for depths equal to and deeper
+C     than the bbl. To reduce computation and storage requirement,
+C     these densities are stored in the dry grid boxes of rhoInSitu.
+C     See BBL_CALC_RHO for details.
            DO k=Nr,1,-1
             CALL BBL_CALC_RHO(
      I                  theta, salt,
@@ -722,6 +766,9 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
            ENDDO
           ENDDO
         ENDIF
+# ifdef ALLOW_AUTODIFF_TAMC
+CADJ STORE rhoInSitu(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
+# endif
 #endif /* ALLOW_AUTODIFF */
 
 #ifdef ALLOW_DEBUG
@@ -735,33 +782,13 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
 C--     Start of diagnostic loop
         DO k=Nr,1,-1
 
-#ifdef ALLOW_AUTODIFF_TAMC
-C? Patrick, is this formula correct now that we change the loop range?
-C? Do we still need this?
-cph kkey formula corrected.
-cph Needed for rhoK, rhoKm1, in the case useGMREDI.
-          kkey = k + (tkey-1)*Nr
-#endif /* ALLOW_AUTODIFF_TAMC */
-
-c#ifdef ALLOW_AUTODIFF_TAMC
-cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
-cCADJ &     kind = isbyte
-cCADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey,
-cCADJ &     kind = isbyte
-c#endif /* ALLOW_AUTODIFF_TAMC */
-
 C--       Calculate gradients of potential density for isoneutral
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
           IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
      &         .OR. usePP81 .OR. useKL10
      &         .OR. useMY82 .OR. useGGL90
      &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
             IF (k.GT.1) THEN
-#ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,kind=isbyte
-CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,kind=isbyte
-CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,kind=isbyte
-#endif /* ALLOW_AUTODIFF_TAMC */
              IF ( usingZCoords ) THEN
               DO j=jMin,jMax
                DO i=iMin,iMax
@@ -840,6 +867,8 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
         ENDDO
 
 #ifdef ALLOW_AUTODIFF_TAMC
+C     Avoid recomputing sigmaR and IVDConvCount in AD routine.
+CADJ STORE sigmaR                   =comlev1_bibj,key=tkey,kind=isbyte
 CADJ STORE IVDConvCount(:,:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
 #endif
 
@@ -873,6 +902,10 @@ SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
      I              myTime, myIter, myThid )
 #endif /* SALT_PLUME_VOLUME */
         ENDIF
+# ifdef ALLOW_AUTODIFF_TAMC
+CADJ STORE saltplumedepth(:,:,bi,bj)= comlev1_bibj,key=tkey,kind=isbyte
+CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
+# endif /* ALLOW_AUTODIFF_TAMC */
 #endif /* ALLOW_SALT_PLUME */
 
 #ifdef ALLOW_DIAGNOSTICS

model/src/forward_step.F

@@ -640,10 +640,6 @@ SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
 # ifdef ALLOW_OBCS
 CADJ STORE salt       = comlev1, key = ikey_dynamics, kind = isbyte
 CADJ STORE totphihyd  = comlev1, key = ikey_dynamics, kind = isbyte
-#  ifdef ALLOW_OBCS_STEVENS
-CADJ STORE gsNm1      = comlev1, key = ikey_dynamics, kind = isbyte
-CADJ STORE gtNm1      = comlev1, key = ikey_dynamics, kind = isbyte
-#  endif
 # endif /* ALLOW_OBCS */
 # ifdef ALLOW_PTRACERS
 CADJ STORE pTracer    = comlev1, key = ikey_dynamics, kind = isbyte