-
Notifications
You must be signed in to change notification settings - Fork 237
/
my82_ri_number.F
142 lines (127 loc) · 4.56 KB
/
my82_ri_number.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#include "MY82_OPTIONS.h"
CBOP
C !ROUTINE: MY82_RI_NUMBER
C !INTERFACE: ==========================================================
subroutine MY82_RI_NUMBER(
I bi, bj, K, iMin, iMax, jMin, jMax,
O RiNumber, buoyFreq, vertShear,
I myTime, myThid )
C !DESCRIPTION: \bv
C /==========================================================\
C | SUBROUTINE MY82_RI_NUMBER |
C | o Compute gradient Richardson number for Mellor and |
C | Yamada (1981) turbulence model |
C \==========================================================/
IMPLICIT NONE
c
c-------------------------------------------------------------------------
c \ev
C !USES: ===============================================================
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "MY82.h"
#include "GRID.h"
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
C buoyFreq - (output) (neg.) buoyancy frequency -N^2
C vertShear - (output) vertical shear of velocity
INTEGER bi, bj, iMin, iMax, jMin, jMax, k
INTEGER myThid
_RL myTime
_RL RiNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL buoyFreq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vertShear (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef ALLOW_MY82
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
INTEGER I,J,Km1
_RL p5 , p125
PARAMETER( p5=0.5D0, p125=0.125D0 )
_RL tempu, tempv
_RL epsilon
PARAMETER ( epsilon = 1.D-10 )
_RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef MY82_SMOOTH_RI
_RL RiTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif /* MY82_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 'MY82_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 _d 0*( 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 _d 0*( 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)
vertShear(I,J) = tempu*tempu+tempv*tempv
C
C calculate g*(drho/dz)/rho0= -N^2 .
C
buoyFreq(I,J) = 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(I,J)/max(vertShear(I,J),epsilon)
ENDDO
ENDDO
#ifdef MY82_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 /* MY82_SMOOTH_RI */
#endif /* ALLOW_MY82 */
RETURN
END