Index: /issm/trunk-jpl/test/MITgcm/code/CPP_EEOPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code/CPP_EEOPTIONS.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code/CPP_EEOPTIONS.h	(revision 26544)
@@ -62,7 +62,4 @@
 C     Use D0 where it means REAL*8 but not where it means REAL*16
 #define D0 d0
-
-C--   Enable some old macro conventions for backward compatibility
-#undef USE_OLD_MACROS_R4R8toRSRL
 
 C=== IO related options ===
Index: /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F	(revision 26544)
@@ -1,4 +1,7 @@
 #include "PACKAGES_CONFIG.h"
 #include "CPP_OPTIONS.h"
+#ifdef ALLOW_MOM_COMMON
+# include "MOM_COMMON_OPTIONS.h"
+#endif
 #ifdef ALLOW_AUTODIFF
 # include "AUTODIFF_OPTIONS.h"
@@ -73,4 +76,6 @@
 C       |-- EXTERNAL_FORCING_SURF
 C       |
+C       |-- OBCS_ADJUST
+C       |
 C       |- k loop (Nr:1):
 C       | - DWNSLP_CALC_RHO
@@ -119,4 +124,6 @@
 C       |
 C       |-- KPP_DO_EXCH
+C       |
+C       |-- GGL90_EXCHANGES
 C       |
 C       |-- DIAGS_RHO_G
@@ -143,6 +150,7 @@
 #ifdef ALLOW_AUTODIFF
 # include "AUTODIFF_MYFIELDS.h"
-# include "tamc.h"
-# include "tamc_keys.h"
+# ifdef ALLOW_AUTODIFF_TAMC
+#  include "tamc.h"
+# endif
 # include "FFIELDS.h"
 # include "SURFACE.h"
@@ -196,5 +204,5 @@
 C     !LOCAL VARIABLES:
 C     == Local variables
-C     rhoK, rhoKm1  :: Density at current level, and level above
+C     rhoKp1,rhoKm1 :: Density at current level, and @ level minus one
 C     iMin, iMax    :: Ranges and sub-block indices on which calculations
 C     jMin, jMax       are applied.
@@ -202,4 +210,5 @@
 C     msgBuf        :: Temp. for building output string
 C     i,j,k         :: loop indices
+C     kSrf          :: surface index
       _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
@@ -210,5 +219,5 @@
       INTEGER jMin, jMax
       INTEGER bi, bj
-      INTEGER i, j, k
+      INTEGER i, j, k, kSrf
       CHARACTER*(MAX_LEN_MBUF) msgBuf
       INTEGER doDiagsRho
@@ -221,4 +230,9 @@
       _RL thetaRef
 #endif /* ALLOW_AUTODIFF */
+#ifdef ALLOW_AUTODIFF_TAMC
+      INTEGER act1, act2, act3, act4
+      INTEGER max1, max2, max3
+      INTEGER kkey, itdkey
+#endif
 CEOP
 
@@ -227,4 +241,7 @@
       itdkey = 1
 #endif /* ALLOW_AUTODIFF_TAMC */
+
+      kSrf = 1
+      IF ( usingPCoords ) kSrf = Nr
 
 #ifdef ALLOW_DEBUG
@@ -282,31 +299,17 @@
 
 #ifdef ALLOW_AUTODIFF
-# ifdef ALLOW_SALT_PLUME
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
+          adjustColdSST_diag(i,j,bi,bj) = 0. _d 0
+# ifdef ALLOW_SALT_PLUME
           saltPlumeDepth(i,j,bi,bj) = 0. _d 0
           saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
+# endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
-# endif
-# ifdef ALLOW_ECCO
-#  ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
-      DO bj=myByLo(myThid),myByHi(myThid)
-       DO bi=myBxLo(myThid),myBxHi(myThid)
-        DO k=1,Nr
-         DO j=1-OLy,sNy+OLy
-          DO i=1-OLx,sNx+OLx
-           sigmaRfield(i,j,k,bi,bj) = 0. _d 0
-          ENDDO
-         ENDDO
-        ENDDO
-       ENDDO
-      ENDDO
-#  endif
-# endif
 #endif /* ALLOW_AUTODIFF */
 
@@ -320,35 +323,22 @@
 #endif /* ALLOW_FRAZIL */
 
-#ifndef OLD_THSICE_CALL_SEQUENCE
 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
       IF ( useThSIce .AND. fluidIsWater ) THEN
 # ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE snowHeight, Tsrf  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE Qice1, Qice2      = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE salt,theta        = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uvel,vvel         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qnet,qsw, empmr   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE atemp,aqh,precip  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE swdown,lwdown     = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+#  ifdef ALLOW_SEAICE
+CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE snowHeight, Tsrf  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE sHeating,snowAge  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE hocemxl, icflxsw  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE salt,theta        = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE qnet,qsw, empmr   = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
 #  ifdef NONLIN_FRSURF
-CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE hFac_surfC        = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 # endif /* ALLOW_AUTODIFF_TAMC */
@@ -363,106 +353,45 @@
       ENDIF
 #endif /* ALLOW_THSICE */
-#endif /* ndef OLD_THSICE_CALL_SEQUENCE */
 
 #ifdef ALLOW_SEAICE
-# ifdef ALLOW_AUTODIFF
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
-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
+# 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)
-CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
-#endif
-      IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
-        CALL SEAICE_FAKE( myTime, myIter, myThid )
-      ENDIF
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
-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
-#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
-CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
-#endif
-# endif /* ALLOW_AUTODIFF */
-#endif /* ALLOW_SEAICE */
-
-#ifdef ALLOW_SEAICE
+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
       IF ( useSEAICE ) THEN
 # ifdef ALLOW_AUTODIFF_TAMC
-cph-adj-test(
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE tices  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE empmr, qnet  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
-cCADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
-cCADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
-cph-adj-test)
-c#ifdef ALLOW_EXF
-CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE evap                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uwind,vwind         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-c#endif
-CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  ifdef SEAICE_CGRID
-CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  endif
-#  ifdef SEAICE_ALLOW_DYNAMICS
-CADJ STORE uice                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE vice                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE dwatn               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#   ifdef SEAICE_ALLOW_EVP
-CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#   endif
+CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
+#  ifdef ALLOW_EXF
+CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uwind,vwind       = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef SEAICE_VARIABLE_SALINITY
-CADJ STORE hsalt               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE hsalt             = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef ATMOSPHERIC_LOADING
-CADJ STORE pload               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE siceload            = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE pload, siceload   = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef NONLIN_FRSURF
-CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE recip_hfacc       = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef ANNUAL_BALANCE
-CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE balance_itcount   = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif /* ANNUAL_BALANCE */
 #  ifdef ALLOW_THSICE
