/[MITgcm]/MITgcm/pkg/pp81/pp81_ri_number.F
ViewVC logotype

Annotation of /MITgcm/pkg/pp81/pp81_ri_number.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Sun Jan 3 19:17:01 2010 UTC (14 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, HEAD
Changes since 1.3: +1 -3 lines
avoid unused variables

1 jmc 1.4 C $Header: /u/gcmpack/MITgcm/pkg/pp81/pp81_ri_number.F,v 1.3 2008/08/11 22:28:06 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     #endif /* ALLOW_AUTODIFF_TAMC */
38    
39     C !INPUT PARAMETERS: ===================================================
40     C Routine arguments
41     C bi, bj - array indices on which to apply calculations
42 jmc 1.3 C iMin, iMax, jMin, jMax
43 mlosch 1.1 C - array boundaries
44     C k - depth level
45     C myTime - Current time in simulation
46     C RiNumber - (output) Richardson number
47     INTEGER bi, bj, iMin, iMax, jMin, jMax, k
48     INTEGER myThid
49     _RL myTime
50     _RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51    
52     #ifdef ALLOW_PP81
53    
54     C !LOCAL VARIABLES: ====================================================
55     C I,J,kUp - loop indices
56     C p0-125 - averaging coefficients
57     C tempu, tempv - temporary variables
58     C rhoK, rhoKm1 - Density below and above current interface
59     C epsilon - small number
60 jmc 1.3 C RiFlux - denominator of Richardson number
61 mlosch 1.1 C BuoyFreq - buoyancy frequency
62     INTEGER I,J,Km1
63 jmc 1.3 _RL p5 , p125
64 mlosch 1.1 PARAMETER( p5=0.5, p125=0.125 )
65     _RL tempu, tempv
66     _RL epsilon
67     PARAMETER ( epsilon = 1.D-10 )
68 jmc 1.3 _RL RiFlux
69 mlosch 1.1 _RL buoyFreq
70     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     #ifdef PP81_SMOOTH_RI
73     _RL RiTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74     #endif /* PP81_SMOOTH_RI */
75 jmc 1.3 CEOP
76 mlosch 1.1
77     Km1 = MAX(1,K-1)
78     C Not sure what is correct for pressure coordinates:
79     C RiFlux is always correct because it quadratic
80     C buoyFreq should also be correct in pressure coordinates:
81     C N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
82     C So the following is wrong:
83     CML IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
84     CML kUp = MAX(1,K-1)
85     CML kDown = K
86     CML ELSEIF ( buoyancyRelation .EQ. 'OCEANIP') THEN
87     CML kUp = K
88     CML kDown = MAX(1,K-1)
89     CML ELSE
90     CML STOP 'PP81_RI_NUMBER: We should never reach this point'
91     CML ENDIF
92     C preparation: find density a kth and k-1st level
93 jmc 1.3 CALL FIND_RHO_2D(
94     I iMin, iMax, jMin, jMax, K,
95     I theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
96 mlosch 1.1 O rhoKm1,
97 jmc 1.3 I Km1, bi, bj, myThid )
98     CALL FIND_RHO_2D(
99     I iMin, iMax, jMin, jMax, K,
100     I theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
101 mlosch 1.1 O rhoK,
102 jmc 1.3 I K, bi, bj, myThid )
103 mlosch 1.1
104     C first step: calculate flux Richardson number.
105     C calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
106     DO J= jMin, jMax
107     DO I = iMin, iMax
108     tempu= .5*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
109     & - (uVel(I,J,K ,bi,bj)+uVel(I+1,J,K ,bi,bj)) )
110     & *recip_drC(K)
111     tempv= .5*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
112     & - (vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) )
113     & *recip_drC(K)
114     RiFlux = tempu*tempu+tempv*tempv
115    
116     C
117     C calculate - g*(drho/dz)/rho0= N^2 .
118     C
119 jmc 1.2 buoyFreq = - gravity*mass2rUnit *
120 mlosch 1.1 & (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
121 jmc 1.3 C
122 mlosch 1.1 C calculates gradient Richardson number, bounded by
123     C a very large negative value.
124 jmc 1.3 C
125 mlosch 1.1 RiNumber(I,J) = buoyFreq/max(RiFlux,epsilon)
126     ENDDO
127     ENDDO
128 jmc 1.3
129 mlosch 1.1 #ifdef PP81_SMOOTH_RI
130     C average Richardson number horizontally
131     DO J=jMin,jMax
132     DO I=iMin,iMax
133     RiTmp(I,J) = RiNumber(I,J)
134     ENDDO
135 jmc 1.3 ENDDO
136 mlosch 1.1 DO J=1-Oly+1,sNy+Oly-1
137     DO I=1-Olx+1,sNx+Olx-1
138 jmc 1.3 RiNumber(I,J) = p5*RiTmp(I,J)
139     & + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
140     & + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
141 mlosch 1.1 ENDDO
142 jmc 1.3 ENDDO
143 mlosch 1.1 #endif /* PP81_SMOOTH_RI */
144    
145     #endif /* ALLOW_PP81 */
146    
147     RETURN
148     END
149    

  ViewVC Help
Powered by ViewVC 1.1.22