1 |
C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_u_bottomdrag.F,v 1.17 2015/01/04 16:16:32 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "MOM_COMMON_OPTIONS.h" |
5 |
#ifdef ALLOW_CTRL |
6 |
# include "CTRL_OPTIONS.h" |
7 |
#endif |
8 |
|
9 |
CBOP |
10 |
C !ROUTINE: MOM_U_BOTDRAG_IMPL |
11 |
|
12 |
C !INTERFACE: ========================================================== |
13 |
SUBROUTINE MOM_U_BOTDRAG_IMPL( |
14 |
I uFld, vFld, kappaRU, |
15 |
U cDrag, |
16 |
I bi, bj, myIter, myThid ) |
17 |
|
18 |
C !DESCRIPTION: |
19 |
C Add contribution from drag due to friction and the no-slip condition at bottom |
20 |
C to matrix main diagonal for implicit momentum solver |
21 |
C \begin{equation*} |
22 |
C C_{drag} = \frac{\Delta t}{\Delta r_f} (r_b + Cd |v| + \nu \frac{2}{\Delta r_c}) |
23 |
C \end{equation*} |
24 |
|
25 |
C !USES: =============================================================== |
26 |
IMPLICIT NONE |
27 |
#include "SIZE.h" |
28 |
#include "EEPARAMS.h" |
29 |
#include "PARAMS.h" |
30 |
#include "GRID.h" |
31 |
#ifdef ALLOW_CTRL |
32 |
# include "CTRL_FIELDS.h" |
33 |
#endif |
34 |
|
35 |
C !INPUT PARAMETERS: =================================================== |
36 |
C uFld :: zonal flow |
37 |
C vFld :: meridional flow |
38 |
C kappaRU :: vertical viscosity |
39 |
C bi,bj :: tile indices |
40 |
C myIter :: current iteration number |
41 |
C myThid :: thread number |
42 |
_RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
43 |
_RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
44 |
_RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) |
45 |
INTEGER bi, bj |
46 |
INTEGER myIter, myThid |
47 |
|
48 |
C !OUTPUT PARAMETERS: ================================================== |
49 |
C cDrag :: drag contribution to matrix main diagnonal |
50 |
_RL cDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
51 |
|
52 |
C !LOCAL VARIABLES: ==================================================== |
53 |
C i,j,k :: loop indices |
54 |
C KE :: Kinetic energy |
55 |
INTEGER i,j,k |
56 |
INTEGER kDown,kLowF,kBottom |
57 |
_RL viscFac, dragFac, uSq |
58 |
_RL recDrC |
59 |
_RL recDrF_bot(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
60 |
_RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
61 |
CEOP |
62 |
|
63 |
C- No-slip BCs impose a drag at bottom |
64 |
viscFac = 0. |
65 |
IF (no_slip_bottom) viscFac = 2. |
66 |
|
67 |
DO k=1,Nr |
68 |
|
69 |
IF ( usingZCoords ) THEN |
70 |
kBottom = Nr |
71 |
kDown = MIN(k+1,Nr) |
72 |
kLowF = k+1 |
73 |
c dragFac = mass2rUnit*rhoConst |
74 |
c dragFac = wUnit2rVel(k+1) |
75 |
dragFac = 1. _d 0 |
76 |
ELSE |
77 |
kBottom = 1 |
78 |
kDown = MAX(k-1,1) |
79 |
kLowF = k |
80 |
dragFac = mass2rUnit*rhoConst |
81 |
c dragFac = wUnit2rVel(k) |
82 |
ENDIF |
83 |
IF ( k.EQ.kBottom ) THEN |
84 |
recDrC = recip_drF(k) |
85 |
DO j=1-OLy,sNy+OLy |
86 |
DO i=1-OLx,sNx+OLx |
87 |
recDrF_bot(i,j) = _recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
88 |
ENDDO |
89 |
ENDDO |
90 |
ELSE |
91 |
recDrC = recip_drC(kLowF) |
92 |
DO j=1-OLy,sNy+OLy |
93 |
DO i=1-OLx,sNx+OLx |
94 |
recDrF_bot(i,j) = _recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
95 |
& * ( 1. _d 0 -_maskW(i,j,kDown,bi,bj) ) |
96 |
ENDDO |
97 |
ENDDO |
98 |
ENDIF |
99 |
|
100 |
C-- Linear bottom drag: |
101 |
DO j=1-OLy,sNy+OLy-1 |
102 |
DO i=1-OLx+1,sNx+OLx-1 |
103 |
cDrag(i,j,k) = cDrag(i,j,k) |
104 |
& + recDrF_bot(i,j) |
105 |
& *( bottomDragLinear*dragFac |
106 |
#ifdef ALLOW_BOTTOMDRAG_CONTROL |
107 |
& + halfRL*( bottomDragFld(i-1,j,bi,bj) |
108 |
& + bottomDragFld(i,j,bi,bj) )*dragFac |
109 |
#endif |
110 |
& )*deltaTMom |
111 |
ENDDO |
112 |
ENDDO |
113 |
|
114 |
C-- Add friction at the bottom (no-slip BC) |
115 |
IF ( no_slip_bottom .AND. bottomVisc_pCell ) THEN |
116 |
C- bottom friction accounts for true distance (including hFac) to the bottom |
117 |
DO j=1-OLy,sNy+OLy-1 |
118 |
DO i=1-OLx+1,sNx+OLx-1 |
119 |
cDrag(i,j,k) = cDrag(i,j,k) |
120 |
& + recDrF_bot(i,j) |
121 |
& *( kappaRU(i,j,kLowF)*recDrC*viscFac |
122 |
& *_recip_hFacW(i,j,k,bi,bj) |
123 |
& )*deltaTMom |
124 |
ENDDO |
125 |
ENDDO |
126 |
ELSEIF ( no_slip_bottom ) THEN |
127 |
C- ignores partial-cell reduction of the distance to the bottom |
128 |
DO j=1-OLy,sNy+OLy-1 |
129 |
DO i=1-OLx+1,sNx+OLx-1 |
130 |
cDrag(i,j,k) = cDrag(i,j,k) |
131 |
& + recDrF_bot(i,j) |
132 |
& *( kappaRU(i,j,kLowF)*recDrC*viscFac |
133 |
& )*deltaTMom |
134 |
ENDDO |
135 |
ENDDO |
136 |
ENDIF |
137 |
|
138 |
C-- Add quadratic bottom drag |
139 |
IF ( selectBotDragQuadr.EQ.0 ) THEN |
140 |
DO j=1-OLy,sNy+OLy-1 |
141 |
DO i=1-OLx,sNx+OLx-1 |
142 |
KE(i,j) = 0.25*( |
143 |
& ( uFld( i , j ,k)*uFld( i , j ,k)*_hFacW(i,j,k,bi,bj) |
144 |
& +uFld(i+1, j ,k)*uFld(i+1, j ,k)*_hFacW(i+1,j,k,bi,bj) ) |
145 |
& + ( vFld( i , j ,k)*vFld( i , j ,k)*_hFacS(i,j,k,bi,bj) |
146 |
& +vFld( i ,j+1,k)*vFld( i ,j+1,k)*_hFacS(i,j+1,k,bi,bj) ) |
147 |
& )*_recip_hFacC(i,j,k,bi,bj) |
148 |
ENDDO |
149 |
ENDDO |
150 |
C- average grid-cell-center KE to get velocity norm @ U.pt |
151 |
DO j=1-OLy,sNy+OLy-1 |
152 |
DO i=1-OLx+1,sNx+OLx-1 |
153 |
IF ( (KE(i,j)+KE(i-1,j)) .GT. 0. ) THEN |
154 |
cDrag(i,j,k) = cDrag(i,j,k) |
155 |
& + recDrF_bot(i,j) |
156 |
& *bottomDragQuadratic*SQRT(KE(i,j)+KE(i-1,j))*dragFac |
157 |
& *deltaTMom |
158 |
ENDIF |
159 |
ENDDO |
160 |
ENDDO |
161 |
ELSEIF ( selectBotDragQuadr.EQ.1 ) THEN |
162 |
C- calculate locally velocity norm @ U.pt (local U & 4 V averaged) |
163 |
DO j=1-OLy,sNy+OLy-1 |
164 |
DO i=1-OLx+1,sNx+OLx-1 |
165 |
uSq = uFld(i,j,k)*uFld(i,j,k) |
166 |
& + ( (vFld(i-1, j ,k)*vFld(i-1, j ,k)*hFacS(i-1, j ,k,bi,bj) |
167 |
& +vFld( i , j ,k)*vFld( i , j ,k)*hFacS( i , j ,k,bi,bj)) |
168 |
& + (vFld(i-1,j+1,k)*vFld(i-1,j+1,k)*hFacS(i-1,j+1,k,bi,bj) |
169 |
& +vFld( i ,j+1,k)*vFld( i ,j+1,k)*hFacS( i ,j+1,k,bi,bj)) |
170 |
& )*recip_hFacW(i,j,k,bi,bj)*0.25 _d 0 |
171 |
IF ( uSq.GT.zeroRL ) THEN |
172 |
cDrag(i,j,k) = cDrag(i,j,k) |
173 |
& + recDrF_bot(i,j) |
174 |
& *bottomDragQuadratic*SQRT(uSq)*dragFac |
175 |
& *deltaTMom |
176 |
ENDIF |
177 |
ENDDO |
178 |
ENDDO |
179 |
ELSEIF ( selectBotDragQuadr.EQ.2 ) THEN |
180 |
C- same as above but using wet-point method to average 4 V |
181 |
DO j=1-OLy,sNy+OLy-1 |
182 |
DO i=1-OLx+1,sNx+OLx-1 |
183 |
uSq = ( hFacS(i-1, j ,k,bi,bj) + hFacS( i , j ,k,bi,bj) ) |
184 |
& + ( hFacS(i-1,j+1,k,bi,bj) + hFacS( i ,j+1,k,bi,bj) ) |
185 |
IF ( uSq.GT.zeroRL ) THEN |
186 |
uSq = uFld(i,j,k)*uFld(i,j,k) |
187 |
& +( (vFld(i-1, j ,k)*vFld(i-1, j ,k)*hFacS(i-1, j ,k,bi,bj) |
188 |
& +vFld( i , j ,k)*vFld( i , j ,k)*hFacS( i , j ,k,bi,bj)) |
189 |
& + (vFld(i-1,j+1,k)*vFld(i-1,j+1,k)*hFacS(i-1,j+1,k,bi,bj) |
190 |
& +vFld( i ,j+1,k)*vFld( i ,j+1,k)*hFacS( i ,j+1,k,bi,bj)) |
191 |
& )/uSq |
192 |
ELSE |
193 |
uSq = uFld(i,j,k)*uFld(i,j,k) |
194 |
ENDIF |
195 |
IF ( uSq.GT.zeroRL ) THEN |
196 |
cDrag(i,j,k) = cDrag(i,j,k) |
197 |
& + recDrF_bot(i,j) |
198 |
& *bottomDragQuadratic*SQRT(uSq)*dragFac |
199 |
& *deltaTMom |
200 |
ENDIF |
201 |
ENDDO |
202 |
ENDDO |
203 |
ELSEIF ( selectBotDragQuadr.NE.-1 ) THEN |
204 |
STOP 'MOM_U_BOTDRAG_IMPL: invalid selectBotDragQuadr value' |
205 |
ENDIF |
206 |
|
207 |
c#ifdef ALLOW_DIAGNOSTICS |
208 |
c IF (useDiagnostics) THEN |
209 |
c CALL DIAGNOSTICS_FILL(uDragTerms,'UBotDrag',k,1,2,bi,bj,myThid) |
210 |
c ENDIF |
211 |
c#endif /* ALLOW_DIAGNOSTICS */ |
212 |
|
213 |
C- end k loop |
214 |
ENDDO |
215 |
|
216 |
RETURN |
217 |
END |