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