-
Notifications
You must be signed in to change notification settings - Fork 237
/
shap_filt_tracer_s4.F
179 lines (148 loc) · 4.48 KB
/
shap_filt_tracer_s4.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
#include "SHAP_FILT_OPTIONS.h"
CBOP
C !ROUTINE: SHAP_FILT_TRACER_S4
C !INTERFACE:
SUBROUTINE SHAP_FILT_TRACER_S4(
U field, tmpFld,
I nShapTr, kSize, myTime, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R SHAP_FILT_TRACER_S4
C | o Applies Shapiro filter to tracer field (cell center).
C | o use filtering function "S4" = [1 - d_xx^n][1- d_yy^n]
C | with no grid spacing (computational Filter)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SHAP_FILT.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments
C field :: cell-centered 2D field on which filter applies
C tmpFld :: working temporary array
C nShapTr :: (total) power of the filter for this tracer
C kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
C myTime :: Current time in simulation
C myThid :: Thread number for this instance of SHAP_FILT_TRACER_S4
INTEGER nShapTr, kSize
_RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
_RL myTime
INTEGER myThid
#ifdef ALLOW_SHAP_FILT
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER bi,bj,K,I,J,N
_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
IF (nShapTr.gt.0) THEN
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,kSize
DO J=1-OLy,sNy+Oly
DO I=1-Olx,sNx+Olx
tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C d_xx^n tmpFld
DO N=1,nShapTr
IF (kSize.EQ.Nr) THEN
_EXCH_XYZ_RL( tmpFld, myThid )
ELSE
_EXCH_XY_RL( tmpFld, myThid )
ENDIF
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,kSize
DO J=1,sNy
DO I=1,sNx
tmpGrd(i,j) = -0.25*(
& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
& *_maskW(i+1,j,k,bi,bj)
& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
& *_maskW(i,j,k,bi,bj) )
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C F <- [1 - d_xx^n *deltaT/tau].F
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,kSize
DO J=1,sNy
DO I=1,sNx
field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C d_yy^n tmpFld
DO N=1,nShapTr
IF (kSize.EQ.1) THEN
_EXCH_XY_RL( tmpFld, myThid )
ELSE
_EXCH_XYZ_RL( tmpFld, myThid )
ENDIF
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,kSize
DO J=1,sNy
DO I=1,sNx
tmpGrd(i,j) = -0.25*(
& ( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
& *_maskS(i,j+1,k,bi,bj)
& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
& *_maskS(i,j,k,bi,bj) )
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C F <- [1 - d_yy^n *deltaT/tau].F
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,kSize
DO J=1,sNy
DO I=1,sNx
field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
IF (kSize.EQ.Nr) THEN
_EXCH_XYZ_RL( field, myThid )
ELSEIF (kSize.EQ.1) THEN
_EXCH_XY_RL( field, myThid )
ELSE
STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'
ENDIF
ENDIF
#endif /* ALLOW_SHAP_FILT */
RETURN
END