--- MITgcm/pkg/pp81/pp81_ri_number.F	2007/08/23 19:13:10	1.2
+++ MITgcm/pkg/pp81/pp81_ri_number.F	2008/08/11 22:28:06	1.3
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/pp81/pp81_ri_number.F,v 1.2 2007/08/23 19:13:10 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/pp81/pp81_ri_number.F,v 1.3 2008/08/11 22:28:06 jmc Exp $
 C $Name:  $
 
 #include "PP81_OPTIONS.h"
@@ -8,7 +8,7 @@
 
 C !INTERFACE: ==========================================================
       subroutine PP81_RI_NUMBER(
-     I     bi, bj, K, iMin, iMax, jMin, jMax, 
+     I     bi, bj, K, iMin, iMax, jMin, jMax,
      O     RiNumber,
      I     myTime, myThid )
 
@@ -41,7 +41,7 @@
 C !INPUT PARAMETERS: ===================================================
 C Routine arguments
 C     bi, bj - array indices on which to apply calculations
-C     iMin, iMax, jMin, jMax 
+C     iMin, iMax, jMin, jMax
 C            - array boundaries
 C     k      - depth level
 C     myTime - Current time in simulation
@@ -59,22 +59,22 @@
 C     tempu, tempv - temporary variables
 C     rhoK, rhoKm1 - Density below and above current interface
 C     epsilon      - small number
-C     RiFlux       - denominator of Richardson number 
+C     RiFlux       - denominator of Richardson number
 C     BuoyFreq     - buoyancy frequency
       INTEGER I,J,Km1
-      _RL        p5    , p125      
+      _RL        p5    , p125
       PARAMETER( p5=0.5, p125=0.125 )
       _RL tempu, tempv
       _RL epsilon
       PARAMETER    (  epsilon = 1.D-10 )
-      _RL RiFlux  
+      _RL RiFlux
       _RL buoyFreq
       _RL rhoKm1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhoK     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 #ifdef PP81_SMOOTH_RI
       _RL RiTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 #endif /* PP81_SMOOTH_RI */
-CEOP     
+CEOP
 
       Km1   = MAX(1,K-1)
 C     Not sure what is correct for pressure coordinates:
@@ -92,16 +92,16 @@
 CML       STOP 'PP81_RI_NUMBER: We should never reach this point'
 CML      ENDIF
 C     preparation: find density a kth and k-1st level
-      CALL FIND_RHO(
-     I     bi, bj, iMin, iMax, jMin, jMax, Km1, K,
-     I     theta, salt,
+      CALL FIND_RHO_2D(
+     I     iMin, iMax, jMin, jMax, K,
+     I     theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
      O     rhoKm1,
-     I     myThid )
-      CALL FIND_RHO(
-     I     bi, bj, iMin, iMax, jMin, jMax, K, K,
-     I     theta, salt,
+     I     Km1, bi, bj, myThid )
+      CALL FIND_RHO_2D(
+     I     iMin, iMax, jMin, jMax, K,
+     I     theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
      O     rhoK,
-     I     myThid )
+     I     K, bi, bj, myThid )
 
 C     first step:  calculate flux Richardson number.
 C     calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
@@ -120,28 +120,28 @@
 C
         buoyFreq = - gravity*mass2rUnit *
      &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
-C     
+C
 C     calculates gradient Richardson number, bounded by
 C     a very large negative value.
-C     
+C
         RiNumber(I,J) = buoyFreq/max(RiFlux,epsilon)
        ENDDO
       ENDDO
-      
+
 #ifdef PP81_SMOOTH_RI
 C     average Richardson number horizontally
       DO J=jMin,jMax
        DO I=iMin,iMax
         RiTmp(I,J) = RiNumber(I,J)
        ENDDO
-      ENDDO	
+      ENDDO
       DO J=1-Oly+1,sNy+Oly-1
        DO I=1-Olx+1,sNx+Olx-1
-        RiNumber(I,J) = p5*RiTmp(I,J) 
-     &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J) 
-     &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1) 
+        RiNumber(I,J) = p5*RiTmp(I,J)
+     &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
+     &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
        ENDDO
-      ENDDO	
+      ENDDO
 #endif /* PP81_SMOOTH_RI */
 
 #endif /* ALLOW_PP81 */