SUBROUTINE SET_FIL 1,4SETFIL1A.36
1 (U,C,ADJUSTMENT_TIMESTEP,ADVECTION_TIMESTEP, SETFIL1A.37
2 SEC_P_LATITUDE,FILTER_WAVE_NUMBER_P_ROWS, SETFIL1A.38
3 FILTER_WAVE_NUMBER_U_ROWS, SETFIL1A.39
4 LONGITUDE_STEP_INVERSE,NORTHERN_FILTERED_P_ROW, SETFIL1A.40
5 SOUTHERN_FILTERED_P_ROW,P_FIELD,U_FIELD, SETFIL1A.41
6 P_LEVELS,ROW_LENGTH, APB0F402.97
*CALL ARGFLDPT 
APB0F402.98
7 FILTERING_SAFETY_FACTOR,TWO_D_GRID_CORRECTION) APB0F402.99
SETFIL1A.45
IMPLICIT NONE SETFIL1A.46
SETFIL1A.47
! All FLDPTR arguments are intent IN APB0F402.100
*CALL TYPFLDPT APB0F402.101
APB0F402.102
INTEGER SETFIL1A.48
* P_FIELD !IN DIMENSION OF FIELDS ON PRESSURE GRID SETFIL1A.49
*, U_FIELD !IN DIMENSION OF FIELDS ON VELOCITY GRID SETFIL1A.50
*, P_LEVELS !IN NUMBER OF MODEL LEVELS. SETFIL1A.51
*, ROW_LENGTH !IN NUMBER OF POINTS ON A ROW. SETFIL1A.52
*, NORTHERN_FILTERED_P_ROW !IN LAST ROW, MOVING EQUATORWARDS, SETFIL1A.53
* ! IN NORTHERN HEMISPHERE SETFIL1A.54
* ! ON WHICH FILTERING IS PERFORMED. SETFIL1A.55
*, SOUTHERN_FILTERED_P_ROW !IN LAST ROW, MOVING EQUATORWARDS, SETFIL1A.56
* ! IN SOUTHERN HEMISPHERE SETFIL1A.57
* ! ON WHICH FILTERING IS PERFORMED. SETFIL1A.58
SETFIL1A.59
INTEGER SETFIL1A.60
& FILTER_WAVE_NUMBER_P_ROWS(GLOBAL_P_FIELD/GLOBAL_ROW_LENGTH) APB0F402.103
& ! LARGEST WAVE NUMBER ON EACH P ROW WHICH IS NOT FILTERED APB0F402.104
&, FILTER_WAVE_NUMBER_U_ROWS(GLOBAL_U_FIELD/GLOBAL_ROW_LENGTH) APB0F402.105
& ! LARGEST WAVE NUMBER ON EACH U ROW WHICH IS NOT FILTERED APB0F402.106
SETFIL1A.69
SETFIL1A.73
REAL SETFIL1A.74
* U(U_FIELD,P_LEVELS) !IN ZONAL WIND COMPONENT. SETFIL1A.75
*,SEC_P_LATITUDE(P_FIELD) !IN 1/(COS(LAT)) AT P POINTS. SETFIL1A.76
*,ADJUSTMENT_TIMESTEP !IN SETFIL1A.77
*,ADVECTION_TIMESTEP !IN SETFIL1A.78
*,C !IN EXTENAL GRAVITY WAVE SPEED. SETFIL1A.79
*,LONGITUDE_STEP_INVERSE !IN 1/(DELTA LAMDA) SETFIL1A.80
*,TWO_D_GRID_CORRECTION(P_FIELD/ROW_LENGTH) !IN TERM REPRESENTING SETFIL1A.81
* ! 2-D EFFECT OF WAVES COMPARED TO 1D SETFIL1A.82
* ! REPRESENTATION USED HERE. SETFIL1A.83
*,FILTERING_SAFETY_FACTOR !IN TERM REPRESENTING SAFETY MARGIN TO SETFIL1A.84
* ! ADD TO GRID_CORRECTION SETFIL1A.85
C*--------------------------------------------------------------------- SETFIL1A.86
*IF DEF,MPP APB7F401.144
! Common blocks and parameters for MPP code APB7F401.145
*CALL PARVARS APB7F401.146
*CALL PARFFTS APB7F401.147
*CALL GCCOM GPB2F402.258
*ENDIF APB7F401.148
SETFIL1A.87
C*L DEFINE ARRAYS AND VARIABLES USED IN THIS ROUTINE----------------- SETFIL1A.88
C DEFINE LOCAL ARRAYS: NONE ARE REQUIRED SETFIL1A.89
SETFIL1A.90
C*--------------------------------------------------------------------- SETFIL1A.91
C DEFINE LOCAL VARIABLES SETFIL1A.92
SETFIL1A.93
REAL SETFIL1A.94
* COURANT !HOLDS MAXIMUM COURANT NUMBER. SETFIL1A.95
*, ADVECTION_COURANT ! HOLDS ADVECTION COURANT NUMBER. SETFIL1A.96
*, ADJUSTMENT_COURANT ! HOLDS ADJUSTMENT COURANT NUMBER. SETFIL1A.97
*, C_TERM ! HOLDS TERM WHICH MULTIPLIES EPSI_J IN SETFIL1A.103
* ! EQUATION 50. SETFIL1A.104
*, C_CONST ! HOLDS CONSTANT PART OF C_TERM SETFIL1A.105
SETFIL1A.106
REAL SETFIL1A.107
* SCALAR1 SETFIL1A.108
*, SCALAR2 SETFIL1A.109
*, ALPHA ! HOLDS DELTA LAMDA / COS(LAT) SETFIL1A.110
REAL APB7F401.149
& U_MAX_INITIAL(P_FIELD/ROW_LENGTH) APB7F401.150
! Holds maximum windspeed along each U row APB7F401.151
&, U_MAX(P_FIELD/ROW_LENGTH) ! Holds maximum windspeed APB7F401.152
! along each P row APB7F401.153
C TWO DIFFERENT NAMES NEEDED WHERE ONE WOULD SUFFICE TO STRICTLY CONFORM SETFIL1A.115
C WITH CRAY'S CASE CONSTRUCT. SEE ALSO CUT_OFF_N AND CUT_OFF_S. SETFIL1A.116
SETFIL1A.117
C COUNT VARIABLES FOR DO LOOPS ETC. SETFIL1A.118
INTEGER SETFIL1A.119
* I,K SETFIL1A.120
&, CUT_OFF ! HOLDS SMALLEST WAVE NUMBER FILTERED. APB7F401.154
*, P_ROWS ! HOLDS NUMBER OF ROWS OF DATA ON P_GRID SETFIL1A.123
*, HALF_P_ROWS ! HOLDS HALF NUMBER OF ROWS OF DATA ON P_GRID SETFIL1A.124
*, ROW ! LOOP COUNTER WHICH IS THE ROW ON WHICH SETFIL1A.131
* ! CALCULATION IS TAKING PLACE. SETFIL1A.132
*IF DEF,MPP APB7F401.155
INTEGER APB7F401.156
& ROW_START,ROW_END ! Start and end points for loop over rows APB7F401.157
&, global_row ! global row number of local row APB7F401.158
&, local_row ! local row number of global row APB7F401.159
&, glob_filterable_rows ! Global number of rows for filtering APB7F401.160
&, approx_ffts_per_proc ! Rough estimate of number of ffts/proc APB7F401.161
&, remainder_ffts ! "Spare" ffts to be added on APB7F401.162
&, filt_row_number ! filterable row number APB7F401.163
&, filt_row_index ! index for a (row,level) to be fft'd APB7F401.164
&, change_pt ! used in calculation of decomposition APB7F401.165
&, proc_row ! Current processor row APB7F401.166
&, proc_dest ! Processor to send a row to APB7F401.167
&, n_send ! No. of fft data blocks I will send APB7F401.168
&, n_recv ! No. of fft data blocks I will receive APB7F401.169
&, fld_type ! indicator of field type (P or U) APB7F401.170
&, north_filt_row ! local vn for fld_type APB7F401.171
&, south_filt_row ! local vn for fld_type APB7F401.172
&, n_glob_rows ! number of global rows for fld_type APB7F401.173
&, oproc_dest ! destination of last previous row APB7F401.174
&, sourcepe ! processor this row comes from APB7F401.175
&, row_number ! row number of this row APB7F401.176
&, fft_row_len ! row length of fft data (glsize(1)+2) GPB2F404.278
&, J ! loop counter APB7F401.177
&, info ! return code for comms APB7F401.178
APB7F401.179
APB7F401.182
LOGICAL APB7F401.183
& sending ! Am I sending this row? APB7F401.184
APB7F401.185
LOGICAL called_before APB7F401.186
DATA called_before/.FALSE./ APB7F401.187
INTEGER old_northern_filt,old_southern_filt APB7F401.188
SAVE called_before,old_northern_filt,old_southern_filt APB7F401.189
APB7F401.190
*ENDIF APB7F401.191
SETFIL1A.133
SETFIL1A.134
C*L EXTERNAL SUBROUTINE CALLS:------------------------------------ SETFIL1A.135
*IF DEF,CRAY SETFIL1A.136
*IF -DEF,MPP APB7F401.192
INTEGER index APB7F401.193
*ENDIF APB7F401.194
INTEGER ISAMAX SETFIL1A.137
EXTERNAL ISAMAX SETFIL1A.138
*ELSE SETFIL1A.141
C NO EXTERNAL SUBROUTINE CALLS SETFIL1A.142
*ENDIF SETFIL1A.143
C*--------------------------------------------------------------------- SETFIL1A.144
CL CALL COMDECK TO GET CONSTANTS USED. SETFIL1A.145
SETFIL1A.146
*CALL C_SETFIL SETFIL1A.147
SETFIL1A.148
CL MAXIMUM VECTOR LENGTH ASSUMED IS ROW_LENGTH SETFIL1A.149
CL--------------------------------------------------------------------- SETFIL1A.150
CL INTERNAL STRUCTURE. SETFIL1A.151
CL--------------------------------------------------------------------- SETFIL1A.152
CL SETFIL1A.153
APB7F401.195
*IF DEF,MPP APB7F401.196
APB7F401.197
IF (.NOT. called_before) THEN APB7F401.198
called_before=.TRUE. APB7F401.199
old_northern_filt=-1 APB7F401.200
old_southern_filt=-1 APB7F401.201
ENDIF APB7F401.202
APB7F401.203
*ENDIF APB7F401.204
APB7F401.205
C_CONST = C* ADJUSTMENT_TIMESTEP * LONGITUDE_STEP_INVERSE / A APB7F401.206
C_CONST = C_CONST * C_CONST APB7F401.207
ALPHA=1.0/LONGITUDE_STEP_INVERSE APB7F401.208
P_ROWS=P_FIELD/ROW_LENGTH APB7F401.209
*IF -DEF,MPP APB7F401.210
HALF_P_ROWS=P_ROWS/2 APB7F401.211
*ELSE APB7F401.212
HALF_P_ROWS=glsize(2)/2 ! half number of global rows APB7F401.213
*ENDIF APB7F401.214
APB7F401.215
! Caculate the local maximum U along each row APB7F401.216
APB7F401.217
*IF -DEF,MPP APB7F401.218
U_MAX_INITIAL(1)=0.0 APB7F401.219
U_MAX_INITIAL(P_ROWS)=0.0 APB7F401.220
DO ROW=2,P_ROWS-1 APB7F401.221
*ELSE APB7F401.222
DO ROW=1,Offy APB7F401.223
U_MAX_INITIAL(ROW)=0.0 APB7F401.224
U_MAX_INITIAL(P_ROWS+1-ROW)=0.0 APB7F401.225
ENDDO APB7F401.226
DO ROW=1+Offy,P_ROWS-Offy ! loop over local rows, missing halos APB7F401.227
*ENDIF APB7F401.228
SCALAR1=0.0 APB7F401.229
DO K=1,P_LEVELS APB7F401.230
*IF DEF,CRAY,AND,-DEF,MPP APB7F401.231
index=ISAMAX(ROW_LENGTH,U((ROW-1)*ROW_LENGTH+1,K),1) APB7F401.232
SCALAR1= APB7F401.233
& MAX(ABS(U((ROW-1)*ROW_LENGTH+index,K)),SCALAR1) APB7F401.234
*ELSE APB7F401.235
*IF -DEF,MPP APB7F401.236
DO I=(ROW-1)*ROW_LENGTH+1,ROW*ROW_LENGTH APB7F401.237
*ELSE APB7F401.238
DO I=(ROW-1)*ROW_LENGTH+1+Offx,ROW*ROW_LENGTH-Offx APB7F401.239
*ENDIF APB7F401.240
SCALAR1=MAX(ABS(U(I,K)),SCALAR1) APB7F401.241
ENDDO APB7F401.242
*ENDIF APB7F401.243
ENDDO ! loop over levels APB7F401.244
U_MAX_INITIAL(ROW)=SCALAR1 APB7F401.245
ENDDO ! loop over rows APB7F401.246
APB7F401.247
*IF DEF,MPP APB7F401.248
! We only have the local maximum so far. We must now do a maximum APB7F401.249
! along each global row. APB7F401.250
APB7F401.251
CALL GCG_RMAX(P_ROWS-2*Offy,gc_proc_row_group,info, APB7F401.252
& U_MAX_INITIAL(1+Offy)) APB7F401.253
APB7F401.254
! We need to know the values of the rows above and below us, so APB7F401.255
! fill in the halos with this information APB7F401.256
CALL SWAPBOUNDS(U_MAX_INITIAL,1,P_ROWS,0,Offy,1) APB7F401.257
*ENDIF APB7F401.258
APB7F401.259
! Now calculate the maximum "wind" along each P_ROW, by taking the APB7F401.260
! maximum wind of the U_ROW either side of it: APB7F401.261
! MAX_WIND(P_ROW)=MAX(U_MAX(P_ROW),U_MAX(P_ROW-1)) APB7F401.262
APB7F401.263
*IF -DEF,MPP APB7F401.264
DO ROW=2,P_ROWS-1 APB7F401.265
*ELSE APB7F401.266
IF (attop) THEN APB7F401.267
ROW_START=Offy+2 APB7F401.268
ELSE APB7F401.269
ROW_START=Offy+1 APB7F401.270
ENDIF APB7F401.271
IF (atbase) THEN APB7F401.272
ROW_END= P_ROWS-Offy-1 APB7F401.273
ELSE APB7F401.274
ROW_END=P_ROWS-Offy APB7F401.275
ENDIF APB7F401.276
APB7F401.277
DO ROW=ROW_START,ROW_END ! loop over local rows APB7F401.278
*ENDIF APB7F401.279
U_MAX(ROW)= APB7F401.280
& MAX(U_MAX_INITIAL(ROW),U_MAX_INITIAL(ROW-1)) APB7F401.281
! and scale U_MAX APB7F401.282
U_MAX(ROW)=U_MAX(ROW)*ADVECTION_TIMESTEP* APB7F401.283
& LONGITUDE_STEP_INVERSE*SEC_P_LATITUDE(ROW*ROW_LENGTH)/A APB7F401.284
& *(TWO_D_GRID_CORRECTION(ROW)+FILTERING_SAFETY_FACTOR) APB7F401.285
! and square it APB7F401.286
U_MAX(ROW)=U_MAX(ROW)*U_MAX(ROW) APB7F401.287
ENDDO APB7F401.288
APB7F401.289
! Now loop round each row, setting the filter_wave_number, and APB7F401.290
! calculating NORTHERN/SOUTHERN_FILTERED_P_ROW. APB7F401.291
*IF DEF,MPP APB7F401.292
! These values are an initial guess, and will have to be MAX/MINed in APB7F401.293
! the North/South direction in order to get the correct global values. APB7F401.294
*ENDIF APB7F401.295
APB7F401.296
NORTHERN_FILTERED_P_ROW=2 APB7F401.297
*IF -DEF,MPP APB7F401.298
SOUTHERN_FILTERED_P_ROW=P_ROWS APB7F401.299
CUT_OFF=ROW_LENGTH/2 APB7F401.300
APB7F401.301
DO ROW=2,P_ROWS-1 APB7F401.302
*ELSE APB7F401.303
SOUTHERN_FILTERED_P_ROW=glsize(2) APB7F401.304
CUT_OFF=glsize(1)/2 APB7F401.305
APB7F401.306
DO ROW=ROW_START,ROW_END ! loop over local rows APB7F401.307
global_row=ROW+datastart(2)-Offy-1 ! global row number APB7F401.308
*ENDIF APB7F401.309
C_TERM=SEC_P_LATITUDE(ROW*ROW_LENGTH)* APB7F401.310
& SEC_P_LATITUDE(ROW*ROW_LENGTH)*C_CONST* APB7F401.311
& TWO_D_GRID_CORRECTION(ROW)* APB7F401.312
& TWO_D_GRID_CORRECTION(ROW) APB7F401.313
APB7F401.314
DO I=2,CUT_OFF APB7F401.315
SCALAR1= SIN(ALPHA*(I+1)) ! EPSI_V APB7F401.316
SCALAR2 = SIN(I*ALPHA*.5) ! EPSI_J APB7F401.317
! First term in eqn 50 - Advection Courant Number APB7F401.318
ADVECTION_COURANT = U_MAX(ROW)*SCALAR1*SCALAR1 APB7F401.319
! Second term in eqn 50 - Adjustment Courant Number APB7F401.320
ADJUSTMENT_COURANT = C_TERM*SCALAR2*SCALAR2 APB7F401.321
COURANT = MAX(ADJUSTMENT_COURANT,ADVECTION_COURANT) APB7F401.322
IF (COURANT .GT. 1.0) THEN ! violated stability criteria APB7F401.323
APB7F401.324
! If Northern hemisphere: APB7F401.325
*IF -DEF,MPP APB7F401.326
IF (ROW .LE. HALF_P_ROWS+1) APB7F401.327
& NORTHERN_FILTERED_P_ROW=ROW APB7F401.328
*ELSE APB7F401.329
IF (global_row .LE. HALF_P_ROWS+1) APB7F401.330
& NORTHERN_FILTERED_P_ROW=global_row APB7F401.331
*ENDIF APB7F401.332
! If Southern hemisphere: APB7F401.333
*IF -DEF,MPP APB7F401.334
IF ((ROW .GE. HALF_P_ROWS+2) .AND. APB7F401.335
& (SOUTHERN_FILTERED_P_ROW .EQ. P_ROWS)) APB7F401.336
& SOUTHERN_FILTERED_P_ROW=ROW APB7F401.337
APB7F401.338
FILTER_WAVE_NUMBER_P_ROWS(ROW) = I-1 APB7F401.339
*ELSE APB7F401.340
IF ((global_row .GE. HALF_P_ROWS+2) .AND. APB7F401.341
& (SOUTHERN_FILTERED_P_ROW .EQ. glsize(2))) APB7F401.342
& SOUTHERN_FILTERED_P_ROW=global_row APB7F401.343
APB7F401.344
FILTER_WAVE_NUMBER_P_ROWS(global_row) = I-1 APB0F402.107
*ENDIF APB7F401.346
APB7F401.347
GOTO 100 APB7F401.348
ENDIF ! if stability criteria violated APB7F401.349
ENDDO ! loop over wave-numbers APB7F401.350
! If stability criteria not violated, set FILTER_WAVE_NUMBER APB7F401.351
! to CUT_OFF APB7F401.352
*IF -DEF,MPP APB7F401.353
FILTER_WAVE_NUMBER_P_ROWS(ROW) = CUT_OFF APB7F401.354
*ELSE APB7F401.355
FILTER_WAVE_NUMBER_P_ROWS(global_row) = CUT_OFF APB0F402.108
APB7F401.357
*ENDIF APB7F401.358
100 CONTINUE APB7F401.359
ENDDO ! loop over local rows APB7F401.360
APB7F401.361
*IF DEF,MPP APB7F401.362
! Work out the correct global values for FILTERED_P_ROWs APB7F401.363
APB7F401.364
CALL GC_IMAX(1,nproc,info,NORTHERN_FILTERED_P_ROW) APB7F401.365
CALL GC_IMIN(1,nproc,info,SOUTHERN_FILTERED_P_ROW) APB7F401.366
APB7F401.367
! Keep a copy of SOUTHERN_FILTERED_P_ROW - we'll use this to decide APB7F401.368
! if the field we're filtering is a U_FIELD or P_FIELD APB7F401.369
south_filt_p_row=SOUTHERN_FILTERED_P_ROW APB7F401.370
APB7F401.371
*ENDIF APB7F401.372
! We have to enforce the condition that the FILTER_WAVE_NUMBER_P_ROWS APB7F401.373
! cannot increase as we move towards the pole. APB7F401.374
*IF DEF,MPP APB7F401.375
! This is easiest if all the data is available on each processor APB7F401.376
APB7F401.377
! First, broadcast my local part of the filt_wave_number array APB7F401.378
APB7F401.379
DO I=0,nproc-1,nproc_x APB7F401.380
CALL GC_IBCAST(I,g_blsizep(2,I),I,nproc,info, APB7F401.381
& FILTER_WAVE_NUMBER_P_ROWS(g_datastart(2,I))) APB0F402.109
ENDDO APB7F401.383
*ENDIF APB7F401.384
APB7F401.385
! Northern hemisphere first: APB7F401.386
DO ROW=HALF_P_ROWS,2,-1 APB7F401.387
FILTER_WAVE_NUMBER_P_ROWS(ROW) = APB7F401.389
& MIN(FILTER_WAVE_NUMBER_P_ROWS(ROW), APB7F401.390
& FILTER_WAVE_NUMBER_P_ROWS(ROW+1)) APB7F401.391
ENDDO APB7F401.396
APB7F401.397
! And similarly for the Southern hemisphere APB7F401.398
*IF -DEF,MPP APB7F401.399
DO ROW=HALF_P_ROWS+2,P_ROWS-1 APB7F401.400
*ELSE APB0F402.110
DO ROW=HALF_P_ROWS+2,glsize(2)-1 APB0F402.111
*ENDIF APB0F402.112
FILTER_WAVE_NUMBER_P_ROWS(ROW) = APB7F401.401
& MIN(FILTER_WAVE_NUMBER_P_ROWS(ROW), APB7F401.402
& FILTER_WAVE_NUMBER_P_ROWS(ROW-1)) APB7F401.403
ENDDO APB7F401.404
APB7F401.411
! And set the values at the poles: APB7F401.412
*IF -DEF,MPP APB7F401.413
FILTER_WAVE_NUMBER_P_ROWS(1) = 0 APB7F401.414
FILTER_WAVE_NUMBER_P_ROWS(P_ROWS) = 0 APB7F401.415
*ELSE APB7F401.416
FILTER_WAVE_NUMBER_P_ROWS(1) = 0 APB0F402.113
FILTER_WAVE_NUMBER_P_ROWS(glsize(2)) = 0 APB0F402.114
*ENDIF APB7F401.419
APB7F401.420
! Make the U_FIELD version APB7F401.422
DO I=1,HALF_P_ROWS APB7F401.423
FILTER_WAVE_NUMBER_U_ROWS(I) = FILTER_WAVE_NUMBER_P_ROWS(I) APB7F401.424
*IF -DEF,MPP APB0F402.115
FILTER_WAVE_NUMBER_U_ROWS(P_ROWS-I) = APB7F401.425
& FILTER_WAVE_NUMBER_P_ROWS(P_ROWS+1-I) APB7F401.426
*ELSE APB0F402.116
FILTER_WAVE_NUMBER_U_ROWS(glsize(2)-I) = APB0F402.117
& FILTER_WAVE_NUMBER_P_ROWS(glsize(2)+1-I) APB0F402.118
*ENDIF APB0F402.119
ENDDO APB7F401.427
GPB3F403.104
IF ((NORTHERN_FILTERED_P_ROW .EQ. HALF_P_ROWS+1) .AND. GPB3F403.105
& (SOUTHERN_FILTERED_P_ROW .EQ. HALF_P_ROWS+2)) THEN GPB3F403.106
WRITE(6,*) '*** SETFIL : WARNING *** ', GPB3F403.107
& 'Filtering entire model domain.' GPB3F403.108
WRITE(6,*) 'Timestep should be reduced.' GPB0F405.44
*IF DEF,MPP GPB0F405.45
WRITE(6,*) 'Filtering has been switched off.' GPB0F405.46
*IF DEF,T3E GPB0F405.47
WRITE(0,*) '*** SETFIL : WARNING *** Filtering turned off' GPB0F405.48
*ENDIF GPB0F405.49
filter_off=.TRUE. GPB0F405.50
GOTO 9999 GPB0F405.51
*ENDIF GPB0F405.52
ENDIF GPB3F403.109
GPB3F403.110
*IF DEF,MPP APB0F402.120
! -------------------------------------------------------------------- APB7F401.458
APB7F401.459
IF ((old_northern_filt .NE. NORTHERN_FILTERED_P_ROW) .OR. APB7F401.460
& (old_southern_filt .NE. SOUTHERN_FILTERED_P_ROW)) THEN APB7F401.461
old_northern_filt = NORTHERN_FILTERED_P_ROW APB7F401.462
old_southern_filt = SOUTHERN_FILTERED_P_ROW APB7F401.463
APB7F401.464
GPB3F403.111
! These diagnostic write statments allow the user to see how GPB3F403.112
! much of the globe is being filtered. GPB3F403.113
GPB3F403.114
! WRITE(6,*) 'SETFIL : Updating filtering area to ', GPB3F403.115
! & '2 - > ',NORTHERN_FILTERED_P_ROW,' and ', GPB3F403.116
! & SOUTHERN_FILTERED_P_ROW,' -> ',glsize(2)-1 GPB3F403.117
! WRITE(6,*) 'SETFIL: Filtering ', GPB3F403.118
! & REAL((glsize(2)-2)- GPB3F403.119
! & (SOUTHERN_FILTERED_P_ROW-NORTHERN_FILTERED_P_ROW-1))/ GPB3F403.120
! & REAL(glsize(2)-2),' % of grid (not counting polar rows)' GPB3F403.121
fft_row_len=glsize(1)+2 GPB2F404.279
DO fld_type=1,2 ! 1=P_FIELD and 2=U_FIELD APB7F401.465
APB7F401.466
north_filt_row=NORTHERN_FILTERED_P_ROW APB7F401.467
IF (fld_type .EQ. 1) THEN APB7F401.468
south_filt_row=SOUTHERN_FILTERED_P_ROW APB7F401.469
n_glob_rows=glsize(2) APB7F401.470
ELSE APB7F401.471
south_filt_row=SOUTHERN_FILTERED_P_ROW-1 APB7F401.472
n_glob_rows=glsize(2)-1 APB7F401.473
ENDIF APB7F401.474
APB7F401.475
! Set up all the information required for decomposing the ffts APB7F401.476
! The information is passed to the filter routine via the APB7F401.477
! PAR_FFT common block APB7F401.478
APB7F401.479
glob_filterable_rows=n_glob_rows-(south_filt_row- APB7F401.480
& north_filt_row-1)-2 APB7F401.481
! The "-2" arises because we don't want to filter the polar APB7F401.482
! rows. APB7F401.483
APB7F401.484
approx_ffts_per_proc=(glob_filterable_rows*P_LEVELS)/nproc GPB2F404.280
remainder_ffts=(glob_filterable_rows*P_LEVELS)- GPB2F404.281
& approx_ffts_per_proc*nproc APB7F401.487
APB7F401.488
filter_off=.FALSE. APB7F401.489
IF (REAL(glob_filterable_rows*P_LEVELS)/REAL(nproc) .GT. GPB2F404.282
& MAX_ROWS_TO_FILTER) THEN APB7F401.491
WRITE(6,*) 'Error : Too many rows to filter' APB7F401.492
WRITE(6,*) 'Increase MAX_ROWS_TO_FILTER in comdeck ', APB7F401.493
& 'PARFFTS.' APB7F401.494
filter_off=.TRUE. APB7F401.495
GOTO 9999 APB7F401.496
ENDIF APB7F401.497
! approx_ffts_per_proc is the number of rows each processor would have APB7F401.498
! to fft, if the total number of fft rows divided exactly between APB7F401.499
! the total number of processors. APB7F401.500
! remainder_ffts is the number of fft rows left unassigned, if each APB7F401.501
! processor takes on approx_ffts_per_proc rows to do. APB7F401.502
APB7F401.503
! The rows are distributed like this: APB7F401.504
! The first "remainder_ffts" processors are given APB7F401.505
! approx_ffts_per_proc+1 rows each to fft APB7F401.506
! All the other processors are given approx_ffts_per_proc rows each. APB7F401.507
! The indexing of the filterable rows is like this: APB7F401.508
! index global_row_number level APB7F401.509
! 1 2 1 APB7F401.510
! 2 2 2 APB7F401.511
! : 2 : APB7F401.512
! : 2 maxlevel APB7F401.513
! : 3 1 APB7F401.514
! etc. APB7F401.515
APB7F401.516
n_send=0 ! number of rows I will send APB7F401.517
n_recv=0 ! number of processor rows I will receive from APB7F401.518
fft_rows(fld_type)=0 ! total number of rows I will fft APB7F401.519
proc_dest = -1 APB7F401.522
change_pt=(approx_ffts_per_proc+1)*remainder_ffts APB7F401.523
APB7F401.524
APB7F401.525
DO ROW=2,n_glob_rows-1 ! loop over global rows, missing poles APB7F401.526
APB7F401.527
IF ((ROW .LE. north_filt_row) .OR. APB7F401.528
& (ROW .GE. south_filt_row)) THEN APB7F401.529
! This row needs to be filtered APB7F401.530
APB7F401.531
IF (ROW .LE. north_filt_row) THEN APB7F401.532
! Northern hemisphere APB7F401.533
filt_row_number=ROW-1 APB7F401.534
ELSE APB7F401.535
! Southern hemisphere APB7F401.536
filt_row_number=ROW- APB7F401.537
& (south_filt_row- APB7F401.538
& north_filt_row-1)-1 APB7F401.539
ENDIF APB7F401.540
APB7F401.541
! Is this row of data on my processor? APB7F401.542
sending=.FALSE. APB7F401.543
IF ((ROW .GE. datastart(2)) .AND. APB7F401.544
& (ROW .LE. datastart(2)+blsizep(2)-1)) sending=.TRUE. APB7F401.545
APB7F401.546
! Which processor row does it live on? APB7F401.547
DO I=0,nproc_y-1 ! loop over processor row APB7F401.548
local_row=ROW-g_datastart(2,I*nproc_x)+1 APB7F401.549
IF ((local_row .LE. g_blsizep(2,I*nproc_x)) .AND. APB7F401.550
& (local_row .GT. 0)) proc_row=I APB7F401.551
ENDDO APB7F401.552
APB7F401.553
DO K=1,P_LEVELS GPB2F404.283
filt_row_index=(K+(filt_row_number-1)*P_LEVELS) GPB2F404.284
APB7F401.556
! Which processor is it going to? APB7F401.557
oproc_dest = proc_dest APB7F401.558
IF (filt_row_index .LE. change_pt) THEN APB7F401.559
! It's a processor with approx_ffts_per_proc+1 rows to do APB7F401.560
proc_dest=(filt_row_index-1)/(approx_ffts_per_proc+1) APB7F401.561
ELSE APB7F401.562
! It's a processor with approx_ffts_per_proc rows to do APB7F401.563
proc_dest=remainder_ffts+ APB7F401.564
& (filt_row_index-1-change_pt)/approx_ffts_per_proc APB7F401.565
ENDIF APB7F401.566
APB7F401.567
! Am I sending it? APB7F401.568
IF (sending) THEN APB7F401.569
if (proc_dest .eq. oproc_dest .and. k .ge. 2) then APB7F401.570
filt_send_max(n_send,fld_type) = APB7F401.571
& filt_send_max(n_send,fld_type) + 1 APB7F401.572
else APB7F401.573
n_send = n_send + 1 APB7F401.574
filt_send_start(n_send,fld_type) = K APB7F401.575
filt_send_max(n_send,fld_type) = 1 APB7F401.576
filt_send_map(S_DESTINATION_PE,n_send,fld_type) = GPB2F402.259
& proc_dest GPB2F402.260
GPB2F402.261
filt_send_map(S_BASE_ADDRESS_IN_SEND_ARRAY, GPB2F402.262
& n_send,fld_type) = GPB2F402.263
& (K-1)*p_field + GPB2F402.264
& (ROW-datastart(2)+offy)*lasize(1) + offx + 1 GPB2F402.265
GPB2F402.266
filt_send_map(S_STRIDE_IN_SEND_ARRAY, GPB2F402.267
& n_send,fld_type) = GPB2F402.268
& p_field GPB2F402.269
GPB2F402.270
filt_send_map(S_ELEMENT_LENGTH,n_send,fld_type) = GPB2F402.271
& blsizep(1) GPB2F402.272
! Calculate the number of this fft within the receiving processor APB7F401.582
! We do this by calculating the first fft the receiving processor does APB7F401.583
! and taking this number away from the filt_row_index APB7F401.584
IF (proc_dest .LE. remainder_ffts-1) THEN APB7F401.585
GPB2F402.273
filt_send_map(S_BASE_ADDRESS_IN_RECV_ARRAY, GPB2F402.274
& n_send,fld_type)= GPB2F402.275
& (filt_row_index GPB2F402.276
& +1-(proc_dest*(approx_ffts_per_proc+1)+1)-1) APB7F401.587
& *fft_row_len + datastart(1) GPB2F404.285
ELSE APB7F401.589
GPB2F402.277
filt_send_map(S_BASE_ADDRESS_IN_RECV_ARRAY, GPB2F402.278
& n_send,fld_type)= GPB2F402.279
& (filt_row_index+ GPB2F402.280
& 1-((remainder_ffts*(approx_ffts_per_proc+1)) APB7F401.591
& +(proc_dest-remainder_ffts) APB7F401.592
& *approx_ffts_per_proc+1) APB7F401.593
& -1)*fft_row_len + datastart(1) GPB2F404.286
ENDIF APB7F401.595
GPB2F402.281
filt_send_map(S_STRIDE_IN_RECV_ARRAY, GPB2F402.282
& n_send,fld_type) = GPB2F402.283
& fft_row_len GPB2F404.288
GPB2F402.285
endif APB7F401.597
ENDIF ! IF sending APB7F401.598
APB7F401.599
! Am I receiving it? APB7F401.600
IF (proc_dest .EQ. mype) THEN ! Is this my row to fft? APB7F401.601
fft_rows(fld_type)=fft_rows(fld_type)+1 APB7F401.602
filt_info(fft_rows(fld_type),fld_type) = row APB7F401.603
filt_level(fft_rows(fld_type),fld_type) = K GPB2F404.287
APB7F401.604
if (proc_dest .eq. oproc_dest .and. k .ge. 2) then APB7F401.605
do j = 0,nproc_x-1 APB7F401.606
filt_recv_max(n_recv-j,fld_type) = APB7F401.607
& filt_recv_max(n_recv-j,fld_type) + 1 APB7F401.608
enddo APB7F401.609
else APB7F401.610
DO J=0,nproc_x-1 ! loop over processor row APB7F401.611
APB7F401.612
sourcepe=proc_row*nproc_x+J APB7F401.613
row_number= APB7F401.614
& ROW-g_datastart(2,proc_row*nproc_x)+Offy+1 APB7F401.615
! The following loop could be replaced with a strided receive APB7F401.616
APB7F401.617
n_recv = n_recv + 1 APB7F401.618
filt_recv_start(n_recv,fld_type) = K APB7F401.619
filt_recv_max(n_recv,fld_type) = 1 APB7F401.620
filt_recv_map(R_SOURCE_PE,n_recv,fld_type) = GPB2F402.286
& sourcepe GPB2F402.287
GPB2F402.288
filt_recv_map(R_BASE_ADDRESS_IN_RECV_ARRAY, GPB2F402.289
& n_recv,fld_type) = GPB2F402.290
& (fft_rows(fld_type)-1)*fft_row_len + GPB2F404.289
& g_datastart(1,sourcepe) GPB2F402.292
GPB2F402.293
filt_recv_map(R_STRIDE_IN_RECV_ARRAY, GPB2F402.294
& n_recv,fld_type) = GPB2F402.295
& fft_row_len GPB2F404.290
GPB2F402.297
filt_recv_map(R_ELEMENT_LENGTH, GPB2F402.298
& n_recv,fld_type) = GPB2F402.299
& g_blsizep(1,sourcepe) GPB2F402.300
GPB2F402.301
filt_recv_map(R_BASE_ADDRESS_IN_SEND_ARRAY, GPB2F402.302
& n_recv,fld_type) = GPB2F402.303
& (K-1)*g_lasize(1,sourcepe)*g_lasize(2,sourcepe) GPB2F402.304
& + 1 + offx + (row_number-1)* GPB2F402.305
& g_lasize(1,sourcepe) GPB2F402.306
GPB2F402.307
filt_recv_map(R_STRIDE_IN_SEND_ARRAY, GPB2F402.308
& n_recv,fld_type) = GPB2F402.309
& g_lasize(1,sourcepe)*g_lasize(2,sourcepe) GPB2F402.310
APB7F401.634
ENDDO ! J APB7F401.635
endif APB7F401.636
APB7F401.637
ENDIF ! IF my row to fft APB7F401.638
APB7F401.639
ENDDO ! loop over levels (K) APB7F401.640
ENDIF ! IF this row needs to be filtered APB7F401.641
ENDDO ! loop over rows (ROW) APB7F401.642
APB7F401.643
n_items_to_send(fld_type) = n_send GJC0F405.36
n_items_to_recv(fld_type) = n_recv GJC0F405.37
APB7F401.646
ENDDO ! loop over field types (fld_type) APB7F401.647
ENDIF ! IF we need to update fft decomposition APB7F401.648
APB7F401.649
! APB7F401.650
! WRITE(6,*) 'Processor ',mype APB7F401.651
! WRITE(6,*) 'Sending ',n_items_to_send,' chunks.' APB7F401.652
! WRITE(6,*) 'Receiving ',n_items_to_recv,' chunks.' APB7F401.653
! WRITE(6,*) 'Doing ',fft_rows,' ffts.' APB7F401.654
! APB7F401.655
! write(6,*) 'filter wave nos:' APB7F401.656
! do i=1,glsize(2) APB7F401.657
! write(6,*) i,filt_wave_number(i) APB7F401.658
! enddo APB7F401.659
APB7F401.660
9999 CONTINUE APB7F401.661
APB7F401.662
*ENDIF APB7F401.663
APB7F401.664
! END OF ROUTINE SET_FIL APB7F401.665
APB7F401.666
RETURN APB7F401.667
END APB7F401.668
*ENDIF SETFIL1A.496