pkg/pp81/pp81_ri_number.F

C $Header: /u/gcmpack/MITgcm/pkg/pp81/pp81_ri_number.F,v 1.2 2007/08/23 19:13:10 jmc Exp $
C $Name:  $

#include "PP81_OPTIONS.h"

CBOP
C !ROUTINE: PP81_RI_NUMBER

C !INTERFACE: ==========================================================
      subroutine PP81_RI_NUMBER(
     I     bi, bj, K, iMin, iMax, jMin, jMax,
     O     RiNumber,
     I     myTime, myThid )

C !DESCRIPTION: \bv
C     /==========================================================\
C     | SUBROUTINE PP81_RI_NUMBER                                |
C     | o Compute gradient Richardson number for Pacanowski and  |
C     |   Philander (1981) mixing scheme                         |
C     \==========================================================/
      IMPLICIT NONE
c
c-------------------------------------------------------------------------

c \ev

C !USES: ===============================================================
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "PP81.h"
#include "GRID.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#else /* ALLOW_AUTODIFF_TAMC */
      integer ikppkey
#endif /* ALLOW_AUTODIFF_TAMC */

C !INPUT PARAMETERS: ===================================================
C Routine arguments
C     bi, bj - array indices on which to apply calculations
C     iMin, iMax, jMin, jMax
C            - array boundaries
C     k      - depth level
C     myTime - Current time in simulation
C     RiNumber - (output) Richardson number
      INTEGER bi, bj, iMin, iMax, jMin, jMax, k
      INTEGER myThid
      _RL     myTime
      _RL     RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

#ifdef ALLOW_PP81

C !LOCAL VARIABLES: ====================================================
C     I,J,kUp      - loop indices
C     p0-125       - averaging coefficients
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     BuoyFreq     - buoyancy frequency
      INTEGER I,J,Km1
      _RL        p5    , p125
      PARAMETER( p5=0.5, p125=0.125 )
      _RL tempu, tempv
      _RL epsilon
      PARAMETER    (  epsilon = 1.D-10 )
      _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

      Km1   = MAX(1,K-1)
C     Not sure what is correct for pressure coordinates:
C     RiFlux is always correct because it quadratic
C     buoyFreq should also be correct in pressure coordinates:
C     N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
C     So the following is wrong:
CML      IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
CML       kUp   = MAX(1,K-1)
CML       kDown = K
CML      ELSEIF  ( buoyancyRelation .EQ. 'OCEANIP') THEN
CML       kUp   = K
CML       kDown = MAX(1,K-1)
CML      ELSE
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_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     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     K, bi, bj, myThid )

C     first step:  calculate flux Richardson number.
C     calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
      DO J= jMin, jMax
       DO I = iMin, iMax
        tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
     &            - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
     &       *recip_drC(K)
        tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
     &            - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
     &       *recip_drC(K)
        RiFlux = tempu*tempu+tempv*tempv

C
C     calculate - g*(drho/dz)/rho0= N^2  .
C
        buoyFreq = - gravity*mass2rUnit *
     &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
C
C     calculates gradient Richardson number, bounded by
C     a very large negative value.
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
      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)
       ENDDO
      ENDDO
#endif /* PP81_SMOOTH_RI */

#endif /* ALLOW_PP81 */

      RETURN
      END