-C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
-CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE Qice1, Qice2  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+C-- store thSIce vars before advection (called from SEAICE_MODEL) updates them:
+CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE snowHeight,hOceMxL= comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif /* ALLOW_THSICE */
 # endif /* ALLOW_AUTODIFF_TAMC */
@@ -473,7 +402,20 @@
         CALL SEAICE_MODEL( myTime, myIter, myThid )
         CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
+# ifdef ALLOW_AUTODIFF_TAMC
+CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE heff  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE area  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uIce  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE vIce  = comlev1, key=ikey_dynamics, kind=isbyte
+# endif
 # ifdef ALLOW_COST
         CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
 # endif
+# ifdef ALLOW_AUTODIFF
+      ELSEIF ( SEAICEadjMODE .EQ. -1 ) THEN
+CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
+        CALL SEAICE_FAKE( myTime, myIter, myThid )
+# endif /* ALLOW_AUTODIFF */
       ENDIF
 #endif /* ALLOW_SEAICE */
@@ -493,41 +435,10 @@
 
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qsw                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE sst, sss          = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE qsw               = comlev1, key=ikey_dynamics, kind=isbyte
 # ifdef ALLOW_SEAICE
-CADJ STORE area               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-# endif
-#endif
-
-#ifdef OLD_THSICE_CALL_SEQUENCE
-#if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
-      IF ( useThSIce .AND. fluidIsWater ) THEN
-# ifdef ALLOW_AUTODIFF_TAMC
-cph(
-#  ifdef NONLIN_FRSURF
-CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  endif
-# endif
-# ifdef ALLOW_DEBUG
-        IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
-# endif
-C--     Step forward Therm.Sea-Ice variables
-C       and modify forcing terms including effects from ice
-        CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
-        CALL THSICE_MAIN( myTime, myIter, myThid )
-        CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
-      ENDIF
-#endif /* ALLOW_THSICE */
-#endif /* OLD_THSICE_CALL_SEQUENCE */
+CADJ STORE area              = comlev1, key=ikey_dynamics, kind=isbyte
+# endif
+#endif
 
 #ifdef ALLOW_CPL_ISSM
@@ -541,8 +452,6 @@
 #endif
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel, vvel        = comlev1, key=ikey_dynamics, kind=isbyte
 #endif
 C     compute temperature and (virtual) salt flux at the
@@ -574,5 +483,5 @@
       IF ( useSALT_PLUME ) THEN
 Catn: exchanging saltPlumeFlux:
-          CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
+        CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
       ENDIF
 #endif /* ALLOW_SALT_PLUME */
@@ -581,6 +490,5 @@
       IF ( allowFreezing ) THEN
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE theta             = comlev1, key=ikey_dynamics, kind=isbyte
 #endif
         CALL FREEZE_SURFACE( myTime, myIter, myThid )
@@ -595,6 +503,5 @@
 C     relaxation terms, etc.
 #ifdef ALLOW_AUTODIFF
-CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
 #else  /* ALLOW_AUTODIFF */
 C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
@@ -608,4 +515,16 @@
      I             iMin, iMax, jMin, jMax,
      I             myTime, myIter, myThid )
+
+#ifdef  ALLOW_OBCS
+      IF (useOBCS) THEN
+C--   After all surface fluxes are known apply balancing fluxes and
+C--   apply tidal forcing to open boundaries
+# ifdef ALLOW_DEBUG
+       IF (debugMode) CALL DEBUG_CALL('OBCS_ADJUST',myThid)
+# endif
+       CALL OBCS_ADJUST(
+     I      myTime+deltaTClock, myIter+1, myThid )
+      ENDIF
+#endif  /* ALLOW_OBCS */
 
 #ifdef ALLOW_AUTODIFF_TAMC
@@ -645,9 +564,11 @@
            sigmaY(i,j,k) = 0. _d 0
            sigmaR(i,j,k) = 0. _d 0
+#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) || defined (ALLOW_LEITH_QG))
+           sigmaRfield(i,j,k,bi,bj) = 0. _d 0
+#endif
           ENDDO
          ENDDO
         ENDDO
 
-#ifdef ALLOW_AUTODIFF
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
@@ -656,4 +577,5 @@
          ENDDO
         ENDDO
+#ifdef ALLOW_AUTODIFF
 cph all the following init. are necessary for TAF
 cph although some of these are re-initialised later.
@@ -697,8 +619,6 @@
            Kwy(i,j,k,bi,bj)  = 0. _d 0
            Kwz(i,j,k,bi,bj)  = 0. _d 0
-#  ifdef GM_NON_UNITY_DIAGONAL
            Kux(i,j,k,bi,bj)  = 0. _d 0
            Kvy(i,j,k,bi,bj)  = 0. _d 0
-#  endif
 #  ifdef GM_EXTRA_DIAGONAL
            Kuz(i,j,k,bi,bj)  = 0. _d 0
