--- MITgcm/model/src/convective_adjustment.F	2002/09/18 16:38:01	1.14
+++ MITgcm/model/src/convective_adjustment.F	2002/11/01 22:00:33	1.15
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/convective_adjustment.F,v 1.14 2002/09/18 16:38:01 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/convective_adjustment.F,v 1.15 2002/11/01 22:00:33 mlosch Exp $
 C $Name:  $
 
 #include "CPP_OPTIONS.h"
@@ -54,7 +54,7 @@
 C     == Local variables ==
 C     rhoKm1, rhoK :: Density at adjacent levels (common ref. level)
 C     ConvectCount :: Convection mixing freq. counter.
-      INTEGER i, j, K
+      INTEGER i, j, K, kTop, kBottom, kDir, kRef
       _RL rhoKm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhoK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL ConvectCount(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
@@ -94,8 +94,32 @@
      &     ) THEN
 #endif  /* ALLOW_KPP */
 
+         IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
+          kTop    = 2
+          kBottom = Nr
+          kDir    = 1
+         ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
+          kTop    = Nr-1
+          kBottom =  1
+          kDir    = -1
+         ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
+          kTop    = Nr-1
+          kBottom =  1
+          kDir    = -1
+         ELSE
+          STOP 'CONVECTIVE_ADJUSTMENT: We should never reach this point'
+         ENDIF	
+
 C--       Loop over all *interior* layers
-          DO K=2,Nr
+          DO K=kTop,kBottom,kDir
+
+           IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
+            kRef = k-1
+           ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
+            kRef = k
+           ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
+            kRef = k
+           ENDIF	
 
 #ifdef ALLOW_AUTODIFF_TAMC
             kkey = (ikey-1)*Nr + k
@@ -104,7 +128,7 @@
 #endif /* ALLOW_AUTODIFF_TAMC */
 C-          Density of K-1 layer (above W(K)) reference to K-1 T-level
             CALL FIND_RHO(
-     I           bi,bj,iMin,iMax,jMin,jMax,K-1,K-1,
+     I           bi,bj,iMin,iMax,jMin,jMax,K-1,KRef,
      I           theta,salt,
      O           rhoKm1,
      I           myThid )
@@ -115,7 +139,7 @@
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
             CALL FIND_RHO(
-     I           bi,bj,iMin,iMax,jMin,jMax,K,K-1,
+     I           bi,bj,iMin,iMax,jMin,jMax,K,KRef,
      I           theta,salt,
      O           rhoK,
      I           myThid )