1	jmc	1.3	C $Header: /u/gcmpack/MITgcm/pkg/pp81/pp81_ri_number.F,v 1.2 2007/08/23 19:13:10 jmc Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "PP81_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: PP81_RI_NUMBER
8
9			C !INTERFACE: ==========================================================
10			subroutine PP81_RI_NUMBER(
11	jmc	1.3	I bi, bj, K, iMin, iMax, jMin, jMax,
12	mlosch	1.1	O RiNumber,
13			I myTime, myThid )
14
15			C !DESCRIPTION: \bv
16			C /==========================================================\
17			C \| SUBROUTINE PP81_RI_NUMBER \|
18			C \| o Compute gradient Richardson number for Pacanowski and \|
19			C \| Philander (1981) mixing scheme \|
20			C \==========================================================/
21			IMPLICIT NONE
22			c
23			c-------------------------------------------------------------------------
24
25			c \ev
26
27			C !USES: ===============================================================
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "DYNVARS.h"
32			#include "PP81.h"
33			#include "GRID.h"
34			#ifdef ALLOW_AUTODIFF_TAMC
35			#include "tamc.h"
36			#include "tamc_keys.h"
37			#else /* ALLOW_AUTODIFF_TAMC */
38			integer ikppkey
39			#endif /* ALLOW_AUTODIFF_TAMC */
40
41			C !INPUT PARAMETERS: ===================================================
42			C Routine arguments
43			C bi, bj - array indices on which to apply calculations
44	jmc	1.3	C iMin, iMax, jMin, jMax
45	mlosch	1.1	C - array boundaries
46			C k - depth level
47			C myTime - Current time in simulation
48			C RiNumber - (output) Richardson number
49			INTEGER bi, bj, iMin, iMax, jMin, jMax, k
50			INTEGER myThid
51			_RL myTime
52			_RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53
54			#ifdef ALLOW_PP81
55
56			C !LOCAL VARIABLES: ====================================================
57			C I,J,kUp - loop indices
58			C p0-125 - averaging coefficients
59			C tempu, tempv - temporary variables
60			C rhoK, rhoKm1 - Density below and above current interface
61			C epsilon - small number
62	jmc	1.3	C RiFlux - denominator of Richardson number
63	mlosch	1.1	C BuoyFreq - buoyancy frequency
64			INTEGER I,J,Km1
65	jmc	1.3	_RL p5 , p125
66	mlosch	1.1	PARAMETER( p5=0.5, p125=0.125 )
67			_RL tempu, tempv
68			_RL epsilon
69			PARAMETER ( epsilon = 1.D-10 )
70	jmc	1.3	_RL RiFlux
71	mlosch	1.1	_RL buoyFreq
72			_RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73			_RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74			#ifdef PP81_SMOOTH_RI
75			_RL RiTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76			#endif /* PP81_SMOOTH_RI */
77	jmc	1.3	CEOP
78	mlosch	1.1
79			Km1 = MAX(1,K-1)
80			C Not sure what is correct for pressure coordinates:
81			C RiFlux is always correct because it quadratic
82			C buoyFreq should also be correct in pressure coordinates:
83			C N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
84			C So the following is wrong:
85			CML IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
86			CML kUp = MAX(1,K-1)
87			CML kDown = K
88			CML ELSEIF ( buoyancyRelation .EQ. 'OCEANIP') THEN
89			CML kUp = K
90			CML kDown = MAX(1,K-1)
91			CML ELSE
92			CML STOP 'PP81_RI_NUMBER: We should never reach this point'
93			CML ENDIF
94			C preparation: find density a kth and k-1st level
95	jmc	1.3	CALL FIND_RHO_2D(
96			I iMin, iMax, jMin, jMax, K,
97			I theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
98	mlosch	1.1	O rhoKm1,
99	jmc	1.3	I Km1, bi, bj, myThid )
100			CALL FIND_RHO_2D(
101			I iMin, iMax, jMin, jMax, K,
102			I theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
103	mlosch	1.1	O rhoK,
104	jmc	1.3	I K, bi, bj, myThid )
105	mlosch	1.1
106			C first step: calculate flux Richardson number.
107			C calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
108			DO J= jMin, jMax
109			DO I = iMin, iMax
110			tempu= .5*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
111			& - (uVel(I,J,K ,bi,bj)+uVel(I+1,J,K ,bi,bj)) )
112			& *recip_drC(K)
113			tempv= .5*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
114			& - (vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) )
115			& *recip_drC(K)
116			RiFlux = temputempu+tempvtempv
117
118			C
119			C calculate - g*(drho/dz)/rho0= N^2 .
120			C
121	jmc	1.2	buoyFreq = - gravitymass2rUnit
122	mlosch	1.1	& (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
123	jmc	1.3	C
124	mlosch	1.1	C calculates gradient Richardson number, bounded by
125			C a very large negative value.
126	jmc	1.3	C
127	mlosch	1.1	RiNumber(I,J) = buoyFreq/max(RiFlux,epsilon)
128			ENDDO
129			ENDDO
130	jmc	1.3
131	mlosch	1.1	#ifdef PP81_SMOOTH_RI
132			C average Richardson number horizontally
133			DO J=jMin,jMax
134			DO I=iMin,iMax
135			RiTmp(I,J) = RiNumber(I,J)
136			ENDDO
137	jmc	1.3	ENDDO
138	mlosch	1.1	DO J=1-Oly+1,sNy+Oly-1
139			DO I=1-Olx+1,sNx+Olx-1
140	jmc	1.3	RiNumber(I,J) = p5*RiTmp(I,J)
141			& + p125RiTmp(I-1,J) + p125RiTmp(I+1,J)
142			& + p125RiTmp(I,J-1) + p125RiTmp(I,J+1)
143	mlosch	1.1	ENDDO
144	jmc	1.3	ENDDO
145	mlosch	1.1	#endif /* PP81_SMOOTH_RI */
146
147			#endif /* ALLOW_PP81 */
148
149			RETURN
150			END
151