@@ -736,22 +656,17 @@
 
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE theta(:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE salt (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE totphihyd(:,:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
 # ifdef ALLOW_KPP
-CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 # ifdef ALLOW_SALT_PLUME
-CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE saltplumedepth(:,:,bi,bj)
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE saltplumeflux(:,:,bi,bj)
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 #endif /* ALLOW_AUTODIFF_TAMC */
@@ -856,12 +771,15 @@
             IF (k.GT.1) THEN
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
-CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
-CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
+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 */
-             CALL FIND_RHO_2D(
+             IF ( usingZCoords ) THEN
+              DO j=jMin,jMax
+               DO i=iMin,iMax
+                rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
+               ENDDO
+              ENDDO
+              CALL FIND_RHO_2D(
      I                 iMin, iMax, jMin, jMax, k,
      I                 theta(1-OLx,1-OLy,k-1,bi,bj),
@@ -869,14 +787,21 @@
      O                 rhoKm1,
      I                 k-1, bi, bj, myThid )
+             ELSE
+              CALL FIND_RHO_2D(
+     I                 iMin, iMax, jMin, jMax, k-1,
+     I                 theta(1-OLx,1-OLy,k,bi,bj),
+     I                 salt (1-OLx,1-OLy,k,bi,bj),
+     O                 rhoKp1,
+     I                 k, bi, bj, myThid )
+              DO j=jMin,jMax
+               DO i=iMin,iMax
+                rhoKm1(i,j) = rhoInSitu(i,j,k-1,bi,bj)
+               ENDDO
+              ENDDO
+             ENDIF
             ENDIF
 #ifdef ALLOW_DEBUG
             IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
 #endif
-cph Avoid variable aliasing for adjoint !!!
-            DO j=jMin,jMax
-             DO i=iMin,iMax
-              rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
-             ENDDO
-            ENDDO
             CALL GRAD_SIGMA(
      I             bi, bj, iMin, iMax, jMin, jMax, k,
@@ -884,6 +809,6 @@
      O             sigmaX, sigmaY, sigmaR,
      I             myThid )
-#ifdef ALLOW_ECCO
-# ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
+
+#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) || defined (ALLOW_LEITH_QG))
             DO j=jMin,jMax
              DO i=iMin,iMax
@@ -891,6 +816,6 @@
              ENDDO
             ENDDO
-# endif
-#endif /* ALLOW_ECCO */
+#endif /* ALLOW_SIGMAR_COST_CONTRIBUTION or ALLOW_LEITH_QG */
+
 #ifdef ALLOW_AUTODIFF
 #ifdef GMREDI_WITH_STABLE_ADJOINT
@@ -929,6 +854,5 @@
 #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE IVDConvCount(:,:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif
 
@@ -936,5 +860,5 @@
         IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
           CALL CALC_OCE_MXLAYER(
-     I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
+     I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
      I              bi, bj, myTime, myIter, myThid )
         ENDIF
@@ -943,5 +867,5 @@
         IF ( useSALT_PLUME ) THEN
           CALL SALT_PLUME_CALC_DEPTH(
-     I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
+     I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
      I              bi, bj, myTime, myIter, myThid )
 #ifdef SALT_PLUME_VOLUME
@@ -950,9 +874,9 @@
 C-- get forcings for kpp
           CALL SALT_PLUME_APPLY(
-     I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
+     I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
      I              theta, 0,
      I              myTime, myIter, myThid )
           CALL SALT_PLUME_APPLY(
-     I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
+     I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
      I              salt, 0,
      I              myTime, myIter, myThid )
@@ -978,18 +902,11 @@
 cph needed for KPP
 CADJ STORE surfaceForcingU(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingV(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingS(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingT(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE surfaceForcingTice(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
 
@@ -1047,9 +964,9 @@
 #ifdef  ALLOW_GGL90
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE GGL90TKE(:,:,:,bi,bj)
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
 C--     Compute GGL90 mixing coefficients
-        IF (useGGL90) THEN
+        IF ( useGGL90 .AND. Nr.GT.1 ) THEN
 #ifdef ALLOW_DEBUG
           IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
@@ -1079,10 +996,7 @@
 cph keep it although TAF says you dont need to.
 cph but I have avoided the #ifdef for now, in case more things change
-CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE sigmaX(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE sigmaY(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE sigmaR(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 #endif /* ALLOW_AUTODIFF_TAMC */
@@ -1166,4 +1080,9 @@
 #endif /* ALLOW_KPP */
 
+#ifdef ALLOW_GGL90
+      IF ( useGGL90 )
+     &  CALL GGL90_EXCHANGES( myThid )
+#endif /* ALLOW_GGL90 */
+
 #ifdef ALLOW_DIAGNOSTICS
       IF ( fluidIsWater .AND. useDiagnostics ) THEN
Index: /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F	(revision 26544)
@@ -40,4 +40,10 @@
       INTEGER myComm
 
+C     !FUNCTIONS:
+c     INTEGER  IFNBLNK
+c     EXTERNAL IFNBLNK
+      INTEGER  ILNBLNK
+      EXTERNAL ILNBLNK
+
 C     !LOCAL VARIABLES:
 C     == Local variables ==
@@ -58,5 +64,7 @@
       LOGICAL doReport
 #if defined(ALLOW_OASIS) || defined(COMPONENT_MODULE)
-      INTEGER mpiMyWid
+      INTEGER mpiMyWId
+#elif defined(ALLOW_NEST2W_COMMON)
+      INTEGER mpiMyWId
 #endif
 #ifdef ALLOW_CPL_ISSM
@@ -67,9 +75,13 @@
 #endif /* ALLOW_CPL_ISSM */
 #if defined(ALLOW_NEST_PARENT) || defined(ALLOW_NEST_CHILD)
-      INTEGER mpiMyWid, color
+      INTEGER mpiMyWId, color
 #endif
 #ifdef USE_PDAF
       INTEGER mpi_task_id
+      CHARACTER*(14) fmtStr
+#else
+      CHARACTER*(6) fmtStr
 #endif /* USE_PDAF */
+      INTEGER iTmp
 #endif /* ALLOW_USE_MPI */
 CEOP
@@ -176,4 +188,15 @@
 #endif /* ALLOW_NEST_PARENT | ALLOW_NEST_CHILD */
 
+#ifdef ALLOW_NEST2W_COMMON
+C--    Case with 2-Ways Nest(ing)
+C-     Set the running directory
+       CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
+       CALL SETDIR( mpiMyWId )
+
+C-     Setup Nesting Execution Environment
+       CALL NEST2W_EEINIT( mpiMyWId )
+       IF ( eeBootError ) GOTO 999
+#endif /* ALLOW_NEST2W_COMMON */
+
 #ifdef ALLOW_CPL_ISSM
 C     add a 1rst preliminary call EESET_PARAMS to set useCoupler
@@ -216,9 +239,14 @@
        ENDIF
        myProcId = mpiMyId
+       iTmp = MAX(4,1 + INT(LOG10(DFLOAT(nPx*nPy))))
 #ifdef USE_PDAF
-       WRITE(myProcessStr,'(I4.4,A1,I4.4)') mpi_task_id,'.',myProcId
+       WRITE(fmtStr,'(4(A,I1),A)')
+     &      '(I',iTmp,'.',iTmp,',A1,I',iTmp,'.',iTmp,')'
+       WRITE(myProcessStr,fmtStr) mpi_task_id,'.',myProcId
 #else
-       WRITE(myProcessStr,'(I4.4)') myProcId
+       WRITE(fmtStr,'(2(A,I1),A)') '(I',iTmp,'.',iTmp,')'
+       WRITE(myProcessStr,fmtStr) myProcId
 #endif /* USE_PDAF */
+       iTmp = ILNBLNK( myProcessStr )
        mpiPidIo = myProcId
        pidIO    = mpiPidIo
@@ -227,13 +255,7 @@
         IF( myProcId .EQ. 0 ) THEN
 #endif
-         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:4)
-#ifdef USE_PDAF
-         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:9)
-#endif
+         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:iTmp)
          OPEN(errorMessageUnit,FILE=fNam,STATUS='unknown')
-         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:4)
-#ifdef USE_PDAF
-         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:9)
-#endif
+         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:iTmp)
          OPEN(standardMessageUnit,FILE=fNam,STATUS='unknown')
 #ifdef SINGLE_DISK_IO
Index: /issm/trunk-jpl/test/MITgcm/code/packages.conf
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code/packages.conf	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code/packages.conf	(revision 26544)
@@ -1,5 +1,6 @@
+#-- list of packages (or group of packages) to compile for this experiment:
 gfd
+obcs
+shelfice
+diagnostics
 timeave
-shelfice
-obcs
-diagnostics
Index: /issm/trunk-jpl/test/MITgcm/code_4003/CPP_EEOPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/CPP_EEOPTIONS.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/CPP_EEOPTIONS.h	(revision 26544)
@@ -62,7 +62,4 @@
 C     Use D0 where it means REAL*8 but not where it means REAL*16
 #define D0 d0
-
-C--   Enable some old macro conventions for backward compatibility
-#undef USE_OLD_MACROS_R4R8toRSRL
 
 C=== IO related options ===
Index: /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h	(revision 26544)
@@ -54,4 +54,8 @@
 #define INCLUDE_PHIHYD_CALCULATION_CODE
 
+C o Include/exclude sound speed calculation code
+C o (Note that this is a diagnostic from Del Grasso algorithm, not derived from EOS)
+#undef INCLUDE_SOUNDSPEED_CALC_CODE
+
 C-- Vertical mixing code options:
 
@@ -74,8 +78,12 @@
 #undef EXCLUDE_PCELL_MIX_CODE
 
+C o Exclude/allow to use isotropic 3-D Smagorinsky viscosity as diffusivity
+C   for tracers (after scaling by constant Prandtl number)
+#undef ALLOW_SMAG_3D_DIFFUSIVITY
+
 C-- Time-stepping code options:
 
 C o Include/exclude combined Surf.Pressure and Drag Implicit solver code
-#undef ALLOW_SOLVE4_PS_AND_DRAG
+#define ALLOW_SOLVE4_PS_AND_DRAG
 
 C o Include/exclude Implicit vertical advection code
@@ -84,4 +92,7 @@
 C o Include/exclude AdamsBashforth-3rd-Order code
 #undef ALLOW_ADAMSBASHFORTH_3
+
+C o Include/exclude Quasi-Hydrostatic Stagger Time-step AdamsBashforth code
+#undef ALLOW_QHYD_STAGGER_TS
 
 C-- Model formulation options:
@@ -94,4 +105,7 @@
 C   this implies that grid-cell thickness (hFactors) varies with time
 #define NONLIN_FRSURF
+C o Disable code for rStar coordinate and/or code for Sigma coordinate
+c#define DISABLE_RSTAR_CODE
+c#define DISABLE_SIGMA_CODE
 
 C o Include/exclude nonHydrostatic code
@@ -103,13 +117,13 @@
 C-- Algorithm options:
 
-C o Use Non Self-Adjoint (NSA) conjugate-gradient solver
+C o Include/exclude code for Non Self-Adjoint (NSA) conjugate-gradient solver
 #undef ALLOW_CG2D_NSA
 
 C o Include/exclude code for single reduction Conjugate-Gradient solver
-#undef ALLOW_SRCG
+#define ALLOW_SRCG
 
 C o Choices for implicit solver routines solve_*diagonal.F
 C   The following has low memory footprint, but not suitable for AD
-#undef SOLVE_DIAGONAL_LOWMEMORY
+#define SOLVE_DIAGONAL_LOWMEMORY
 C   The following one suitable for AD but does not vectorize
 #undef SOLVE_DIAGONAL_KINNER
@@ -117,8 +131,39 @@
 C-- Retired code options:
 
+C o ALLOW isotropic scaling of harmonic and bi-harmonic terms when
+C   using an locally isotropic spherical grid with (dlambda) x (dphi*cos(phi))
+C *only for use on a lat-lon grid*
+C   Setting this flag here affects both momentum and tracer equation unless
+C   it is set/unset again in other header fields (e.g., GAD_OPTIONS.h).
+C   The definition of the flag is commented to avoid interference with
+C   such other header files.
+C   The preferred method is specifying a value for viscAhGrid or viscA4Grid
+C   in data which is then automatically scaled by the grid size;
+C   the old method of specifying viscAh/viscA4 and this flag is provided
+C   for completeness only (and for use with the adjoint).
+c#define ISOTROPIC_COS_SCALING
+
+C o This flag selects the form of COSINE(lat) scaling of bi-harmonic term.
+C *only for use on a lat-lon grid*
+C   Has no effect if ISOTROPIC_COS_SCALING is undefined.
+C   Has no effect on vector invariant momentum equations.
+C   Setting this flag here affects both momentum and tracer equation unless
+C   it is set/unset again in other header fields (e.g., GAD_OPTIONS.h).
+C   The definition of the flag is commented to avoid interference with
+C   such other header files.
+c#define COSINEMETH_III
+
+C o Use "OLD" UV discretisation near boundaries (*not* recommended)
+C   Note - only works with pkg/mom_fluxform and "no_slip_sides=.FALSE."
+C          because the old code did not have no-slip BCs
+#undef OLD_ADV_BCS
+
 C o Use LONG.bin, LATG.bin, etc., initialization for ini_curviliear_grid.F
 C   Default is to use "new" grid files (OLD_GRID_IO undef) but OLD_GRID_IO
 C   is still useful with, e.g., single-domain curvilinear configurations.
 #undef OLD_GRID_IO
+
+C o Use old EXTERNAL_FORCING_U,V,T,S subroutines (for backward compatibility)
+#undef USE_OLD_EXTERNAL_FORCING
 
 C-- Other option files:
Index: /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h	(revision 26544)
@@ -23,4 +23,5 @@
       PARAMETER( diagSt_size = 10*Nr )
 
+
 CEH3 ;;; Local Variables: ***
 CEH3 ;;; mode:fortran ***
Index: /issm/trunk-jpl/test/MITgcm/code_4003/DIAG_OPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/DIAG_OPTIONS.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/DIAG_OPTIONS.h	(revision 26544)
@@ -2,13 +2,12 @@
 #define DIAG_OPTIONS_H
 #include "PACKAGES_CONFIG.h"
-#ifdef ALLOW_DIAGNOSTICS
-
 #include "CPP_OPTIONS.h"
 
-CEH3 package-specific options go here
+#ifdef ALLOW_DIAGNOSTICS
+C     Package-specific Options & Macros go here
 
 C allow to define specific regions and the corresponding mask ;
 C  used to perform regional statistics over a limited area
-#undef DIAGSTATS_REGION_MASK
+#define DIAGSTATS_REGION_MASK
 
 C allow to stop & restart at any time (i.e. not at a multiple of
@@ -19,11 +18,6 @@
 #undef  DIAGNOSTICS_HAS_PICKUP
 
-C for NetCDF (mnc) output: define a missing value (default is UNSET_RL)
-C and fill land points of (so far only) scalar fields with it
-#undef DIAGNOSTICS_MISSING_VALUE
-
 #endif /* ALLOW_DIAGNOSTICS */
 #endif /* DIAG_OPTIONS_H */
-
 
 CEH3 ;;; Local Variables: ***
Index: /issm/trunk-jpl/test/MITgcm/code_4003/OBCS_OPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/OBCS_OPTIONS.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/OBCS_OPTIONS.h	(revision 26544)
@@ -1,7 +1,4 @@
 C CPP options file for OBCS package
-C
 C Use this file for selecting options within the OBCS package
-C
-C OBCS is enabled with ALLOW_OBCS in CPP_OPTIONS.h
 
 #ifndef OBCS_OPTIONS_H
@@ -11,4 +8,5 @@
 
 #ifdef ALLOW_OBCS
+C Package-specific Options & Macros go here
 
 C Enable individual open boundaries
@@ -25,11 +23,27 @@
 #define ALLOW_OBCS_PRESCRIBE
 
-C This includes hooks to sponge layer treatment of uvel, vvel
+C Enable OB conditions following Stevens (1990)
+#undef ALLOW_OBCS_STEVENS
+
+C Allow sponge layer treatment of open boundary conditions
 #define ALLOW_OBCS_SPONGE
 
+C Include hooks to sponge layer treatment of pkg/seaice variables
+#undef ALLOW_OBCS_SEAICE_SPONGE
+
 C balance barotropic velocity
-C KS16 turn this on
 #define ALLOW_OBCS_BALANCE
 
+C add tidal contributions to normal OB flow
+C NOTE that at the moment tidal forcing is applied
+C only to "normal" flow.
+#undef ALLOW_OBCS_TIDES
+
+C Use older implementation of obcs in seaice-dynamics
+C note: most of the "experimental" options listed below have not yet
+C       been implementated in new version.
+#undef OBCS_UVICE_OLD
+
+#ifdef OBCS_UVICE_OLD
 C     The following five CPP options are experimental and aim to deal
 C     with artifacts due to the low-frequency specification of sea-ice
@@ -46,10 +60,11 @@
 #undef OBCS_SEAICE_SMOOTH_UVICE_PAR
 
+C     Compute rather than specify seaice velocities at the edges.
+#undef OBCS_SEAICE_COMPUTE_UVICE
+#endif /* OBCS_UVICE_OLD */
+
 C     Smooth the tracer sea-ice variables near the edges.
 #undef OBCS_SEAICE_SMOOTH_EDGE
 
-C     Compute rather than specify seaice velocities at the edges.
-#undef OBCS_SEAICE_COMPUTE_UVICE
-
 #endif /* ALLOW_OBCS */
 #endif /* OBCS_OPTIONS_H */
Index: /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h	(revision 26544)
@@ -43,10 +43,10 @@
       INTEGER Nr
       PARAMETER (
-     &           sNx =  3,
-     &           sNy =  200,
+     &           sNx =   1,
+     &           sNy =  50,
      &           OLx =   3,
      &           OLy =   3,
      &           nSx =   1,
-     &           nSy =   1,
+     &           nSy =   4,
      &           nPx =   1,
      &           nPy =   1,
Index: sm/trunk-jpl/test/MITgcm/code_4003/STREAMICE_OPTIONS.h
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/STREAMICE_OPTIONS.h	(revision 26543)
+++ 	(revision )
@@ -1,32 +1,0 @@
-C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-
-C CPP options file for MYPACKAGE
-C
-C Use this file for selecting options within package "streamice"
-
-#ifndef STREAMICE_OPTIONS_H
-#define STREAMICE_OPTIONS_H
-#include "PACKAGES_CONFIG.h"
-#ifdef ALLOW_STREAMICE
-
-#include "CPP_OPTIONS.h"
-
-C Place CPP define/undef flag here
-
-#define STREAMICE_CONSTRUCT_MATRIX
-#define STREAMICE_HYBRID_STRESS
-#define USE_ALT_RLOW
-#define STREAMICE_GEOM_FILE_SETUP
-!#define STREAMICE_SMOOTH_FLOATATION
-!#define ALLOW_PETSC
-#undef ALLOW_STREAMICE_2DTRACER
-#undef STREAMICE_TRACER_AB
-#undef STREAMICE_SERIAL_TRISOLVE
-#define STREAMICE_FLOWLINE_BUTTRESS
-
-#endif /* ALLOW_MYPACKAGE */
-#endif /* MYPACKAGE_OPTIONS_H */
-
-CEH3 ;;; Local Variables: ***
-CEH3 ;;; mode:fortran ***
-CEH3 ;;; End: ***
Index: /issm/trunk-jpl/test/MITgcm/code_4003/do_oceanic_phys.F
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/do_oceanic_phys.F	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/do_oceanic_phys.F	(revision 26544)
@@ -1,4 +1,7 @@
 #include "PACKAGES_CONFIG.h"
 #include "CPP_OPTIONS.h"
+#ifdef ALLOW_MOM_COMMON
+# include "MOM_COMMON_OPTIONS.h"
+#endif
 #ifdef ALLOW_AUTODIFF
 # include "AUTODIFF_OPTIONS.h"
@@ -73,4 +76,6 @@
 C       |-- EXTERNAL_FORCING_SURF
 C       |
+C       |-- OBCS_ADJUST
+C       |
 C       |- k loop (Nr:1):
 C       | - DWNSLP_CALC_RHO
@@ -119,4 +124,6 @@
 C       |
 C       |-- KPP_DO_EXCH
+C       |
+C       |-- GGL90_EXCHANGES
 C       |
 C       |-- DIAGS_RHO_G
@@ -143,6 +150,7 @@
 #ifdef ALLOW_AUTODIFF
 # include "AUTODIFF_MYFIELDS.h"
-# include "tamc.h"
-# include "tamc_keys.h"
+# ifdef ALLOW_AUTODIFF_TAMC
+#  include "tamc.h"
+# endif
 # include "FFIELDS.h"
 # include "SURFACE.h"
@@ -196,5 +204,5 @@
 C     !LOCAL VARIABLES:
 C     == Local variables
-C     rhoK, rhoKm1  :: Density at current level, and level above
+C     rhoKp1,rhoKm1 :: Density at current level, and @ level minus one
 C     iMin, iMax    :: Ranges and sub-block indices on which calculations
 C     jMin, jMax       are applied.
@@ -202,4 +210,5 @@
 C     msgBuf        :: Temp. for building output string
 C     i,j,k         :: loop indices
+C     kSrf          :: surface index
       _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
@@ -210,5 +219,5 @@
       INTEGER jMin, jMax
       INTEGER bi, bj
-      INTEGER i, j, k
+      INTEGER i, j, k, kSrf
       CHARACTER*(MAX_LEN_MBUF) msgBuf
       INTEGER doDiagsRho
@@ -221,4 +230,9 @@
       _RL thetaRef
 #endif /* ALLOW_AUTODIFF */
+#ifdef ALLOW_AUTODIFF_TAMC
+      INTEGER act1, act2, act3, act4
+      INTEGER max1, max2, max3
+      INTEGER kkey, itdkey
+#endif
 CEOP
 
@@ -227,4 +241,7 @@
       itdkey = 1
 #endif /* ALLOW_AUTODIFF_TAMC */
+
+      kSrf = 1
+      IF ( usingPCoords ) kSrf = Nr
 
 #ifdef ALLOW_DEBUG
@@ -282,31 +299,17 @@
 
 #ifdef ALLOW_AUTODIFF
-# ifdef ALLOW_SALT_PLUME
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
+          adjustColdSST_diag(i,j,bi,bj) = 0. _d 0
+# ifdef ALLOW_SALT_PLUME
           saltPlumeDepth(i,j,bi,bj) = 0. _d 0
           saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
+# endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
-# endif
-# ifdef ALLOW_ECCO
-#  ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
-      DO bj=myByLo(myThid),myByHi(myThid)
-       DO bi=myBxLo(myThid),myBxHi(myThid)
-        DO k=1,Nr
-         DO j=1-OLy,sNy+OLy
-          DO i=1-OLx,sNx+OLx
-           sigmaRfield(i,j,k,bi,bj) = 0. _d 0
-          ENDDO
-         ENDDO
-        ENDDO
-       ENDDO
-      ENDDO
-#  endif
-# endif
 #endif /* ALLOW_AUTODIFF */
 
@@ -320,35 +323,22 @@
 #endif /* ALLOW_FRAZIL */
 
-#ifndef OLD_THSICE_CALL_SEQUENCE
 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
       IF ( useThSIce .AND. fluidIsWater ) THEN
 # ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE snowHeight, Tsrf  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE Qice1, Qice2      = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE salt,theta        = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uvel,vvel         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qnet,qsw, empmr   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE atemp,aqh,precip  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE swdown,lwdown     = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+#  ifdef ALLOW_SEAICE
+CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
+#  endif
+CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE snowHeight, Tsrf  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE sHeating,snowAge  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE hocemxl, icflxsw  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE salt,theta        = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE qnet,qsw, empmr   = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
 #  ifdef NONLIN_FRSURF
-CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE hFac_surfC        = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 # endif /* ALLOW_AUTODIFF_TAMC */
@@ -363,106 +353,45 @@
       ENDIF
 #endif /* ALLOW_THSICE */
-#endif /* ndef OLD_THSICE_CALL_SEQUENCE */
 
 #ifdef ALLOW_SEAICE
-# ifdef ALLOW_AUTODIFF
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
-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
+# 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)
-CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
-#endif
-      IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
-        CALL SEAICE_FAKE( myTime, myIter, myThid )
-      ENDIF
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
-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
-#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
-CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
-#endif
-# endif /* ALLOW_AUTODIFF */
-#endif /* ALLOW_SEAICE */
-
-#ifdef ALLOW_SEAICE
+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
       IF ( useSEAICE ) THEN
 # ifdef ALLOW_AUTODIFF_TAMC
-cph-adj-test(
-CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE tices  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE empmr, qnet  = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
-CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
-cCADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
-cCADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
-cph-adj-test)
-c#ifdef ALLOW_EXF
-CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE evap                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uwind,vwind         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-c#endif
-CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  ifdef SEAICE_CGRID
-CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  endif
-#  ifdef SEAICE_ALLOW_DYNAMICS
-CADJ STORE uice                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE vice                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE dwatn               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#   ifdef SEAICE_ALLOW_EVP
-CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#   endif
+CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
+#  ifdef ALLOW_EXF
+CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uwind,vwind       = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef SEAICE_VARIABLE_SALINITY
-CADJ STORE hsalt               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE hsalt             = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef ATMOSPHERIC_LOADING
-CADJ STORE pload               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE siceload            = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE pload, siceload   = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef NONLIN_FRSURF
-CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE recip_hfacc       = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif
 #  ifdef ANNUAL_BALANCE
-CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE balance_itcount   = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif /* ANNUAL_BALANCE */
 #  ifdef ALLOW_THSICE
-C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
-CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE Qice1, Qice2  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+C-- store thSIce vars before advection (called from SEAICE_MODEL) updates them:
+CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE snowHeight,hOceMxL= comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
 #  endif /* ALLOW_THSICE */
 # endif /* ALLOW_AUTODIFF_TAMC */
@@ -473,7 +402,20 @@
         CALL SEAICE_MODEL( myTime, myIter, myThid )
         CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
+# ifdef ALLOW_AUTODIFF_TAMC
+CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE heff  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE area  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uIce  = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE vIce  = comlev1, key=ikey_dynamics, kind=isbyte
+# endif
 # ifdef ALLOW_COST
         CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
 # endif
+# ifdef ALLOW_AUTODIFF
+      ELSEIF ( SEAICEadjMODE .EQ. -1 ) THEN
+CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
+        CALL SEAICE_FAKE( myTime, myIter, myThid )
+# endif /* ALLOW_AUTODIFF */
       ENDIF
 #endif /* ALLOW_SEAICE */
@@ -493,41 +435,10 @@
 
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qsw                = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE sst, sss          = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE qsw               = comlev1, key=ikey_dynamics, kind=isbyte
 # ifdef ALLOW_SEAICE
-CADJ STORE area               = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-# endif
-#endif
-
-#ifdef OLD_THSICE_CALL_SEQUENCE
-#if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
-      IF ( useThSIce .AND. fluidIsWater ) THEN
-# ifdef ALLOW_AUTODIFF_TAMC
-cph(
-#  ifdef NONLIN_FRSURF
-CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-#  endif
-# endif
-# ifdef ALLOW_DEBUG
-        IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
-# endif
-C--     Step forward Therm.Sea-Ice variables
-C       and modify forcing terms including effects from ice
-        CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
-        CALL THSICE_MAIN( myTime, myIter, myThid )
-        CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
-      ENDIF
-#endif /* ALLOW_THSICE */
-#endif /* OLD_THSICE_CALL_SEQUENCE */
+CADJ STORE area              = comlev1, key=ikey_dynamics, kind=isbyte
+# endif
+#endif
 
 #ifdef ALLOW_CPL_ISSM
@@ -541,8 +452,6 @@
 #endif
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
-CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
+CADJ STORE uvel, vvel        = comlev1, key=ikey_dynamics, kind=isbyte
 #endif
 C     compute temperature and (virtual) salt flux at the
@@ -574,5 +483,5 @@
       IF ( useSALT_PLUME ) THEN
 Catn: exchanging saltPlumeFlux:
-          CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
+        CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
       ENDIF
 #endif /* ALLOW_SALT_PLUME */
@@ -581,6 +490,5 @@
       IF ( allowFreezing ) THEN
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE theta             = comlev1, key=ikey_dynamics, kind=isbyte
 #endif
         CALL FREEZE_SURFACE( myTime, myIter, myThid )
@@ -595,6 +503,5 @@
 C     relaxation terms, etc.
 #ifdef ALLOW_AUTODIFF
-CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
-CADJ &     kind = isbyte
+CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
 #else  /* ALLOW_AUTODIFF */
 C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
@@ -608,4 +515,16 @@
      I             iMin, iMax, jMin, jMax,
      I             myTime, myIter, myThid )
+
+#ifdef  ALLOW_OBCS
+      IF (useOBCS) THEN
+C--   After all surface fluxes are known apply balancing fluxes and
+C--   apply tidal forcing to open boundaries
+# ifdef ALLOW_DEBUG
+       IF (debugMode) CALL DEBUG_CALL('OBCS_ADJUST',myThid)
+# endif
+       CALL OBCS_ADJUST(
+     I      myTime+deltaTClock, myIter+1, myThid )
+      ENDIF
+#endif  /* ALLOW_OBCS */
 
 #ifdef ALLOW_AUTODIFF_TAMC
@@ -645,9 +564,11 @@
            sigmaY(i,j,k) = 0. _d 0
            sigmaR(i,j,k) = 0. _d 0
+#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) || defined (ALLOW_LEITH_QG))
+           sigmaRfield(i,j,k,bi,bj) = 0. _d 0
+#endif
           ENDDO
          ENDDO
         ENDDO
 
-#ifdef ALLOW_AUTODIFF
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
@@ -656,4 +577,5 @@
          ENDDO
         ENDDO
+#ifdef ALLOW_AUTODIFF
 cph all the following init. are necessary for TAF
 cph although some of these are re-initialised later.
@@ -697,8 +619,6 @@
            Kwy(i,j,k,bi,bj)  = 0. _d 0
            Kwz(i,j,k,bi,bj)  = 0. _d 0
-#  ifdef GM_NON_UNITY_DIAGONAL
            Kux(i,j,k,bi,bj)  = 0. _d 0
            Kvy(i,j,k,bi,bj)  = 0. _d 0
-#  endif
 #  ifdef GM_EXTRA_DIAGONAL
            Kuz(i,j,k,bi,bj)  = 0. _d 0
@@ -736,22 +656,17 @@
 
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE theta(:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE salt (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE totphihyd(:,:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
 # ifdef ALLOW_KPP
-CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 # ifdef ALLOW_SALT_PLUME
-CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE saltplumedepth(:,:,bi,bj)
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE saltplumeflux(:,:,bi,bj)
+CADJ &                         = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 #endif /* ALLOW_AUTODIFF_TAMC */
@@ -856,12 +771,15 @@
             IF (k.GT.1) THEN
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
-CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
-CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
-CADJ &     kind = isbyte
+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 */
-             CALL FIND_RHO_2D(
+             IF ( usingZCoords ) THEN
+              DO j=jMin,jMax
+               DO i=iMin,iMax
+                rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
+               ENDDO
+              ENDDO
+              CALL FIND_RHO_2D(
      I                 iMin, iMax, jMin, jMax, k,
      I                 theta(1-OLx,1-OLy,k-1,bi,bj),
@@ -869,14 +787,21 @@
      O                 rhoKm1,
      I                 k-1, bi, bj, myThid )
+             ELSE
+              CALL FIND_RHO_2D(
+     I                 iMin, iMax, jMin, jMax, k-1,
+     I                 theta(1-OLx,1-OLy,k,bi,bj),
+     I                 salt (1-OLx,1-OLy,k,bi,bj),
+     O                 rhoKp1,
+     I                 k, bi, bj, myThid )
+              DO j=jMin,jMax
+               DO i=iMin,iMax
+                rhoKm1(i,j) = rhoInSitu(i,j,k-1,bi,bj)
+               ENDDO
+              ENDDO
+             ENDIF
             ENDIF
 #ifdef ALLOW_DEBUG
             IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
 #endif
-cph Avoid variable aliasing for adjoint !!!
-            DO j=jMin,jMax
-             DO i=iMin,iMax
-              rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
-             ENDDO
-            ENDDO
             CALL GRAD_SIGMA(
      I             bi, bj, iMin, iMax, jMin, jMax, k,
@@ -884,6 +809,6 @@
      O             sigmaX, sigmaY, sigmaR,
      I             myThid )
-#ifdef ALLOW_ECCO
-# ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
+
+#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) || defined (ALLOW_LEITH_QG))
             DO j=jMin,jMax
              DO i=iMin,iMax
@@ -891,6 +816,6 @@
              ENDDO
             ENDDO
-# endif
-#endif /* ALLOW_ECCO */
+#endif /* ALLOW_SIGMAR_COST_CONTRIBUTION or ALLOW_LEITH_QG */
+
 #ifdef ALLOW_AUTODIFF
 #ifdef GMREDI_WITH_STABLE_ADJOINT
@@ -929,6 +854,5 @@
 #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE IVDConvCount(:,:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif
 
@@ -936,5 +860,5 @@
         IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
           CALL CALC_OCE_MXLAYER(
-     I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
+     I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
      I              bi, bj, myTime, myIter, myThid )
         ENDIF
@@ -943,5 +867,5 @@
         IF ( useSALT_PLUME ) THEN
           CALL SALT_PLUME_CALC_DEPTH(
-     I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
+     I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
      I              bi, bj, myTime, myIter, myThid )
 #ifdef SALT_PLUME_VOLUME
@@ -950,9 +874,9 @@
 C-- get forcings for kpp
           CALL SALT_PLUME_APPLY(
-     I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
+     I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
      I              theta, 0,
      I              myTime, myIter, myThid )
           CALL SALT_PLUME_APPLY(
-     I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
+     I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
      I              salt, 0,
      I              myTime, myIter, myThid )
@@ -978,18 +902,11 @@
 cph needed for KPP
 CADJ STORE surfaceForcingU(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingV(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingS(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 CADJ STORE surfaceForcingT(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE surfaceForcingTice(:,:,bi,bj)
-CADJ &     = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
 
@@ -1047,9 +964,9 @@
 #ifdef  ALLOW_GGL90
 #ifdef ALLOW_AUTODIFF_TAMC
-CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE GGL90TKE(:,:,:,bi,bj)
+CADJ &     = comlev1_bibj, key=itdkey, kind=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
 C--     Compute GGL90 mixing coefficients
-        IF (useGGL90) THEN
+        IF ( useGGL90 .AND. Nr.GT.1 ) THEN
 #ifdef ALLOW_DEBUG
           IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
@@ -1079,10 +996,7 @@
 cph keep it although TAF says you dont need to.
 cph but I have avoided the #ifdef for now, in case more things change
-CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
-CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
-CADJ &     kind = isbyte
+CADJ STORE sigmaX(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE sigmaY(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
+CADJ STORE sigmaR(:,:,:)       = comlev1_bibj, key=itdkey, kind=isbyte
 # endif
 #endif /* ALLOW_AUTODIFF_TAMC */
@@ -1166,4 +1080,9 @@
 #endif /* ALLOW_KPP */
 
+#ifdef ALLOW_GGL90
+      IF ( useGGL90 )
+     &  CALL GGL90_EXCHANGES( myThid )
+#endif /* ALLOW_GGL90 */
+
 #ifdef ALLOW_DIAGNOSTICS
       IF ( fluidIsWater .AND. useDiagnostics ) THEN
Index: /issm/trunk-jpl/test/MITgcm/code_4003/eeboot_minimal.F
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/eeboot_minimal.F	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/eeboot_minimal.F	(revision 26544)
@@ -40,4 +40,10 @@
       INTEGER myComm
 
+C     !FUNCTIONS:
+c     INTEGER  IFNBLNK
+c     EXTERNAL IFNBLNK
+      INTEGER  ILNBLNK
+      EXTERNAL ILNBLNK
+
 C     !LOCAL VARIABLES:
 C     == Local variables ==
@@ -58,5 +64,7 @@
       LOGICAL doReport
 #if defined(ALLOW_OASIS) || defined(COMPONENT_MODULE)
-      INTEGER mpiMyWid
+      INTEGER mpiMyWId
+#elif defined(ALLOW_NEST2W_COMMON)
+      INTEGER mpiMyWId
 #endif
 #ifdef ALLOW_CPL_ISSM
@@ -67,9 +75,13 @@
 #endif /* ALLOW_CPL_ISSM */
 #if defined(ALLOW_NEST_PARENT) || defined(ALLOW_NEST_CHILD)
-      INTEGER mpiMyWid, color
+      INTEGER mpiMyWId, color
 #endif
 #ifdef USE_PDAF
       INTEGER mpi_task_id
+      CHARACTER*(14) fmtStr
+#else
+      CHARACTER*(6) fmtStr
 #endif /* USE_PDAF */
+      INTEGER iTmp
 #endif /* ALLOW_USE_MPI */
 CEOP
@@ -176,4 +188,15 @@
 #endif /* ALLOW_NEST_PARENT | ALLOW_NEST_CHILD */
 
+#ifdef ALLOW_NEST2W_COMMON
+C--    Case with 2-Ways Nest(ing)
+C-     Set the running directory
+       CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
+       CALL SETDIR( mpiMyWId )
+
+C-     Setup Nesting Execution Environment
+       CALL NEST2W_EEINIT( mpiMyWId )
+       IF ( eeBootError ) GOTO 999
+#endif /* ALLOW_NEST2W_COMMON */
+
 #ifdef ALLOW_CPL_ISSM
 C     add a 1rst preliminary call EESET_PARAMS to set useCoupler
@@ -216,9 +239,14 @@
        ENDIF
        myProcId = mpiMyId
+       iTmp = MAX(4,1 + INT(LOG10(DFLOAT(nPx*nPy))))
 #ifdef USE_PDAF
-       WRITE(myProcessStr,'(I4.4,A1,I4.4)') mpi_task_id,'.',myProcId
+       WRITE(fmtStr,'(4(A,I1),A)')
+     &      '(I',iTmp,'.',iTmp,',A1,I',iTmp,'.',iTmp,')'
+       WRITE(myProcessStr,fmtStr) mpi_task_id,'.',myProcId
 #else
-       WRITE(myProcessStr,'(I4.4)') myProcId
+       WRITE(fmtStr,'(2(A,I1),A)') '(I',iTmp,'.',iTmp,')'
+       WRITE(myProcessStr,fmtStr) myProcId
 #endif /* USE_PDAF */
+       iTmp = ILNBLNK( myProcessStr )
        mpiPidIo = myProcId
        pidIO    = mpiPidIo
@@ -227,13 +255,7 @@
         IF( myProcId .EQ. 0 ) THEN
 #endif
-         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:4)
-#ifdef USE_PDAF
-         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:9)
-#endif
+         WRITE(fNam,'(A,A)') 'STDERR.', myProcessStr(1:iTmp)
          OPEN(errorMessageUnit,FILE=fNam,STATUS='unknown')
-         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:4)
-#ifdef USE_PDAF
-         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:9)
-#endif
+         WRITE(fNam,'(A,A)') 'STDOUT.', myProcessStr(1:iTmp)
          OPEN(standardMessageUnit,FILE=fNam,STATUS='unknown')
 #ifdef SINGLE_DISK_IO
Index: /issm/trunk-jpl/test/MITgcm/code_4003/packages.conf
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/packages.conf	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/packages.conf	(revision 26544)
@@ -1,4 +1,5 @@
+#-- list of packages (or group of packages) to compile for this experiment:
 gfd
+obcs
 shelfice
-obcs
 diagnostics
Index: /issm/trunk-jpl/test/MITgcm/code_4003/shelfice_step_icemass.F
===================================================================
--- /issm/trunk-jpl/test/MITgcm/code_4003/shelfice_step_icemass.F	(revision 26543)
+++ /issm/trunk-jpl/test/MITgcm/code_4003/shelfice_step_icemass.F	(revision 26544)
@@ -45,6 +45,8 @@
 C     i,j, bi,bj  :: loop indices
       INTEGER bi,bj,i,j
+#ifdef ALLOW_CPL_ISSM
       _RL  issm_density
       issm_density = 917.
+#endif /* ALLOW_CPL_ISSM */
 
       IF ( SHELFICEMassStepping ) THEN