SUBROUTINE FILTR( 12,2@DYALLOC.4151
*CALL ARGSIZE 
@DYALLOC.4152
*CALL ARGOCFIL @DYALLOC.4153
& FTARR, ORH1F304.105
* S,IM,MM,N,ISS) ! ############################################### @DYALLOC.4154
C OFILTR.4
C======================================================================= OFILTR.5
C === OFILTR.6
C FILTER FOURIER ANALYSES THE ARRAYS OF VARIOUS === OFILTR.7
C PHYSICAL QUANTITIES, THEN TRUNCATES THE SERIES AND === OFILTR.8
C RESYNTHESIZES THE FILTERED QUANTITIES WHERE: === OFILTR.9
C S =THE STRING TO BE FILTERED === OFILTR.10
C IM =THE LENGTH OF S === OFILTR.11
C MM =1 (COSINE SERIES, DERIV AT BNDRY PTS=0) === OFILTR.12
C =2 (SINE SERIES, BNDRY PTS=0) === OFILTR.13
C =3 (FULL SERIES, CYCLIC) === OFILTR.14
C N =NUMBER OF WAVES TO KEEP === OFILTR.15
C ISS=0 (CANT USE FOURIER COEFS FROM PREVIOUS CALL) === OFILTR.16
C ISS>0 (CAN USE FOURIER COEFS FROM PREVIOUS CALL) === OFILTR.17
CLL History: ORH1F304.100
CLL Version Date Comment ORH1F304.101
CLL ------- ---- ------- ORH1F304.102
CLL 3.4 01.09.94 put array FTARR in argument list explicitly ORH1F304.103
CLL R. Hill. ORH1F304.104
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.4233
C === OFILTR.18
C======================================================================= OFILTR.19
C OFILTR.20
C--------------------------------------------------------------------- OFILTR.21
C DEFINE GLOBAL DATA OFILTR.22
C--------------------------------------------------------------------- OFILTR.23
C OFILTR.24
*CALL OARRYSIZ ORH6F401.21
*CALL TYPSIZE @DYALLOC.4155
*CALL TYPOCFIL @DYALLOC.4156
*CALL CNTLOCN ORH1F305.4234
*CALL OTIMER ORH1F305.4236
C OFILTR.26
*CALL COCNINDX ORH1F403.147
C--------------------------------------------------------------------- OFILTR.27
C DEFINE LOCAL DATA AND DIMENSION ARGUMENT ARRAYS OFILTR.28
C--------------------------------------------------------------------- OFILTR.29
C OFILTR.30
C COSSAV MUST REMAIN FULL PRECISION IF MOST OF FILTER IS MADE HALF-P @DYALLOC.4157
C OFILTR.35
REAL FTARR(IMTIMT) ! Coef used in filtering routine ORH1F304.106
DIMENSION TEMP(IMT*4) @DYALLOC.4158
DIMENSION CIRCLE(4) OFILTR.41
DIMENSION INDX(IMT*8),COF(IMT*8) @DYALLOC.4159
DIMENSION COSINE(IMT*8) @DYALLOC.4160
DIMENSION DENOM(IMT*4) @DYALLOC.4161
DIMENSION S(IMT),SPRIME(IMT) OFILTR.45
DIMENSION DIFF(IM+1*2) @DYALLOC.4162
C OFILTR.51
DATA PI/3.141592653589793/, CIRCLE/0.,-1.,0.,1./ OFILTR.52
C OFILTR.59
C--------------------------------------------------------------------- OFILTR.60
C BEGIN EXECUTABLE CODE OFILTR.61
C--------------------------------------------------------------------- OFILTR.62
C OFILTR.63
IF (L_OTIMER) THEN ORH1F305.4237
CALL TIMER('FILTR ',103) GPB8F405.109
ENDIF ORH1F305.4239
IF (IM.LT.1 .OR. MM.LT.1 .OR. MM.GT.3 .OR. N.LT.0 .OR. ISS.LT.0) OFILTR.69
* GOTO 4000 OFILTR.70
C OFILTR.71
IF (INITDN) GOTO 90 OFILTR.72
C OFILTR.73
C THIS SECTION SETS UP TABLES FOR FILTER; IT MUST BE CALLED ONCE PER OFILTR.74
C EXECUTION OF OCEAN OFILTR.75
C OFILTR.76
C NOTE: LQMSUM IS THE SUM OF (IM-1)/2 FOR IM=1,IMTP1 OFILTR.77
C LHSUM IS THE SUM OF IM-1 FOR IM=1,IMTP1 OFILTR.78
C OFILTR.79
IMSAVE=IM OFILTR.80
C OFILTR.81
C ASSEMBLE INDEX ARRAY OFILTR.82
C OFILTR.83
DO 10 I=1,IMTX8 OFILTR.84
IND(I)=I OFILTR.85
10 CONTINUE OFILTR.86
C OFILTR.87
C CALCULATE AND SAVE ALL COSINES WHICH WILL BE NEEDED OFILTR.88
C OFILTR.89
IBASE=0 OFILTR.90
JBASE=0 OFILTR.91
DO 50 IM=1,IMTP1 OFILTR.92
FIMR=1.0/FLOAT(IM) OFILTR.93
C OFILTR.94
IMM1=IM-1 OFILTR.95
IF (IMM1.EQ.0) GOTO 25 OFILTR.96
DO 20 I=1,IMM1 OFILTR.97
DENMSV(IBASE+I)=1.0/(1.0-COS(PI*FLOAT(I)*FIMR)) OFILTR.98
20 CONTINUE OFILTR.99
25 IDBASE(IM)=IBASE OFILTR.100
IBASE=IBASE+IMM1 OFILTR.101
C OFILTR.102
IMQC=(IM-1)/2 OFILTR.103
IF (IMQC.EQ.0) GOTO 35 OFILTR.104
DO 30 I=1,IMQC OFILTR.105
COSSAV(JBASE+I)=COS(PI*FLOAT(I)*FIMR) OFILTR.106
30 CONTINUE OFILTR.107
35 ICBASE(IM)=JBASE OFILTR.108
JBASE=JBASE+IMQC OFILTR.109
C OFILTR.110
50 CONTINUE OFILTR.111
C OFILTR.112
C CALCULATE ADJUSTMENTS FOR GENERAL FOURIER CASE IF IM=2*N OFILTR.113
C OFILTR.114
DO 60 IM=1,IMT OFILTR.115
COSNPI(IM)=CIRCLE(MOD(IM-1,4)+1) OFILTR.116
60 CONTINUE OFILTR.117
INITDN=.TRUE. OFILTR.118
IM=IMSAVE OFILTR.119
C OFILTR.120
C OFILTR.121
C CALCULATE SOME USEFUL CONSTANTS OFILTR.122
C OFILTR.123
90 IF(MM.EQ.2 .AND. N.EQ.0) THEN OFILTR.124
DO 92 I=1,IM OFILTR.125
S(I)=0.0 OFILTR.126
92 CONTINUE OFILTR.127
GO TO 3950 OFILTR.128
ENDIF OFILTR.129
IF (MM.EQ.1) THEN OFILTR.130
NMAX=N-1 OFILTR.131
ELSE OFILTR.132
NMAX=N OFILTR.133
ENDIF OFILTR.134
NMAXP1=NMAX+1 OFILTR.135
C1=0.5*FLOAT(NMAX)+0.25 OFILTR.136
C2=FLOAT(NMAX)+0.5 OFILTR.137
C OFILTR.138
IF (MM.EQ.2) THEN OFILTR.139
LCY=2*(IM+1) OFILTR.140
FNORM=2.0/FLOAT(IM+1) OFILTR.141
ELSE OFILTR.142
LCY=2*IM OFILTR.143
FNORM=2.0/FLOAT(IM) OFILTR.144
ENDIF OFILTR.145
LH=LCY/2 OFILTR.146
LHM1=LH-1 OFILTR.147
LQM=(LH-1)/2 OFILTR.148
L2CY=2*LCY OFILTR.149
LCYM1=LCY-1 OFILTR.150
LCYP1=LCY+1 OFILTR.151
C OFILTR.152
IMX2=IM*2 OFILTR.153
IMX4=IM*4 OFILTR.154
IMX8=IM*8 OFILTR.155
C OFILTR.156
C AVERAGE INCOMING ARRAY OFILTR.157
C OFILTR.158
SSUM=0.0 OFILTR.159
DO 100 I=1,IM OFILTR.160
100 SSUM=SSUM + S(I) OFILTR.161
C OFILTR.162
C MM = 1 DERIVATIVE MUST BE ZERO AT BOUNDARIES (COSINE) OFILTR.163
C MM = 2 VALUE MUST BE ZERO AT BOUNDARIES (SINE) OFILTR.164
C MM = 3 CYCLIC BOUNDARY CONDITIONS (GENERAL FOURIER SERIES) OFILTR.165
C OFILTR.166
FIM=FLOAT(IM) OFILTR.167
FIMR=1.0/FIM OFILTR.168
STEMP=SSUM*FIMR OFILTR.169
IF (N.GT.1 .OR. MM.NE.1) GO TO 400 OFILTR.170
DO 300 I=1,IM OFILTR.171
300 S(I)=STEMP OFILTR.172
GO TO 3950 OFILTR.173
400 CONTINUE OFILTR.174
IF(MM.NE.2) THEN OFILTR.175
DO 450 I=1,IM OFILTR.176
S(I)=S(I)-STEMP OFILTR.177
450 CONTINUE OFILTR.178
ENDIF OFILTR.179
IF (ISS.GT.0) GO TO 3000 OFILTR.180
C OFILTR.181
C OFILTR.182
C ASSEMBLE APPROPRIATE 1-CYCLE (2*PI) COSINE ARRAY OFILTR.183
C OFILTR.184
C USE STORED 1/4 CYCLE TO CALCULATE FIRST 1/2 CYCLE OFILTR.185
JBASE=ICBASE(LH) OFILTR.186
DO 700 I=1,LQM OFILTR.187
COSINE(I)=COSSAV(JBASE+I) OFILTR.188
700 CONTINUE OFILTR.189
DO 701 I=1,LQM OFILTR.190
701 COSINE(LH-I)=-COSSAV(JBASE+I) OFILTR.191
C FILL IN COS(PI/2) IF LH IS EVEN OFILTR.192
IF (2*(LQM+1).EQ.LH) COSINE(LQM+1)=0.0 OFILTR.193
C FILL IN COS(PI) IN ANY CASE OFILTR.194
COSINE(LH)=-1.0 OFILTR.195
C FILL IN REST OF CYCLE OFILTR.196
DO 710 I=1,LH OFILTR.197
COSINE(LH+I)=-COSINE(I) OFILTR.198
710 CONTINUE OFILTR.199
C OFILTR.200
C ASSEMBLE DENOMINATOR ARRAY OFILTR.201
C OFILTR.202
IBASE=IDBASE(LH) OFILTR.203
FXA=0.25 OFILTR.204
DO 720 I=1,LHM1 OFILTR.205
DENOM(I)=FXA*DENMSV(IBASE+I) OFILTR.206
720 CONTINUE OFILTR.207
DENOM(LH)=0.125 OFILTR.208
DO 721 I=1,LHM1 OFILTR.209
TEMP(I)=DENOM(LH-I) OFILTR.210
721 CONTINUE OFILTR.211
DO 722 I=1,LHM1 OFILTR.212
722 DENOM(LH+I)=TEMP(I) OFILTR.213
NPRINT=0 OFILTR.214
DENOM(LCY)=0.0 OFILTR.215
DO 730 I=LCYP1,IMX4 OFILTR.216
DENOM(I)=DENOM(I-LCY) OFILTR.217
730 CONTINUE OFILTR.218
C OFILTR.219
C ASSEMBLE APPROPRIATE SUBSCRIPT ARRAYS OFILTR.220
C OFILTR.221
C CALCULATE NEEDED INDICES OFILTR.222
C OFILTR.223
IF (MM.EQ.3) THEN OFILTR.224
FACT1=2*NMAX OFILTR.225
FACT2=2*NMAXP1 OFILTR.226
ELSE OFILTR.227
FACT1=NMAX OFILTR.228
FACT2=NMAXP1 OFILTR.229
ENDIF OFILTR.230
DO 740 I=1,IMX4 OFILTR.231
INDX(I)=IND(I)*FACT1 OFILTR.232
740 CONTINUE OFILTR.233
DO 741 I=1,IMX4 OFILTR.234
741 INDX(IMX4+I)=IND(I)*FACT2 OFILTR.235
C CALCULATE PARAMETERS FOR REDUCING INDICES OFILTR.236
MAXIND=IMX4*FACT2 OFILTR.237
NCYC=(MAXIND-1)/LCY + 1 OFILTR.238
MAXNDX=LCY OFILTR.239
IF (MAXNDX.GE.MAXIND) GOTO 790 OFILTR.240
DO 750 NPWR=1,NCYC+2 OFILTR.241
MAXNDX=2*MAXNDX OFILTR.242
IF (MAXNDX.GE.MAXIND) GOTO 760 OFILTR.243
750 CONTINUE OFILTR.244
STOP 'ERROR: FELL THROUGH DO-LOOP TERMINATION AT 750' OFILTR.245
760 DO 770 NP=1,NPWR OFILTR.246
MAXNDX=MAXNDX/2 OFILTR.247
DO 765 I=1,IMX8 OFILTR.248
IF(INDX(I).GT.MAXNDX) INDX(I)=INDX(I)-MAXNDX OFILTR.249
765 CONTINUE OFILTR.250
770 CONTINUE OFILTR.251
790 CONTINUE OFILTR.252
C OFILTR.253
C GATHER COEFFICIENTS OFILTR.254
C OFILTR.255
800 DO 810 J=1,IMX8 OFILTR.256
COF(J)=COSINE(INDX(J)) OFILTR.257
810 CONTINUE OFILTR.258
C OFILTR.259
C OFILTR.260
C ASSEMBLE TRANSFORMATION ARRAY WHICH WILL FILTER S OFILTR.261
C OFILTR.262
C OFILTR.263
IF(MM.EQ.1)THEN OFILTR.264
C COSINE TRANSFORM OFILTR.265
IOFF1=LCY OFILTR.266
IOFF2=LCY+IMX4 OFILTR.267
FXA=0.5 OFILTR.268
DO 1200 J=1,IM OFILTR.269
JOFF=(J-1)*IMT OFILTR.270
DO 1100 I=1,IM OFILTR.271
FTARR(JOFF+I)=(COF(I-J+IOFF1)-COF(I-J+IOFF2))*DENOM(I-J+IOFF1) OFILTR.272
* +(COF(I+J-1)-COF(IMX4+I+J-1))*DENOM(I+J-1) - FXA OFILTR.273
1100 CONTINUE OFILTR.274
1200 CONTINUE OFILTR.275
DO 1201 J=1,IM OFILTR.276
1201 FTARR(J*IMTP1-IMT)=FTARR(J*IMTP1-IMT)+C1 OFILTR.277
ELSE IF(MM.EQ.2) THEN OFILTR.278
C OFILTR.279
C SINE TRANSFORM OFILTR.280
IOFF1=LCY OFILTR.281
IOFF2=LCY+IMX4 OFILTR.282
DO 1500 J=1,IM OFILTR.283
JOFF=(J-1)*IMT OFILTR.284
DO 1400 I=1,IM OFILTR.285
FTARR(JOFF+I)=(COF(I-J+IOFF1)-COF(I-J+IOFF2))*DENOM(I-J+IOFF1) OFILTR.286
* -(COF(I+J)-COF(IMX4+I+J))*DENOM(I+J) OFILTR.287
1400 CONTINUE OFILTR.288
1500 CONTINUE OFILTR.289
DO 1501 J=1,IM OFILTR.290
1501 FTARR(J*IMTP1-IMT)=FTARR(J*IMTP1-IMT)+C1 OFILTR.291
ELSE IF(MM.EQ.3) THEN OFILTR.292
C OFILTR.293
C GENERAL FOURIER TRANSFORM OFILTR.294
IF(2*N.EQ.IM)THEN OFILTR.295
GENADJ=0.5 OFILTR.296
ELSE OFILTR.297
GENADJ=0.0 OFILTR.298
ENDIF OFILTR.299
IOFF1=LCY OFILTR.300
IOFF2=LCY+IMX4 OFILTR.301
FXA=2.0 OFILTR.302
FXB=0.5 OFILTR.303
DO 1800 J=1,IM OFILTR.304
JOFF=(J-1)*IMT OFILTR.305
DO 1700 I=1,IM OFILTR.306
FTARR(JOFF+I)= OFILTR.307
* (FXA*(COF(I-J+IOFF1)-COF(I-J+IOFF2)))*DENOM(2*I-2*J+IOFF1) OFILTR.308
* -FXB - GENADJ*COSNPI(I)*COSNPI(J) OFILTR.309
1700 CONTINUE OFILTR.310
1800 CONTINUE OFILTR.311
DO 1801 J=1,IM OFILTR.312
1801 FTARR(J*IMTP1-IMT)=FTARR(J*IMTP1-IMT)+C2 OFILTR.313
ENDIF OFILTR.314
C OFILTR.315
C FILTER S OFILTR.316
C OFILTR.317
3000 CONTINUE @DYALLOC.4164
DO 3010 I=1,IM @DYALLOC.4165
SPRIME(I)=0. OFILTR.319
3010 CONTINUE OFILTR.320
DO 3100 I=1,IM OFILTR.321
IOFF=(I-1)*IMT OFILTR.322
DO 3100 J=1,IM OFILTR.323
C NOTE THAT FTARR(J,I)=FTARR(I,J), SO FOLLOWING IS LEGAL OFILTR.324
SPRIME(J)=SPRIME(J)+S(I)*FTARR(IOFF+J) OFILTR.325
3100 CONTINUE OFILTR.326
DO 3110 I=1,IM OFILTR.327
SPRIME(I)=FNORM*SPRIME(I) OFILTR.328
3110 CONTINUE OFILTR.329
IF(MM.EQ.2) THEN OFILTR.330
DO 3150 I=1,IM OFILTR.331
S(I)=SPRIME(I) OFILTR.332
3150 CONTINUE OFILTR.333
GO TO 3950 OFILTR.334
ENDIF OFILTR.335
C OFILTR.336
3700 SSM=0.0 OFILTR.337
DO 3800 I=1,IM OFILTR.338
SSM=SSM+SPRIME(I) OFILTR.339
3800 CONTINUE OFILTR.340
SSM=(SSUM-SSM)*FIMR OFILTR.341
DO 3900 I=1,IM OFILTR.342
S(I)=SSM+SPRIME(I) OFILTR.343
3900 CONTINUE OFILTR.344
3950 CONTINUE OFILTR.345
IF (L_OTIMER) THEN ORH1F305.4240
CALL TIMER('FILTR ',104) GPB8F405.110
ENDIF ORH1F305.4242
RETURN OFILTR.351
4000 PRINT 4001, IM,MM,N,ISS OFILTR.352
4001 FORMAT (' BAD ARGUMENT(S) IN CALL TO FILTR'/' IM,MM,N,ISS = ', OFILTR.353
* 4I10) OFILTR.354
STOP 'BAD ARGUMENT(S) IN CALL TO FILTR' OFILTR.355
END OFILTR.356
*ENDIF @DYALLOC.4167