/[MITgcm]/MITgcm/pkg/shap_filt/shap_filt_uv_s2.F
ViewVC logotype

Contents of /MITgcm/pkg/shap_filt/shap_filt_uv_s2.F

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

Revision 1.6 - (show annotations) (download)
Mon Mar 4 01:32:55 2002 UTC (22 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47e_post, checkpoint46l_post, checkpoint46g_pre, checkpoint47c_post, checkpoint50c_post, checkpoint46f_post, checkpoint48e_post, checkpoint50c_pre, checkpoint44f_post, checkpoint46b_post, checkpoint48i_post, checkpoint46l_pre, checkpoint51, checkpoint50, checkpoint50d_post, checkpoint50b_pre, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47a_post, checkpoint48d_pre, checkpoint47i_post, checkpoint47d_post, checkpoint46d_pre, checkpoint48d_post, checkpoint48f_post, checkpoint45d_post, checkpoint46j_pre, checkpoint44h_pre, checkpoint48h_post, checkpoint51b_pre, checkpoint46a_post, checkpoint47g_post, checkpoint46j_post, checkpoint46k_post, checkpoint48a_post, checkpoint45a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint47j_post, branch-exfmods-tag, checkpoint44g_post, checkpoint46e_pre, checkpoint48c_post, checkpoint45b_post, checkpoint46b_pre, checkpoint51b_post, checkpoint51c_post, checkpoint46c_pre, checkpoint46, checkpoint47b_post, checkpoint46h_pre, checkpoint46m_post, checkpoint46a_pre, checkpoint50g_post, checkpoint45c_post, checkpoint44h_post, checkpoint46g_post, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, checkpoint47f_post, checkpoint50e_post, checkpoint46i_post, checkpoint46c_post, checkpoint50d_pre, checkpoint46e_post, checkpoint47, checkpoint45, checkpoint48, checkpoint49, checkpoint46h_post, checkpoint48g_post, checkpoint47h_post, checkpoint44f_pre, checkpoint46d_post, checkpoint50b_post, checkpoint51a_post
Branch point for: branch-exfmods-curt
Changes since 1.5: +51 -19 lines 
o parameter Shap_noSlip replace CPP option NO_SLIP_SHAP
o add filter time scale (missing in some Shap_filt S/R)
o add exchange at the end (missing in some Shapi_filt S/R)
o working arrays (in common block SHAP_FILT.h) become argument
  of Shap_filt S/R.  => Enable to filter 2D fields.
1 C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_uv_s2.F,v 1.5 2002/02/08 21:37:05 jmc Exp $
2 C $Name: $
3
4 #include "SHAP_FILT_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: SHAP_FILT_UV_S2
8 C !INTERFACE:
9 SUBROUTINE SHAP_FILT_UV_S2(
10 U uFld, vFld, tmpFldU, tmpFldV,
11 I kSize, myTime, myThid )
12 C !DESCRIPTION: \bv
13 C *==========================================================*
14 C | S/R SHAP_FILT_UV_S2
15 C | o Applies Shapiro filter to velocity field (u & v).
16 C | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
17 C | o Options for computational filter (no grid spacing)
18 C | or physical space filter (with grid spacing) or both.
19 C *==========================================================*
20 C \ev
21
22 C !USES:
23 IMPLICIT NONE
24
25 C == Global variables ===
26 #include "SIZE.h"
27 #include "EEPARAMS.h"
28 #include "PARAMS.h"
29 #include "GRID.h"
30 #include "SHAP_FILT.h"
31
32 C !INPUT/OUTPUT PARAMETERS:
33 C == Routine arguments
34 C uFld :: velocity field (U component) on which filter applies
35 C vFld :: velocity field (V component) on which filter applies
36 C tmpFldU :: working temporary array
37 C tmpFldV :: working temporary array
38 C kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
39 C myTime :: Current time in simulation
40 C myThid :: Thread number for this instance of SHAP_FILT_UV_S2
41 INTEGER kSize
42 _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
43 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
44 _RL tmpFldU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
45 _RL tmpFldV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
46 _RL myTime
47 INTEGER myThid
48
49 #ifdef ALLOW_SHAP_FILT
50
51 C------
52 C Combine computational Filter of Div & Vorticity
53 C and Physical Filter of U,V field
54 C nShapUVPhys = 0 ==> use only computational Filter
55 C nShapUVPhys = 1 ==> compute Div & Vort. with Grid factors,
56 C Filter Div & Vort. Numerically (power nShapUV-1)
57 C and return filtered U.V in physical space
58 C nShapUVPhys = nShapUV ==> Filter in Physical space only (power nShapUV)
59 C------
60
61 C !LOCAL VARIABLES:
62 C == Local variables ==
63 INTEGER bi,bj,k,i,j,N
64 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RS maskZ
69 CEOP
70
71 IF (nShapUV.GT.0 .AND. Shap_uvtau.GT.0.) THEN
72
73 DO bj=myByLo(myThid),myByHi(myThid)
74 DO bi=myBxLo(myThid),myBxHi(myThid)
75 DO K=1,kSize
76 DO J=1-Oly,sNy+Oly
77 DO I=1-Olx,sNx+Olx
78 tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
79 & *_maskW(i,j,k,bi,bj)
80 tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
81 & *_maskS(i,j,k,bi,bj)
82 ENDDO
83 ENDDO
84 ENDDO
85 ENDDO
86 ENDDO
87
88 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89
90 C [d_xx+d_yy]^n tmpFld
91
92 DO N=1,nShapUV
93
94 IF (kSize.EQ.Nr) THEN
95 CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
96 ELSE
97 CALL EXCH_UV_XY_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
98 ENDIF
99
100 DO bj=myByLo(myThid),myByHi(myThid)
101 DO bi=myBxLo(myThid),myBxHi(myThid)
102 DO K=1,kSize
103
104 C [d_xx+d_yy] tmpFld
105 IF (N.LE.nShapUVPhys) THEN
106 CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
107 CALL MOM_VI_CALC_HDIV(bi,bj,k,
108 I tmpFldU(1-OLx,1-OLy,k,bi,bj),
109 I tmpFldV(1-OLx,1-OLy,k,bi,bj),
110 & hDiv,myThid)
111 CALL MOM_VI_CALC_RELVORT3(bi,bj,k,
112 I tmpFldU(1-OLx,1-OLy,k,bi,bj),
113 I tmpFldV(1-OLx,1-OLy,k,bi,bj),
114 & hFacZ,vort3,myThid)
115 ELSE
116 C- replace Physical calc Div & Vort by computational one :
117 DO J=0,sNy+1
118 DO I=0,sNx+1
119 hDiv(i,j)=tmpFldU(i+1,j,k,bi,bj)-tmpFldU(i,j,k,bi,bj)
120 & +tmpFldV(i,j+1,k,bi,bj)-tmpFldV(i,j,k,bi,bj)
121 ENDDO
122 ENDDO
123 DO J=1,sNy+1
124 DO I=1,sNx+1
125 vort3(i,j)=(tmpFldV(i,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj)
126 & -tmpFldU(i,j,k,bi,bj)+tmpFldU(i,j-1,k,bi,bj)
127 & )
128 maskZ = (maskW(i,j,k,bi,bj)+maskW(i,j-1,k,bi,bj))
129 & *(maskS(i,j,k,bi,bj)+maskS(i-1,j,k,bi,bj))
130 IF (maskZ.LT.1.) vort3(i,j)=0.
131 ENDDO
132 ENDDO
133
134 C Special stuff for Cubed Sphere
135 IF (useCubedSphereExchange) THEN
136 c---
137 I=1
138 J=1
139 maskZ = maskW(i,j,k,bi,bj)+maskW(i,j-1,k,bi,bj)
140 & +maskS(i,j,k,bi,bj)
141 IF (maskZ.GE.2.) THEN
142 vort3(I,J)=
143 & tmpFldV(I,J,k,bi,bj)
144 & -tmpFldU(I,J,k,bi,bj)
145 & +tmpFldU(I,J-1,k,bi,bj)
146 vort3(I,J)=vort3(I,J)*4.d0/3.d0
147 ELSE
148 vort3(I,J)=0.
149 ENDIF
150 c---
151 I=sNx+1
152 J=1
153 maskZ = maskW(i,j,k,bi,bj)+maskW(i,j-1,k,bi,bj)
154 & +maskS(i-1,j,k,bi,bj)
155 IF (maskZ.GE.2.) THEN
156 vort3(I,J)=
157 & -tmpFldV(I-1,J,k,bi,bj)
158 & -tmpFldU(I,J,k,bi,bj)
159 & +tmpFldU(I,J-1,k,bi,bj)
160 vort3(I,J)=vort3(I,J)*4.d0/3.d0
161 ELSE
162 vort3(I,J)=0.
163 ENDIF
164 c---
165 I=1
166 J=sNy+1
167 maskZ = maskW(i,j,k,bi,bj)+maskW(i,j-1,k,bi,bj)
168 & +maskS(i,j,k,bi,bj)
169 IF (maskZ.GE.2.) THEN
170 vort3(I,J)=
171 & tmpFldV(I,J,k,bi,bj)
172 & -tmpFldU(I,J,k,bi,bj)
173 & +tmpFldU(I,J-1,k,bi,bj)
174 vort3(I,J)=vort3(I,J)*4.d0/3.d0
175 ELSE
176 vort3(I,J)=0.
177 ENDIF
178 c---
179 I=sNx+1
180 J=sNy+1
181 maskZ = maskW(i,j,k,bi,bj)+maskW(i,j-1,k,bi,bj)
182 & +maskS(i-1,j,k,bi,bj)
183 IF (maskZ.GE.2.) THEN
184 vort3(I,J)=
185 & -tmpFldV(I-1,J,k,bi,bj)
186 & -tmpFldU(I,J,k,bi,bj)
187 & +tmpFldU(I,J-1,k,bi,bj)
188 vort3(I,J)=vort3(I,J)*4.d0/3.d0
189 ELSE
190 vort3(I,J)=0.
191 ENDIF
192 c---
193 ENDIF
194 ENDIF
195
196 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
197
198 IF (N.GT.nShapUV-nShapUVPhys) THEN
199 CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
200 CALL MOM_VI_DEL2UV(
201 I bi,bj,k,hDiv,vort3,hFacZ,
202 O tmpFldU(1-OLx,1-OLy,k,bi,bj),
203 O tmpFldV(1-OLx,1-OLy,k,bi,bj),
204 I myThid)
205 IF (Shap_uvLength.EQ.0.) THEN
206 DO J=1,sNy+1
207 DO I=1,sNx+1
208 tmpFldU(i,j,k,bi,bj) = -0.125*tmpFldU(i,j,k,bi,bj)
209 & *rAw(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
210 tmpFldV(i,j,k,bi,bj) = -0.125*tmpFldV(i,j,k,bi,bj)
211 & *rAs(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
212 ENDDO
213 ENDDO
214 ELSE
215 DO J=1,sNy+1
216 DO I=1,sNx+1
217 tmpFldU(i,j,k,bi,bj) = -0.125*tmpFldU(i,j,k,bi,bj)
218 & *Shap_uvLength*Shap_uvLength*_maskW(i,j,k,bi,bj)
219 tmpFldV(i,j,k,bi,bj) = -0.125*tmpFldV(i,j,k,bi,bj)
220 & *Shap_uvLength*Shap_uvLength*_maskS(i,j,k,bi,bj)
221 ENDDO
222 ENDDO
223 ENDIF
224 ELSE
225 DO J=1,sNy
226 DO I=1,sNx+1
227 tmpFldU(i,j,k,bi,bj) = -0.125*
228 & ( hDiv(i,j)-hDiv(i-1,j)
229 & -vort3(i,j+1)+vort3(i,j)
230 & )*maskW(i,j,k,bi,bj)
231 ENDDO
232 ENDDO
233 DO J=1,sNy+1
234 DO I=1,sNx
235 tmpFldV(i,j,k,bi,bj) = -0.125*
236 & ( vort3(i+1,j)-vort3(i,j)
237 & +hDiv(i,j)-hDiv(i,j-1)
238 & )*maskS(i,j,k,bi,bj)
239 ENDDO
240 ENDDO
241
242 ENDIF
243
244 ENDDO
245 ENDDO
246 ENDDO
247 C end loop N=1,nShapUV
248 ENDDO
249
250 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
251
252 C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
253 DO bj=myByLo(myThid),myByHi(myThid)
254 DO bi=myBxLo(myThid),myBxHi(myThid)
255 DO K=1,kSize
256 DO J=1,sNy+1
257 DO I=1,sNx
258 uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
259 & -tmpFldU(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
260 ENDDO
261 ENDDO
262 DO J=1,sNy+1
263 DO I=1,sNx
264 vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
265 & -tmpFldV(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
266 ENDDO
267 ENDDO
268 ENDDO
269 ENDDO
270 ENDDO
271
272 IF (kSize.EQ.Nr) THEN
273 CALL EXCH_UV_XYZ_RL(uFld,vFld,.TRUE.,myThid)
274 ELSEIF (kSize.EQ.1) THEN
275 CALL EXCH_UV_XY_RL(uFld,vFld,.TRUE.,myThid)
276 ELSE
277 STOP 'S/R SHAP_FILT_UV_S2: kSize is wrong'
278 ENDIF
279
280 ENDIF
281 #endif /* ALLOW_SHAP_FILT */
282
283 RETURN
284 END
  ViewVC Help
Powered by ViewVC 1.1.22