/[MITgcm]/MITgcm/pkg/my82/my82_ri_number.F
ViewVC logotype

Annotation of /MITgcm/pkg/my82/my82_ri_number.F

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

Revision 1.5 - (hide annotations) (download)
Tue Jan 20 00:26:04 2009 UTC (15 years, 4 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, checkpoint62, 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, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, HEAD
Changes since 1.4: +3 -10 lines 
add missing "_d 0"
1 jmc 1.5 C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_ri_number.F,v 1.4 2008/08/11 22:28:06 jmc Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "MY82_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: MY82_RI_NUMBER
8    
9     C !INTERFACE: ==========================================================
10     subroutine MY82_RI_NUMBER(
11 jmc 1.4 I bi, bj, K, iMin, iMax, jMin, jMax,
12 mlosch 1.1 O RiNumber, buoyFreq, vertShear,
13     I myTime, myThid )
14    
15     C !DESCRIPTION: \bv
16     C /==========================================================\
17     C | SUBROUTINE MY82_RI_NUMBER |
18     C | o Compute gradient Richardson number for Mellor and |
19     C | Yamada (1981) turbulence model |
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 "MY82.h"
33     #include "GRID.h"
34    
35     C !INPUT PARAMETERS: ===================================================
36     C Routine arguments
37     C bi, bj - array indices on which to apply calculations
38 jmc 1.4 C iMin, iMax, jMin, jMax
39 mlosch 1.1 C - array boundaries
40     C k - depth level
41     C myTime - Current time in simulation
42     C RiNumber - (output) Richardson number
43     C buoyFreq - (output) (neg.) buoyancy frequency -N^2
44     C vertShear - (output) vertical shear of velocity
45     INTEGER bi, bj, iMin, iMax, jMin, jMax, k
46     INTEGER myThid
47     _RL myTime
48     _RL RiNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49     _RL buoyFreq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50     _RL vertShear (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51    
52     #ifdef ALLOW_MY82
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     INTEGER I,J,Km1
61 jmc 1.4 _RL p5 , p125
62 mlosch 1.2 PARAMETER( p5=0.5D0, p125=0.125D0 )
63 mlosch 1.1 _RL tempu, tempv
64     _RL epsilon
65     PARAMETER ( epsilon = 1.D-10 )
66     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     #ifdef MY82_SMOOTH_RI
69     _RL RiTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     #endif /* MY82_SMOOTH_RI */
71 jmc 1.4 CEOP
72 mlosch 1.1
73     Km1 = MAX(1,K-1)
74     C Not sure what is correct for pressure coordinates:
75     C RiFlux is always correct because it quadratic
76     C buoyFreq should also be correct in pressure coordinates:
77     C N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
78     C So the following is wrong:
79     CML IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
80     CML kUp = MAX(1,K-1)
81     CML kDown = K
82     CML ELSEIF ( buoyancyRelation .EQ. 'OCEANIP') THEN
83     CML kUp = K
84     CML kDown = MAX(1,K-1)
85     CML ELSE
86     CML STOP 'MY82_RI_NUMBER: We should never reach this point'
87     CML ENDIF
88     C preparation: find density a kth and k-1st level
89 jmc 1.4 CALL FIND_RHO_2D(
90     I iMin, iMax, jMin, jMax, K,
91     I theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
92 mlosch 1.1 O rhoKm1,
93 jmc 1.4 I Km1, bi, bj, myThid )
94     CALL FIND_RHO_2D(
95     I iMin, iMax, jMin, jMax, K,
96     I theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
97 mlosch 1.1 O rhoK,
98 jmc 1.4 I K, bi, bj, myThid )
99 mlosch 1.1
100     C first step: calculate flux Richardson number.
101     C calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
102     DO J= jMin, jMax
103     DO I = iMin, iMax
104 jmc 1.5 tempu= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
105 mlosch 1.1 & - (uVel(I,J,K ,bi,bj)+uVel(I+1,J,K ,bi,bj)) )
106     & *recip_drC(K)
107 jmc 1.5 tempv= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
108 mlosch 1.1 & - (vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) )
109     & *recip_drC(K)
110     vertShear(I,J) = tempu*tempu+tempv*tempv
111    
112     C
113     C calculate g*(drho/dz)/rho0= -N^2 .
114     C
115 jmc 1.3 buoyFreq(I,J) = gravity*mass2rUnit *
116 mlosch 1.1 & (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
117 jmc 1.4 C
118 mlosch 1.1 C calculates gradient Richardson number, bounded by
119     C a very large negative value.
120 jmc 1.4 C
121 mlosch 1.1 RiNumber(I,J) = -buoyFreq(I,J)/max(vertShear(I,J),epsilon)
122 jmc 1.4
123 mlosch 1.1 ENDDO
124     ENDDO
125 jmc 1.4
126 mlosch 1.1 #ifdef MY82_SMOOTH_RI
127     C average Richardson number horizontally
128     DO J=jMin,jMax
129     DO I=iMin,iMax
130     RiTmp(I,J) = RiNumber(I,J)
131     ENDDO
132 jmc 1.4 ENDDO
133 mlosch 1.1 DO J=1-Oly+1,sNy+Oly-1
134     DO I=1-Olx+1,sNx+Olx-1
135 jmc 1.4 RiNumber(I,J) = p5*RiTmp(I,J)
136     & + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
137     & + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
138 mlosch 1.1 ENDDO
139 jmc 1.4 ENDDO
140 mlosch 1.1 #endif /* MY82_SMOOTH_RI */
141    
142     #endif /* ALLOW_MY82 */
143    
144     RETURN
145     END
  ViewVC Help
Powered by ViewVC 1.1.22