1 |
C $Header: $ |
2 |
|
3 |
#include "ZONAL_FILT_OPTIONS.h" |
4 |
|
5 |
SUBROUTINE ZONAL_FILT_PRESMOOTH( |
6 |
I holeMask, |
7 |
U field, |
8 |
O avgField, |
9 |
I lField, |
10 |
I myThid ) |
11 |
C /==========================================================\ |
12 |
C | S/R ZONAL_FILT_PRESMOOTH | |
13 |
C | o Smooth data with holes ready for FFT. | |
14 |
C |==========================================================| |
15 |
C | FFT routines act on a series of data points. Having | |
16 |
C | arbitrary values in land introduces unrealistic noise in | |
17 |
C | the series. A better approach is to linearly ramp data | |
18 |
C | through the missing points. The mean of the field is also| |
19 |
C | removed. The mean is restored in FFT_POSTSMOOTH. This | |
20 |
C | step ensures no bias is introduced by the FFT process - | |
21 |
C | strictly it isnt necessary, but it can help improve | |
22 |
C | numerical conditioning. | |
23 |
C \==========================================================/ |
24 |
IMPLICIT NONE |
25 |
|
26 |
C == Global data == |
27 |
#include "SIZE.h" |
28 |
#include "EEPARAMS.h" |
29 |
#include "PARAMS.h" |
30 |
#include "GRID.h" |
31 |
|
32 |
C == Routine arguments == |
33 |
C holeMask - Array with 0 for holes and != 0 for valid data. |
34 |
C lField - Length of field to smooth (assumed periodic) |
35 |
C field - Field smoothed. |
36 |
C myThid - Thread number of this instance of FFT_PRESMOOTH_IN_X |
37 |
INTEGER lField |
38 |
Real*8 holeMask(lField) |
39 |
Real*8 field(lField) |
40 |
Real*8 avgField |
41 |
INTEGER myThid |
42 |
|
43 |
#ifdef ALLOW_ZONAL_FILT |
44 |
|
45 |
C == Local variables ==== |
46 |
C I - Loop counter |
47 |
C lbuf - Size of buffer arrays |
48 |
C hBase - Coord for last valid data point. |
49 |
C hHead - Coord for next valid data point. |
50 |
C hLo - Index for last valid data point. |
51 |
C hHi - Index for next valid data point. |
52 |
C nValid - Count of valid data points. |
53 |
C dist, len - Temps for interpolating. |
54 |
C frac Interpolation is simple linear ramp |
55 |
C between end-point in grid-point space. |
56 |
C e.g. for a series of points |
57 |
C Index 1 2 3 4 5 6 7 8 9 10 11 12 |
58 |
C Data M V M M M M V V V V V V |
59 |
C where M indicates missing data and V valid |
60 |
C and 1 2 3 .... are indexes in field. Values |
61 |
C for the M points are found by interpolating |
62 |
C between the two V values that bracket a series |
63 |
C of M point. For index=1 |
64 |
C dist=1, len=2 and frac=dist/len=1/2 so that |
65 |
C the interpolated value at index=1 is |
66 |
C V(index=12)+frac*( V(index=2)-V(index=12) ). |
67 |
C For index=5 dist=3, len=5 and frac=dist/len=3/5 |
68 |
C so interpolated value at index=5 is |
69 |
C V(index=2)+frac*( V(index=7)-V(index=2) ). |
70 |
C lastGood - Temp for tracking of last good data point index. |
71 |
C nValid - Temp for counting number of valid points. |
72 |
C |
73 |
INTEGER lBuf |
74 |
PARAMETER ( lBuf = sNx ) |
75 |
INTEGER hBase(lBuf) |
76 |
INTEGER hHead(lBuf) |
77 |
INTEGER hLo(lBuf) |
78 |
INTEGER hHi(lBuf) |
79 |
INTEGER lastGood |
80 |
_RL dist |
81 |
_RL len |
82 |
_RL frac |
83 |
INTEGER nValid |
84 |
INTEGER I, iLo, iHi |
85 |
|
86 |
|
87 |
C |
88 |
IF ( lField .GT. lBuf ) THEN |
89 |
STOP 'S/R FFT_PRESMOOTH_1D: lField .GT. lBuf' |
90 |
ENDIF |
91 |
|
92 |
CcnhDebugStarts |
93 |
C WRITE(0,*) 'Before FFT pre-smooth' |
94 |
C WRITE(0,'(A,A)') 'Mask ', 'Data' |
95 |
C DO I=1,lField |
96 |
C WRITE(0,'(A,I4,1X,F3.1,1X,1PE35.25)') 'I= ',I,holeMask(I), field(I) |
97 |
C ENDDO |
98 |
CcnhDebugEnds |
99 |
|
100 |
C Count number of valid data points |
101 |
nValid = 0 |
102 |
avgField = 0. |
103 |
DO I=1,lField |
104 |
IF ( holeMask(I) .NE. 0. ) THEN |
105 |
nValid = nValid+1 |
106 |
avgField = avgField+field(I) |
107 |
ENDIF |
108 |
ENDDO |
109 |
|
110 |
IF ( lField .GT. 1 .AND. nValid .GT. 0 ) THEN |
111 |
C Get lists of hole starts, ends and extents |
112 |
C ( use periodic wrap around ). |
113 |
|
114 |
C 1. hLo - Index of valid point at start of run of holes |
115 |
C hBase - Coord of hLo (used to get offset when interpolating). |
116 |
C Note: The mean is also subtracted from field here. |
117 |
lastGood = -1 |
118 |
avgField = avgField/FLOAT(nValid) |
119 |
DO I=1,lField |
120 |
IF ( holeMask(I) .EQ. 0. ) THEN |
121 |
C A hole |
122 |
hLo (I) = lastGood |
123 |
hBase(I) = lastGood |
124 |
ELSE |
125 |
C Data |
126 |
hLo(I) = 0 |
127 |
hBase(I) = 0 |
128 |
lastGood = I |
129 |
field(I) = field(I)-avgField |
130 |
ENDIF |
131 |
ENDDO |
132 |
DO I=1,lField |
133 |
IF ( hLo(I) .EQ. -1 ) THEN |
134 |
hLo(I) = lastGood |
135 |
hBase(I) = lastGood-lField |
136 |
ENDIF |
137 |
ENDDO |
138 |
|
139 |
C 2. hHi - Coord of valid point at end of run of holes. |
140 |
lastGood = -1 |
141 |
DO I=lField,1,-1 |
142 |
IF ( holeMask(I) .EQ. 0. ) THEN |
143 |
C A hole |
144 |
hHi(I) = lastGood |
145 |
hHead(I) = lastGood |
146 |
ELSE |
147 |
C Data |
148 |
hHi(I) = 0 |
149 |
lastGood = I |
150 |
ENDIF |
151 |
ENDDO |
152 |
DO I=lField,1,-1 |
153 |
IF ( hHi(I) .EQ. -1 ) THEN |
154 |
hHi(I) = lastGood |
155 |
hHead(I) = lastGood+lField |
156 |
ENDIF |
157 |
ENDDO |
158 |
|
159 |
CcnhDebugStarts |
160 |
C WRITE(0,*) 'During FFT pre-smooth' |
161 |
C WRITE(0,'(A,A,A,A,A,A)') 'I ','mask(I)','hHi(I)','hLo(I)','hBase(I)','hHead(I)' |
162 |
C DO I=1,lField |
163 |
C WRITE(0,'(6(I4,1X))') |
164 |
C & I,INT(holeMask(I)),hHi(I),hLo(I),hBase(I),hHead(I) |
165 |
C ENDDO |
166 |
CcnhDebugEnds |
167 |
|
168 |
C 3. Interpolate |
169 |
DO I=1,lField |
170 |
IF ( holeMask(I) .EQ. 0. ) THEN |
171 |
C A hole |
172 |
iLo = hLo(I) |
173 |
iHi = hHi(I) |
174 |
dist = I-hBase(I) |
175 |
len = hHead(I) - hBase(I) |
176 |
CcnhDebugStarts |
177 |
C WRITE(0,'(A,1X,I4,1X,1PE35.25,1X,1PE35.25,)') 'I= ',I,dist, len |
178 |
C IF ( dist .LT. 0 ) STOP ' DIST .LT. 0 ' |
179 |
C IF ( len .LT. 0 ) STOP ' LEN .LT. 0 ' |
180 |
C IF ( dist .GT. lField ) STOP ' dist .GT. lField ' |
181 |
C IF ( len .GT. lField ) STOP ' len .GT. lField ' |
182 |
C IF ( dist .GT. len ) STOP ' dist .GT. len ' |
183 |
CcnhDebugStarts |
184 |
frac = dist/len |
185 |
field(I) = field(iLo) |
186 |
& +(field(iHi)-field(iLo))*frac |
187 |
ENDIF |
188 |
ENDDO |
189 |
|
190 |
CcnhDebugStarts |
191 |
C WRITE(0,*) 'After FFT pre-smooth' |
192 |
C WRITE(0,'(A,A)') 'Mask ', 'Data' |
193 |
C DO I=1,lField |
194 |
C WRITE(0,'(A,I4,1X,F3.1,1X,1PE35.25)') 'I= ',I,holeMask(I), field(I) |
195 |
C ENDDO |
196 |
CcnhDebugEnds |
197 |
|
198 |
ENDIF |
199 |
C |
200 |
#endif /* ALLOW_ZONAL_FILT */ |
201 |
|
202 |
RETURN |
203 |
END |