-
Notifications
You must be signed in to change notification settings - Fork 237
/
mom_uv_boundary.F
160 lines (147 loc) · 5.37 KB
/
mom_uv_boundary.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
#include "MOM_FLUXFORM_OPTIONS.h"
CBOP
C !ROUTINE: MOM_UV_BOUNDARY
C !INTERFACE: ==========================================================
SUBROUTINE MOM_UV_BOUNDARY (
I bi,bj,k,
I uFld, vFld,
O uBnd, vBnd,
I myTime, myIter, myThid )
C !DESCRIPTION:
C Set velocity at a boundary for a momentum conserving advection
C Note: really conserve momentum when "steps" (vertical plane)
C or coastline (horizontal plane) are only 1 grid-point wide.
C !USES: ===============================================================
C == Global variables ==
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
C !INPUT PARAMETERS: ===================================================
C bi,bj :: tile indices
C k :: vertical level
C uFld :: zonal velocity
C vFld :: meridional velocity
C myTime :: current time
C myIter :: current iteration number
C myThid :: My Thread Id. number
INTEGER bi,bj
INTEGER k
_RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL myTime
INTEGER myIter
INTEGER myThid
C !OUTPUT PARAMETERS: ==================================================
C uBnd :: zonal velocity extended to boundaries
C vBnd :: meridional velocity extended to boundaries
_RL uBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef MOM_BOUNDARY_CONSERVE
C !LOCAL VARIABLES: ====================================================
C i,j :: loop indices
INTEGER i,j
INTEGER km1,kp1
_RL maskM1, maskP1
_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL aBndU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL aBndV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL aBndW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tmpVar
LOGICAL useMomBndConserve
PARAMETER ( useMomBndConserve = .TRUE. )
CEOP
C Initialise output array
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
uBnd(i,j) = uFld(i,j,k,bi,bj)
vBnd(i,j) = vFld(i,j,k,bi,bj)
ENDDO
ENDDO
IF ( useMomBndConserve ) THEN
C- Initialise intermediate arrays:
km1 = MAX( k-1, 1 )
kp1 = MIN( k+1, Nr )
maskM1 = 1.
maskP1 = 1.
IF ( k.EQ.1 ) maskM1 = 0.
IF ( k.EQ.Nr ) maskP1 = 0.
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
aBndU(i,j) = 0.
aBndV(i,j) = 0.
aBndW(i,j) = 0.
ENDDO
ENDDO
C- Calculate Divergence in 3 directions:
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
uTrans(i,j) = uFld(i,j,k,bi,bj)
& * dyG(i,j,bi,bj)*deepFacC(k)
& * drF(k)*hFacW(i,j,k,bi,bj)*rhoFacC(k)
vTrans(i,j) = vFld(i,j,k,bi,bj)
& * dxG(i,j,bi,bj)*deepFacC(k)
& * drF(k)*hFacS(i,j,k,bi,bj)*rhoFacC(k)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
aBndU(i,j) = uTrans(i+1,j)-uTrans(i,j)
aBndV(i,j) = vTrans(i,j+1)-vTrans(i,j)
aBndW(i,j) = ABS(aBndU(i,j)+aBndV(i,j))
aBndU(i,j) = ABS(aBndU(i,j))
aBndV(i,j) = ABS(aBndV(i,j))
ENDDO
ENDDO
C- Normalise by the sum:
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
tmpVar = aBndU(i,j)+aBndV(i,j)+aBndW(i,j)
IF ( tmpVar.GT.0. ) THEN
tmpVar = 1. _d 0 / tmpVar
aBndU(i,j) = aBndU(i,j)*tmpVar
aBndV(i,j) = aBndV(i,j)*tmpVar
aBndW(i,j) = aBndW(i,j)*tmpVar
ENDIF
ENDDO
ENDDO
C- At a boundary, replace uFld,vFld by a weighted average
C Note: multiply by 2 to cancel the 1/2 factor in advections S/R
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF (maskW(i,j,k,bi,bj).EQ.0.) THEN
C Note: only 1 set of aBnd_U,V,W is non-zero (either i-1 or i)
C and only 1 uFld is non-zero (either i-1 or i+1)
C and only 1 uFld is non-zero (either k-1 or k+1)
uBnd(i,j) = (
& (aBndU(i-1,j)+aBndU(i,j))
& *(uFld(i-1,j,k,bi,bj)+uFld(i+1,j,k,bi,bj))
& +(aBndV(i-1,j)+aBndV(i,j))
& *(uFld(i,j-1,k,bi,bj)+uFld(i,j+1,k,bi,bj))
& +(aBndW(i-1,j)+aBndW(i,j))
& *(uFld(i,j,km1,bi,bj)*maskM1
& +uFld(i,j,kp1,bi,bj)*maskP1)
& )*2. _d 0
ENDIF
IF (maskS(i,j,k,bi,bj).EQ.0.) THEN
C Note: only 1 set of aBnd_U,V,W is non-zero (either j-1 or j)
C and only 1 vFld is non-zero (either j-1 or j+1)
C and only 1 vFld is non-zero (either k-1 or k+1)
vBnd(i,j) = (
& (aBndU(i,j-1)+aBndU(i,j))
& *(vFld(i-1,j,k,bi,bj)+vFld(i+1,j,k,bi,bj))
& +(aBndV(i,j-1)+aBndV(i,j))
& *(vFld(i,j-1,k,bi,bj)+vFld(i,j+1,k,bi,bj))
& +(aBndW(i,j-1)+aBndW(i,j))
& *(vFld(i,j,km1,bi,bj)*maskM1
& +vFld(i,j,kp1,bi,bj)*maskP1)
& )*2. _d 0
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* MOM_BOUNDARY_CONSERVE */
RETURN
END