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