MITgcm |
|
|
|
|||
File indexing completed on 2018-03-02 18:44:41 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC42c525bfb4 Alis*0001 #include "ZONAL_FILT_OPTIONS.h" 0002 0003 SUBROUTINE ZONAL_FILT_PRESMOOTH( 0004 I holeMask, 0005 U field, 0006 O avgField, 0007 I lField, 0008 I myThid ) 0009 C /==========================================================\ 0010 C | S/R ZONAL_FILT_PRESMOOTH | 0011 C | o Smooth data with holes ready for FFT. | 0012 C |==========================================================| 0013 C | FFT routines act on a series of data points. Having | 0014 C | arbitrary values in land introduces unrealistic noise in | 0015 C | the series. A better approach is to linearly ramp data | 0016 C | through the missing points. The mean of the field is also| 0017 C | removed. The mean is restored in FFT_POSTSMOOTH. This | 0018 C | step ensures no bias is introduced by the FFT process - | 0019 C | strictly it isnt necessary, but it can help improve | 0020 C | numerical conditioning. | 0021 C \==========================================================/ 0022 IMPLICIT NONE 0023 0024 C == Global data == 0025 #include "SIZE.h" 0026 #include "EEPARAMS.h" 0027 #include "PARAMS.h" 0028 #include "GRID.h" 0029 0030 C == Routine arguments == 0031 C holeMask - Array with 0 for holes and != 0 for valid data. 0032 C lField - Length of field to smooth (assumed periodic) 0033 C field - Field smoothed. 0034 C myThid - Thread number of this instance of FFT_PRESMOOTH_IN_X 0035 INTEGER lField 0036 Real*8 holeMask(lField) 0037 Real*8 field(lField) 0038 Real*8 avgField 0039 INTEGER myThid 0040 0041 #ifdef ALLOW_ZONAL_FILT 0042 0043 C == Local variables ==== 0044 C I - Loop counter 0045 C lbuf - Size of buffer arrays 0046 C hBase - Coord for last valid data point. 0047 C hHead - Coord for next valid data point. 0048 C hLo - Index for last valid data point. 0049 C hHi - Index for next valid data point. 0050 C nValid - Count of valid data points. 0051 C dist, len - Temps for interpolating. 0052 C frac Interpolation is simple linear ramp 0053 C between end-point in grid-point space. 0054 C e.g. for a series of points 0055 C Index 1 2 3 4 5 6 7 8 9 10 11 12 0056 C Data M V M M M M V V V V V V 0057 C where M indicates missing data and V valid 0058 C and 1 2 3 .... are indexes in field. Values 0059 C for the M points are found by interpolating 0060 C between the two V values that bracket a series 0061 C of M point. For index=1 0062 C dist=1, len=2 and frac=dist/len=1/2 so that 0063 C the interpolated value at index=1 is 0064 C V(index=12)+frac*( V(index=2)-V(index=12) ). 0065 C For index=5 dist=3, len=5 and frac=dist/len=3/5 0066 C so interpolated value at index=5 is 0067 C V(index=2)+frac*( V(index=7)-V(index=2) ). 0068 C lastGood - Temp for tracking of last good data point index. 0069 C nValid - Temp for counting number of valid points. 0070 C 0071 INTEGER lBuf 0072 PARAMETER ( lBuf = sNx ) 0073 INTEGER hBase(lBuf) 0074 INTEGER hHead(lBuf) 0075 INTEGER hLo(lBuf) 0076 INTEGER hHi(lBuf) 0077 INTEGER lastGood 0078 _RL dist 0079 _RL len 0080 _RL frac 0081 INTEGER nValid 0082 INTEGER I, iLo, iHi 0083 0084 0085 C 0086 IF ( lField .GT. lBuf ) THEN 0087 STOP 'S/R FFT_PRESMOOTH_1D: lField .GT. lBuf' 0088 ENDIF 0089 0090 CcnhDebugStarts 2cfc9d59a2 Patr*0091 C WRITE(*,*) 'Before FFT pre-smooth' 0092 C WRITE(*,'(A,A)') 'Mask ', 'Data' 42c525bfb4 Alis*0093 C DO I=1,lField 2cfc9d59a2 Patr*0094 C WRITE(*,'(A,I4,1X,F3.1,1X,1PE35.25)') 'I= ',I,holeMask(I), field(I) 42c525bfb4 Alis*0095 C ENDDO 0096 CcnhDebugEnds 0097 0098 C Count number of valid data points 0099 nValid = 0 0100 avgField = 0. 0101 DO I=1,lField 0102 IF ( holeMask(I) .NE. 0. ) THEN 0103 nValid = nValid+1 0104 avgField = avgField+field(I) 0105 ENDIF 0106 ENDDO 0107 0108 IF ( lField .GT. 1 .AND. nValid .GT. 0 ) THEN 0109 C Get lists of hole starts, ends and extents 0110 C ( use periodic wrap around ). 0111 0112 C 1. hLo - Index of valid point at start of run of holes 0113 C hBase - Coord of hLo (used to get offset when interpolating). 0114 C Note: The mean is also subtracted from field here. 0115 lastGood = -1 0116 avgField = avgField/FLOAT(nValid) 0117 DO I=1,lField 0118 IF ( holeMask(I) .EQ. 0. ) THEN 0119 C A hole 0120 hLo (I) = lastGood 0121 hBase(I) = lastGood 0122 ELSE 0123 C Data 0124 hLo(I) = 0 0125 hBase(I) = 0 0126 lastGood = I 0127 field(I) = field(I)-avgField 0128 ENDIF 0129 ENDDO 0130 DO I=1,lField 0131 IF ( hLo(I) .EQ. -1 ) THEN 0132 hLo(I) = lastGood 0133 hBase(I) = lastGood-lField 0134 ENDIF 0135 ENDDO 0136 0137 C 2. hHi - Coord of valid point at end of run of holes. 0138 lastGood = -1 0139 DO I=lField,1,-1 0140 IF ( holeMask(I) .EQ. 0. ) THEN 0141 C A hole 0142 hHi(I) = lastGood 0143 hHead(I) = lastGood 0144 ELSE 0145 C Data 0146 hHi(I) = 0 0147 lastGood = I 0148 ENDIF 0149 ENDDO 0150 DO I=lField,1,-1 0151 IF ( hHi(I) .EQ. -1 ) THEN 0152 hHi(I) = lastGood 0153 hHead(I) = lastGood+lField 0154 ENDIF 0155 ENDDO 0156 0157 CcnhDebugStarts 2cfc9d59a2 Patr*0158 C WRITE(*,*) 'During FFT pre-smooth' 0159 C WRITE(*,'(A,A,A,A,A,A)') 'I ','mask(I)','hHi(I)','hLo(I)','hBase(I)','hHead(I)' 42c525bfb4 Alis*0160 C DO I=1,lField 2cfc9d59a2 Patr*0161 C WRITE(*,'(6(I4,1X))') 42c525bfb4 Alis*0162 C & I,INT(holeMask(I)),hHi(I),hLo(I),hBase(I),hHead(I) 0163 C ENDDO 0164 CcnhDebugEnds 0165 0166 C 3. Interpolate 0167 DO I=1,lField 0168 IF ( holeMask(I) .EQ. 0. ) THEN 0169 C A hole 0170 iLo = hLo(I) 0171 iHi = hHi(I) 0172 dist = I-hBase(I) 0173 len = hHead(I) - hBase(I) 0174 CcnhDebugStarts 2cfc9d59a2 Patr*0175 C WRITE(*,'(A,1X,I4,1X,1PE35.25,1X,1PE35.25,)') 'I= ',I,dist, len 42c525bfb4 Alis*0176 C IF ( dist .LT. 0 ) STOP ' DIST .LT. 0 ' 0177 C IF ( len .LT. 0 ) STOP ' LEN .LT. 0 ' 0178 C IF ( dist .GT. lField ) STOP ' dist .GT. lField ' 0179 C IF ( len .GT. lField ) STOP ' len .GT. lField ' 0180 C IF ( dist .GT. len ) STOP ' dist .GT. len ' 0181 CcnhDebugStarts 0182 frac = dist/len 0183 field(I) = field(iLo) 0184 & +(field(iHi)-field(iLo))*frac 0185 ENDIF 0186 ENDDO 0187 0188 CcnhDebugStarts 2cfc9d59a2 Patr*0189 C WRITE(*,*) 'After FFT pre-smooth' 0190 C WRITE(*,'(A,A)') 'Mask ', 'Data' 42c525bfb4 Alis*0191 C DO I=1,lField 2cfc9d59a2 Patr*0192 C WRITE(*,'(A,I4,1X,F3.1,1X,1PE35.25)') 'I= ',I,holeMask(I), field(I) 42c525bfb4 Alis*0193 C ENDDO 0194 CcnhDebugEnds 0195 0196 ENDIF 0197 C 0198 #endif /* ALLOW_ZONAL_FILT */ 0199 0200 RETURN 0201 END
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine. The LXR team |
|