forked from darwinproject/darwin3
-
Notifications
You must be signed in to change notification settings - Fork 1
/
mom_vi_coriolis.F
193 lines (182 loc) · 7.49 KB
/
mom_vi_coriolis.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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
#include "MOM_VECINV_OPTIONS.h"
CBOP
C !ROUTINE: MOM_VI_CORIOLIS
C !INTERFACE: ==========================================================
SUBROUTINE MOM_VI_CORIOLIS(
I bi, bj, k,
I uFld, vFld, hFacZ, r_hFacZ,
O uCoriolisTerm, vCoriolisTerm,
I myThid )
C !DESCRIPTION:
C Calculates the 2 horizontal components of Coriolis acceleration
C !USES: ===============================================================
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"
#include "PARAMS.h"
C !INPUT PARAMETERS: ===================================================
C bi, bj :: current tile indices
C k :: current vertical level
C uFld, vFld :: local copy of horizontal velocity (u & v components)
C hFacZ :: hFac thickness factor at corner location
C r_hFacZ :: reciprocal hFac thickness factor at corner location
C myThid :: my Thread Id number
INTEGER bi,bj,k
_RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER myThid
C !OUTPUT PARAMETERS: ==================================================
C uCoriolisTerm :: Coriolis tendency for u-component momentum
C vCoriolisTerm :: Coriolis tendency for v-component momentum
_RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C !LOCAL VARIABLES: ====================================================
C i, j :: loop indices
C uBarXY, vBarXY :: averaged (in X & Y direction) of u & v velocity
C uBarYm, uBarYp :: u velocity, Y direction averaged (at i and i+1)
C vBarXm, vBarXp :: v velocity, X direction averaged (at j and j+1)
C epsil :: small number
INTEGER i,j
_RL uBarXY, vBarXY
_RL uBarYm, uBarYp, vBarXm, vBarXp
_RS epsil
CEOP
epsil = 1. _d -9
IF ( selectCoriScheme .EQ. 0 ) THEN
C- Simple average, no hFac :
DO j=1-OLy,sNy+OLy-1
DO i=2-OLx,sNx+OLx
vBarXY=0.25*(
& (vFld( i , j )*dxG( i , j ,bi,bj)
& +vFld(i-1, j )*dxG(i-1, j ,bi,bj))
& +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)
& +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj))
& )
uCoriolisTerm(i,j)=
& +0.5*( fCoriG(i,j,bi,bj)+fCoriG(i,j+1,bi,bj)
& )*vBarXY*recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 1 ) THEN
C- Partial-cell generalization of the Wet-point average method :
C (formerly called: useJamartWetPoints)
DO j=1-OLy,sNy+OLy-1
DO i=2-OLx,sNx+OLx
vBarXY=(
& (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
& +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
& +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
& +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj)))
& / MAX( epsil,(_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
& +(_hFacS(i,j+1,k,bi,bj)+_hFacS(i-1,j+1,k,bi,bj)) )
uCoriolisTerm(i,j)=
& +0.5*( fCoriG(i,j,bi,bj)+fCoriG(i,j+1,bi,bj)
& )*vBarXY*recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 2 ) THEN
C- hFac weighted average :
DO j=1-OLy,sNy+OLy-1
DO i=2-OLx,sNx+OLx
vBarXY=0.25*(
& (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
& +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
& +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
& +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj))
& )
uCoriolisTerm(i,j)=
& +0.5*( fCoriG(i,j,bi,bj)+fCoriG(i,j+1,bi,bj)
& )*vBarXY*recip_dxC(i,j,bi,bj)*_recip_hFacW(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 3 ) THEN
C- Energy-conserving discretisation with hFac weighted average :
DO j=1-OLy,sNy+OLy-1
DO i=2-OLx,sNx+OLx
vBarXm = halfRL *(
& vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
& +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
vBarXp = halfRL *(
& vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
& +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
uCoriolisTerm(i,j) = +0.5 _d 0
& *( vBarXm*fCoriG(i, j ,bi,bj)
& +vBarXp*fCoriG(i,j+1,bi,bj)
& )*recip_dxC(i,j,bi,bj)*_recip_hFacW(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
STOP 'MOM_VI_CORIOLIS: invalid selectCoriScheme'
ENDIF
IF ( selectCoriScheme .EQ. 0 ) THEN
C- Simple average, no hFac :
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
uBarXY=0.25*(
& (uFld( i , j )*dyG( i , j ,bi,bj)
& +uFld( i ,j-1)*dyG( i ,j-1,bi,bj))
& +(uFld(i+1, j )*dyG(i+1, j ,bi,bj)
& +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj))
& )
vCoriolisTerm(i,j)=
& -0.5*( fCoriG(i,j,bi,bj)+fCoriG(i+1,j,bi,bj)
& )*uBarXY*recip_dyC(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 1 ) THEN
C- Partial-cell generalization of the Wet-point average method :
C (formerly called: useJamartWetPoints)
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
uBarXY=(
& (uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj)
& +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj))
& +(uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj)
& +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj)))
& / MAX( epsil,(_hFacW( i ,j,k,bi,bj)+_hFacW( i ,j-1,k,bi,bj))
& +(_hFacW(i+1,j,k,bi,bj)+_hFacW(i+1,j-1,k,bi,bj)) )
vCoriolisTerm(i,j)=
& -0.5*( fCoriG(i,j,bi,bj)+fCoriG(i+1,j,bi,bj)
& )*uBarXY*recip_dyC(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 2 ) THEN
C- hFac weighted average :
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
uBarXY=0.25*(
& (uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj)
& +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj))
& +(uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj)
& +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj))
& )
vCoriolisTerm(i,j)=
& -0.5*( fCoriG(i,j,bi,bj)+fCoriG(i+1,j,bi,bj)
& )*uBarXY*recip_dyC(i,j,bi,bj)*_recip_hFacS(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectCoriScheme .EQ. 3 ) THEN
C- Energy-conserving discretisation with hFac weighted average :
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
uBarYm = halfRL *(
& uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj)
& +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) )
uBarYp = halfRL *(
& uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj)
& +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) )
vCoriolisTerm(i,j) = -0.5 _d 0
& *( uBarYm*fCoriG( i ,j,bi,bj)
& +uBarYp*fCoriG(i+1,j,bi,bj)
& )*recip_dyC(i,j,bi,bj)*_recip_hFacS(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
STOP 'MOM_VI_CORIOLIS: invalid selectCoriScheme'
ENDIF
RETURN
END