*IF DEF,C84_1A,OR,DEF,C98_1A,OR,DEF,FLDIO,OR,DEF,UTILIO UIE3F404.5
C ******************************COPYRIGHT****************************** GTS2F400.1027
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.1028
C GTS2F400.1029
C Use, duplication or disclosure of this code is subject to the GTS2F400.1030
C restrictions as set forth in the contract. GTS2F400.1031
C GTS2F400.1032
C Meteorological Office GTS2F400.1033
C London Road GTS2F400.1034
C BRACKNELL GTS2F400.1035
C Berkshire UK GTS2F400.1036
C RG12 2SZ GTS2F400.1037
C GTS2F400.1038
C If no contract has been raised with this copy of the code, the use, GTS2F400.1039
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.1040
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.1041
C Modelling at the above address. GTS2F400.1042
C ******************************COPYRIGHT****************************** GTS2F400.1043
C GTS2F400.1044
CLL SUBROUTINE COEX,COEX2,CMPS,XPND,INSTIN,EXTRIN ----------------- COEX1A.3
CLL COEX1A.4
CLL PURPOSE: PACK TO AND UNPACK FROM WGDOS FORMAT COEX1A.5
CLL COEX1A.6
CLL (Note that for optimal performance the routines INSTIN and COEX1A.7
CLL EXTRIN are inline expanded (enable option 8 on fpp) COEX1A.8
CLL COEX1A.9
CLL Model Modification history from model version 3.0: COEX1A.10
CLL version Date COEX1A.11
CLL 3.2 19/04/93 Code for ne real missing data indicator (TCJ). TJ050593.77
CLL 4.2 Nov. 96 T3E migration: WHENNE,WHENEQ replaced by GSS9F402.3
CLL portable fortran code. S.J.Swarbrick GSS9F402.4
CLL 4.3 Apr. 97 T3E migration: Calling the functions CRI2IBM/IBM2CRI UIE2F403.1
CLL directly replaces the calls to the (now unsupported) UIE2F403.2
CLL routines USICTI, USICTC, USSCTC, USSCTI. Ian Edmond UIE2F403.3
CLL COEX1A.12
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4, COEX1A.13
CLL STANDARD B, VERSION 2, DATED 18/01/90 COEX1A.14
CLL COEX1A.15
CLL Logical component number: S72 COEX1A.16
CLL COEX1A.17
CLL SYSTEM TASK: P7 COEX1A.18
CLL OWNED BY P J SMITH COEX1A.19
CLL COEX1A.20
CLL DOCUMENTATION: ?????? COEX1A.21
CLLEND------------------------------------------------------------- COEX1A.22
COEX1A.23
SUBROUTINE COEX(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI,LWORD) 15,2TJ050593.78
INTEGER N,ICOMP(N),M,IX,IY,NUM,ISC,LWORD COEX1A.25
REAL FIELD(M),RMDI TJ050593.79
LOGICAL OCO COEX1A.27
INTEGER JJ,IFAIL COEX1A.28
C INTEGER IRC(IY) ! this was cause of errors as IY not set COEX1A.29
EXTERNAL COEX2 UIE2F403.4
C COEX1A.31
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.32
C COEX1A.33
C COEX1A.34
C OCO=.TRUE. OCO=.FALSE. COEX1A.35
C COEX1A.36
C FIELD = FIELD FOR COMPRESSING RETURNED EXPANDED DATA COEX1A.37
C M = SIZE OF FIELD SIZE OF FIELD COEX1A.38
C ICOMP = RETURNED COMPRESSED DATA FIELD FOR EXPANSION COEX1A.39
C N = SIZE OF COMP ------- COEX1A.40
C IX = X DIMENSION OF FIELD RETURNED X DIMENSION COEX1A.41
C IY = Y DIMENSION OF FIELD RETURNED Y DIMENSION COEX1A.42
C NUM = TOTAL NO. OF COMPRESSED ------- COEX1A.43
C (32 BIT) WORDS RETURNED COEX1A.44
C ISC = ACCURACY IN POWER OF 2 COEX1A.45
C OCO = .TRUE. FOR COMPRESSION .FALSE. FOR EXPANSION COEX1A.46
C COEX1A.47
C USERS OF THIS ROUTINE SHOULD ENSURE THAT THE ARRAY 'COMP' COEX1A.48
C WILL BE BIG ENOUGH TO HOLD THE COMPRESSED DATA. COEX1A.49
C IF INPUT ARRAY IS ONE DIMENSIONAL PUT IY=1 AND IX=DIMENSION, COEX1A.50
C WHEN USING WITH PRINTFILE FIELDS USE IX=192,IY=121 NOT IX=23232, COEX1A.51
C IY=1 THIS WILL MAKE THE COMPRESSION MORE EFFICIENT. COEX1A.52
C FOR MOST PRINTFILE FIELDS USING AN ACCURACY IF AN EIGTH (ISC=-3) COEX1A.53
C A COMPRESSION FACTOR OF 4 CAN BE ACHIEVED. FOR DDDFFF FIELDS COEX1A.54
C USERS ARE ADVISED TO SPLIT THE FIELD INTO U AND V COEX1A.55
C COMPONENTS. COEX1A.56
C COEX1A.57
C CRAY ROUTINE 64 BIT WORDS - CRAY FULLWORD COEX1A.58
C COEX1A.59
IF (OCO) THEN COEX1A.60
C WRITE(6,*) 'CRAY COEX - PACKING' GIE0F403.116
C COEX1A.62
C COMPRESSION OF DATA COEX1A.63
C COEX1A.64
C COEX1A.65
C INSERT SCALE FACTOR, COLS AND ROWS INTO HEADER COEX1A.66
C COEX1A.67
IER=CRI2IBM(2,1,ICOMP(1),32,ISC,1,64,32) UIE2F403.5
IER=CRI2IBM(2,1,ICOMP(2),0,IX,1,64,16) UIE2F403.6
IER=CRI2IBM(2,1,ICOMP(2),16,IY,1,64,16) UIE2F403.7
C COEX1A.71
C CALL AUTOTASKED/DYNAMIC ARRAY PART OF COEX COEX1A.72
C COEX1A.73
CALL COEX2
(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI) TJ050593.80
C CALL COEX2(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI,IRC) TJ050593.81
C COEX1A.76
C CHECK RETURN CODES COEX1A.77
C COEX1A.78
C IFAIL=0 COEX1A.79
C DO JJ = 1,IY COEX1A.80
C IF (IRC(JJ).NE.0) IFAIL=IFAIL+1 COEX1A.81
C ENDDO COEX1A.82
C IF(IFAIL.GT.0)THEN COEX1A.83
C DO JJ = 1,IY COEX1A.84
C IF (IRC(JJ).NE.0) THEN COEX1A.85
C WRITE(6,*)'RETURN CODE',IRC(JJ),'FROM COEX FOR ROW',JJ GIE0F403.117
C ENDIF COEX1A.87
C ENDDO COEX1A.88
C ENDIF COEX1A.89
C COEX1A.90
C END OF COMPRESSION SECTION COEX1A.91
C COEX1A.92
ELSE COEX1A.93
C WRITE(6,*) 'CRAY COEX - UNPACKING' GIE0F403.118
C COEX1A.95
C EXPANSION SECTION COEX1A.96
C COEX1A.97
C COEX1A.98
C EXTRACT SCALE FACTOR, COLS AND ROWS FROM HEADER COEX1A.99
C COEX1A.100
IER=IBM2CRI(2,1,ICOMP(1),32,ISC,1,64,32) UIE2F403.8
IER=IBM2CRI(2,1,ICOMP(2),0,IX,1,64,16) UIE2F403.9
IER=IBM2CRI(2,1,ICOMP(2),16,IY,1,64,16) UIE2F403.10
C COEX1A.104
C CALL AUTOTASKED/DYNAMIC ARRAY PART OF COEX COEX1A.105
C COEX1A.106
CALL COEX2
(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI) TJ050593.82
C CALL COEX2(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI,IRC) TJ050593.83
C COEX1A.109
C CHECK RETURN CODES (NOT YET IMPLEMENTED) COEX1A.110
C COEX1A.111
C DO JJ = 1,IY COEX1A.112
C IF (IRC(JJ).NE.0) THEN COEX1A.113
C WRITE(6,*)' NON-ZERO RETURN CODE FOR ROW ',JJ,IRC(JJ) GIE0F403.119
C ENDIF COEX1A.115
C ENDDO COEX1A.116
END IF COEX1A.117
C COEX1A.118
C COEX1A.119
RETURN COEX1A.120
COEX1A.121
END COEX1A.122
SUBROUTINE COEX2(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI) 2,2TJ050593.84
C SUBROUTINE COEX2(FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI,IRC) TJ050593.85
INTEGER N,ICOMP(N),M,IX,IY,NUM,ISC COEX1A.125
C INTEGER IRC(IY) COEX1A.126
REAL FIELD(M),RMDI TJ050593.86
LOGICAL OCO COEX1A.128
INTEGER JJ,KK,ILEN,IST,ICX,JCX,JBIT,JLEN,NOB,NOC,IER2 COEX1A.129
INTEGER IC(IY),ICB(IY),NOP(IY),JC(IY),IBIT(IY),ISTART(IY) COEX1A.131
INTEGER JCOMP(IX,IY),IERR1(IY),IERR2(IY),IERR3(IY) COEX1A.132
REAL ACC,APREC COEX1A.133
REAL BASE(IY) COEX1A.134
EXTERNAL CMPS,XPND,STRMOV UIE2F403.11
C COEX1A.136
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.137
C CRAY ROUTINE 64 BIT WORDS - CRAY FULLWORD COEX1A.138
C COEX1A.139
C OCO=.TRUE. OCO=.FALSE. COEX1A.140
C COEX1A.141
C FIELD = FIELD FOR COMPRESSING RETURNED EXPANDED DATA COEX1A.142
C M = SIZE OF FIELD SIZE OF FIELD COEX1A.143
C ICOMP = RETURNED COMPRESSED DATA FIELD FOR EXPANSION COEX1A.144
C N = SIZE OF COMP ------- COEX1A.145
C IX = X DIMENSION OF FIELD X DIMENSION OF FIELD COEX1A.146
C IY = Y DIMENSION OF FIELD Y DIMENSION OF FIELD COEX1A.147
C NUM = TOTAL NO. OF COMPRESSED ------- COEX1A.148
C (32 BIT) WORDS RETURNED COEX1A.149
C ISC = ACCURACY IN POWER OF 2 ACCURACY IN POWER OF 2 COEX1A.150
C OCO = .TRUE. FOR COMPRESSION .FALSE. FOR EXPANSION COEX1A.151
C IRC = RETURN CODE FOR EACH ROW RETURN CODE FOR EACH ROW COEX1A.152
C COEX1A.153
C INITIALISE TEMP. VARIABLES/ARRAYS COEX1A.154
C COEX1A.155
IF (.NOT. OCO) THEN GBCQF405.457
DO JJ=1,IY GBCQF405.458
NOP(JJ)=0 GBCQF405.459
IBIT(JJ)=0 GBCQF405.460
ENDDO GBCQF405.461
DO JJ=1,IY GBCQF405.462
BASE(JJ)=0.0 GBCQF405.463
IBIT(JJ)=0.0 GBCQF405.464
ENDDO GBCQF405.465
ENDIF GBCQF405.466
GBCQF405.467
DO JJ=1,IY COEX1A.173
DO JCX=1,IX COEX1A.174
JCOMP(JCX,JJ) = 0 COEX1A.175
END DO COEX1A.176
END DO COEX1A.177
IF (OCO) THEN COEX1A.178
C COEX1A.179
C COMPRESSION OF DATA COEX1A.180
C COEX1A.181
IC(1) = 1 COEX1A.182
ACC = ISC COEX1A.183
APREC = 2.**ACC COEX1A.184
C COEX1A.185
C PUT PACKED DATA FOR EACH ROW INTO TEMP ARRAY JCOMP COEX1A.186
C COEX1A.187
CFPP$ CNCALL COEX1A.188
DO 10 JJ = 1,IY COEX1A.189
C IP = POSITION IN INPUT ARRAY COEX1A.190
IP = (JJ-1)*IX + 1 COEX1A.191
CALL CMPS
(IX,FIELD(IP),JCOMP(2,JJ),NOP(JJ),APREC, COEX1A.192
+ IBIT(JJ),BASE(JJ),RMDI) TJ050593.87
10 CONTINUE COEX1A.194
C COEX1A.195
C ADD BASE VALUE, NUMBER OF BITS, AND LENGTH COEX1A.196
C COEX1A.197
CFPP$ CNCALL COEX1A.198
DO 11 JJ = 1,IY COEX1A.199
IERR1(JJ)=CRI2IBM(3,1, JCOMP(1,JJ),0, BASE(JJ),1,64,32) UIE2F403.12
IERR2(JJ)=CRI2IBM(2,1,JCOMP(1,JJ),32,IBIT(JJ),1,64,16) UIE2F403.13
IERR3(JJ)=CRI2IBM(2,1,JCOMP(1,JJ),48,NOP(JJ),1,64,16) UIE2F403.14
11 CONTINUE COEX1A.203
C COEX1A.204
C CHECK ROW HEADER AND SET RETURN CODES COEX1A.205
C COEX1A.206
CFPP$ NOCONCUR COEX1A.207
C DO 12 JJ = 1,IY COEX1A.208
C IF (IERR1(JJ).NE.0) IRC(JJ) = IRC(JJ) + 1 COEX1A.209
C IF (IERR2(JJ).NE.0) IRC(JJ) = IRC(JJ) + 2 COEX1A.210
C IF (IERR3(JJ).NE.0) IRC(JJ) = IRC(JJ) + 4 COEX1A.211
C IF (JCOMP(1,JJ).EQ.0) THEN COEX1A.212
C IF (IBIT(JJ).NE.0) IRC(JJ) = IRC(JJ) + 8 COEX1A.213
C IF (NOP(JJ) .NE.0) IRC(JJ) = IRC(JJ) + 16 COEX1A.214
C ENDIF COEX1A.215
C 12 CONTINUE COEX1A.216
C COEX1A.217
C CALCULATE POSITIONS IN OUTPUT ARRAY FOR PACKED ROWS COEX1A.218
C (FIRST PACKED ROW STARTS AT WORD 1; BIT 31) COEX1A.219
C COEX1A.220
IC(1) = 2 COEX1A.221
ICB(1) = -1 COEX1A.222
ISTART(1) = 5 COEX1A.223
CFPP$ NOCONCUR COEX1A.224
DO 20 JJ = 2,IY COEX1A.225
IF (MOD(NOP(JJ-1),2).EQ.1) THEN COEX1A.226
IC(JJ ) = IC(JJ-1) + NOP(JJ-1)/2 + 1 COEX1A.227
ICB(JJ) = -ICB(JJ-1) COEX1A.228
IF (ICB(JJ).GT.0) IC(JJ) = IC(JJ) + 1 COEX1A.229
ELSE COEX1A.230
IC(JJ) = IC(JJ-1) + (NOP(JJ-1)+1)/2 + 1 COEX1A.231
ICB(JJ) = ICB(JJ-1) COEX1A.232
ENDIF COEX1A.233
ISTART(JJ) = 5 COEX1A.234
IF(ICB(JJ).EQ.1) ISTART(JJ) = 1 COEX1A.235
20 CONTINUE COEX1A.236
C COEX1A.237
C MOVE TEMP. ARRAY INTO OUTPUT ARRAY COEX1A.238
C COEX1A.239
CFPP$ NOCONCUR COEX1A.240
DO 30 JJ=1,IY COEX1A.241
NOB = NOP(JJ)*4 + 8 COEX1A.242
C CHECK IF PACKED FIELD GREATER THAN UN PACKED FIELD COEX1A.243
C IF(NOB.GT.IX*8)IRC(JJ)=IRC(JJ)+32 COEX1A.244
IST = ISTART(JJ) COEX1A.245
ICX = IC(JJ) COEX1A.246
CALL STRMOV(
JCOMP(1,JJ),1,NOB,ICOMP(ICX),IST) COEX1A.247
30 CONTINUE COEX1A.248
C COEX1A.249
C INSERT TOTAL LENGTH OF THIS FIELD COEX1A.250
NUM = IC(IY)*2 + NOP(IY) COEX1A.251
IF (ICB(IY).LT.0) NUM = NUM + 1 COEX1A.252
IER2=CRI2IBM(2,1,ICOMP(1),0,NUM,1,64,32) UIE2F403.15
C COEX1A.254
C END OF COMPRESSION SECTION COEX1A.255
C COEX1A.256
ELSE COEX1A.257
C COEX1A.258
C EXPANSION SECTION COEX1A.259
C COEX1A.260
ACC = ISC COEX1A.261
APREC = 2.**ACC COEX1A.262
ICX = 2 COEX1A.263
JCX = -1 COEX1A.264
CFPP$ NOCONCUR COEX1A.265
DO 40 JJ = 1,IY COEX1A.266
C COEX1A.267
C MOVE PACKED ROWS INTO TEMP ARRAYS COEX1A.268
C COEX1A.269
IF (JCX.LT.0) THEN COEX1A.270
C COEX1A.271
C EXTRACT BASE, NO. BITS, NO 32 BIT WORDS COEX1A.272
C COEX1A.273
IERR=IBM2CRI(3,1,ICOMP(ICX),32,BASE(JJ),1,64,32) UIE2F403.16
IERR=IBM2CRI(2,1,ICOMP(ICX+1),0,IBIT(JJ),1,64,16) UIE2F403.17
IERR=IBM2CRI(2,1,ICOMP(ICX+1),16,NOP(JJ),1,64,16) UIE2F403.18
C SAVE START POSITION OF ROW COEX1A.277
IC(JJ) = ICX COEX1A.278
ISTART(JJ) = 5 COEX1A.279
ELSE COEX1A.280
IERR=IBM2CRI(3,1,ICOMP(ICX),0,BASE(JJ),1,64,32) UIE2F403.19
IERR=IBM2CRI(2,1,ICOMP(ICX),32,IBIT(JJ),1,64,16) UIE2F403.20
IERR=IBM2CRI(2,1,ICOMP(ICX),48,NOP(JJ),1,64,16) UIE2F403.21
C SAVE START POSITION OF ROW COEX1A.284
IC(JJ) = ICX COEX1A.285
ISTART(JJ) = 1 COEX1A.286
END IF COEX1A.287
C COEX1A.288
C CALCULATE START OF NEXT ROW COEX1A.289
C COEX1A.290
IF (MOD(NOP(JJ),2).EQ.1) THEN COEX1A.291
ICX = ICX + NOP(JJ)/2 + 1 COEX1A.292
JCX = -JCX COEX1A.293
IF (JCX.GT.0) ICX = ICX + 1 COEX1A.294
ELSE COEX1A.295
ICX = ICX + (NOP(JJ)+1)/2 + 1 COEX1A.296
END IF COEX1A.297
40 CONTINUE COEX1A.298
C COEX1A.299
C MOVE EACH PACKED ROW INTO TEMP ARRAY JCOMP COEX1A.300
C COEX1A.301
CFPP$ NOCONCUR COEX1A.302
DO 50 JJ = 1,IY COEX1A.303
ICX = IC(JJ) COEX1A.304
IST = ISTART(JJ) COEX1A.305
NOB = NOP(JJ)*4 + 8 COEX1A.306
CALL STRMOV(
ICOMP(ICX),IST,NOB,JCOMP(1,JJ),1) COEX1A.307
50 CONTINUE COEX1A.308
C COEX1A.309
C CALCULATE START OF EACH ROW IN FIELD COEX1A.310
C COEX1A.311
ICX = 1 COEX1A.312
DO 60 JJ = 1,IY COEX1A.313
IC(JJ) = ICX COEX1A.314
ICX = ICX + IX COEX1A.315
60 CONTINUE COEX1A.316
C COEX1A.317
C UNPACK DATA INTO FIELD COEX1A.318
C COEX1A.319
CFPP$ CNCALL COEX1A.320
DO 70 JJ=1,IY COEX1A.321
ICX = IC(JJ) COEX1A.322
CALL XPND
(IX,JCOMP(2,JJ),FIELD(ICX),APREC,IBIT(JJ), COEX1A.323
+ BASE(JJ),NOP(JJ),RMDI) TJ050593.88
70 CONTINUE COEX1A.325
END IF COEX1A.326
RETURN COEX1A.327
COEX1A.328
END COEX1A.329
SUBROUTINE CMPS(IX,FIELD,ICOMP,NOP,APREC,IBIT,BASE,RMDI) 2,1TJ050593.89
CFPP$ NOCONCUR R COEX1A.331
CFPP$ EXPAND (INSTIN) COEX1A.332
INTEGER IX,NOP,ICOMP(IX-1),IBIT COEX1A.333
REAL FIELD(IX),APREC,BASE,RMDI TJ050593.90
INTEGER I,IBASE,II,IMAX,I2,JBIT,JJ,JWORD,NMIS,NPOINT,NZERO COEX1A.335
INTEGER IMAP(IX),IMIS(IX),IZERO(IX),ITEMP(IX) GBCQF405.468
INTEGER IGATH1(IX),IGATH2(IX) COEX1A.337
REAL ATEMP(IX),BPREC COEX1A.338
LOGICAL OBTMIS,OBTZER COEX1A.339
C COEX1A.340
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.341
C COEX1A.342
C IX = LENTH OF FIELD COEX1A.343
C FIELD = FIELD FOR COMPRESSING COEX1A.344
C ICOMP = RETURNED COMPRESSED DATA COEX1A.345
C NOP = NUMBER OF WORDS OF COMP FILLED BY THIS CALL COEX1A.346
C APREC = PRECISION COEX1A.347
C IBIT = NUMBER OF BITS INTO WHICH DATA PACKED COEX1A.348
C BASE = REFERENCE (MIN) VALUE FOR PACKED DATA COEX1A.349
C RMDI = MISSING DATA INDICATOR TJ050593.91
C COEX1A.351
C INITIALISE VARIABLES/TEMP ARRAYS COEX1A.352
C COEX1A.353
OBTMIS = .FALSE. COEX1A.354
OBTZER = .FALSE. COEX1A.355
BPREC = 1./APREC COEX1A.356
DO I = 1,IX COEX1A.357
IMAP(I) = 1 COEX1A.358
END DO COEX1A.364
DO I = 1,IX-1 COEX1A.365
ICOMP(I) = 0 COEX1A.366
END DO COEX1A.367
C COEX1A.368
C SCAN FIELD FOR MISSING DATA AND STORE RESULT IN IMIS, COEX1A.369
C SCAN FIELD FOR ZERO VALUES AND STORE RESULT IN IZERO, COEX1A.370
C SCAN FIELD FOR MINIMUM VALUE (IGNORE MISSING DATA INDICATOR) COEX1A.371
C COEX1A.372
BASE = 999999. COEX1A.373
NMIS = 0 COEX1A.374
NZERO = 0 COEX1A.375
JJ = 0 COEX1A.376
COEX1A.377
DO I = 1,IX COEX1A.378
IF (FIELD(I).NE.RMDI) THEN GBCQF405.469
IMIS(I)=0 GBCQF405.470
ELSE GBCQF405.471
IMIS(I)=1 GBCQF405.472
ENDIF GBCQF405.473
END DO COEX1A.380
COEX1A.381
DO I = 1,IX COEX1A.382
IF (FIELD(I).NE.0.0) THEN GBCQF405.474
IZERO(I)=0 GBCQF405.475
ELSE GBCQF405.476
IZERO(I)=1 GBCQF405.477
ENDIF GBCQF405.478
END DO COEX1A.384
COEX1A.385
C COEX1A.386
C GET NO. OF NON-RMDI POINTS + COMPRESS INDEX TO REMOVE THEM TJ050593.93
C GSS9F402.6
JJ =0 GSS9F402.7
DO I=1,IX GSS9F402.8
IF (FIELD(I).NE.RMDI) THEN GSS9F402.9
JJ =JJ+1 GSS9F402.10
IGATH1(JJ)=I GSS9F402.11
END IF GSS9F402.12
END DO GSS9F402.13
C GSS9F402.14
NMIS=IX-JJ COEX1A.389
C COEX1A.390
IF(JJ.NE.0)THEN COEX1A.391
C REMOVE MISSING DATA COEX1A.392
CDIR$ IVDEP COEX1A.393
DO I =1,JJ COEX1A.394
ATEMP(I)=FIELD(IGATH1(I)) COEX1A.395
END DO COEX1A.396
C COEX1A.397
C GET NO. OF NON-ZERO (NON-RMDI) POINTS + COMPRESS INDEX TO REMOVE THEM TJ050593.95
C GSS9F402.16
NPOINT=0 GSS9F402.17
DO I =1,JJ GSS9F402.18
IF (ATEMP(I).NE.0.0) THEN GSS9F402.19
NPOINT =NPOINT+1 GSS9F402.20
IGATH2(NPOINT)=I GSS9F402.21
END IF GSS9F402.22
END DO GSS9F402.23
C GSS9F402.24
NZERO=JJ-NPOINT COEX1A.400
C COEX1A.401
C SET BASE VALUE COEX1A.402
DO I =1,JJ COEX1A.403
IF(ATEMP(I).LT.BASE) BASE=ATEMP(I) COEX1A.404
END DO COEX1A.405
ENDIF COEX1A.406
C COEX1A.407
C CHECK IF A BITMAP FOR MISSING DATA IS NEEDED, COEX1A.408
C COEX1A.409
IF (NMIS.GT.0) THEN COEX1A.410
OBTMIS = .TRUE. COEX1A.411
ELSE COEX1A.412
OBTMIS = .FALSE. COEX1A.413
END IF COEX1A.414
C COEX1A.415
C ROUND BASE TO PRECISION REQUIRED COEX1A.416
C COEX1A.417
IF (BASE.NE.0.0) THEN COEX1A.418
BASE = BASE*BPREC COEX1A.419
IBASE = NINT(BASE) COEX1A.420
BASE = IBASE*APREC COEX1A.421
ELSE COEX1A.422
IBASE=0 COEX1A.423
END IF COEX1A.424
! Fujitsu vectorization directive GRB0F405.199
!OCL NOVREC GRB0F405.200
C COEX1A.425
C FIND DIFFERENCE FROM BASE AND SCALE COEX1A.426
C FIND MAXIMUM DIFFERENCE COEX1A.427
C COEX1A.428
IMAX = 0 COEX1A.429
DO 20 I = 1,JJ COEX1A.430
ATEMP(I) = ATEMP(I)*BPREC COEX1A.431
ITEMP(I) = NINT(ATEMP(I)) - IBASE COEX1A.432
IF(ITEMP(I).LT.0) ITEMP(I) = 0 COEX1A.433
IF (IMAX.LT.ITEMP(I)) IMAX = ITEMP(I) COEX1A.434
20 CONTINUE COEX1A.435
C COEX1A.436
C FIND NUMBER OF BITS REQUIRED TO STORE MAX DIFFERENCE COEX1A.437
C COEX1A.438
IBIT = 0 COEX1A.439
IF (IMAX.GT.0) THEN COEX1A.440
I2 = 1 COEX1A.441
DO WHILE(IMAX.GE.I2) COEX1A.442
IBIT = IBIT + 1 COEX1A.443
I2 = I2*2 COEX1A.444
ENDDO COEX1A.445
ENDIF COEX1A.446
C COEX1A.447
C SET START POSITION IN OUTPUT ARRAY COEX1A.448
C COEX1A.449
JWORD = 1 COEX1A.450
JBIT = 63 COEX1A.451
C COEX1A.452
C IF BIT-MAPPING FOR MISSING DATA THEN INSERT BITMAP COEX1A.453
C COEX1A.454
IF (OBTMIS) THEN COEX1A.455
DO 30 I = 1,IX COEX1A.456
IF (IMIS(I).EQ.1) THEN COEX1A.457
ICOMP(JWORD) = IBSET(ICOMP(JWORD),JBIT) COEX1A.458
IMAP(I) = 0 COEX1A.459
END IF COEX1A.460
IF (JBIT.EQ.0) THEN COEX1A.461
JBIT = 63 COEX1A.462
JWORD = JWORD + 1 COEX1A.463
ELSE COEX1A.464
JBIT = JBIT - 1 COEX1A.465
END IF COEX1A.466
30 CONTINUE COEX1A.467
END IF COEX1A.468
C COEX1A.469
C IF WORTHWHILE USE BIT MAP AND COMPRESS OUT ZEROS. COEX1A.470
C COEX1A.471
IF (IBIT.GT.0) THEN COEX1A.472
IF (NZERO.GT.IX/IBIT) THEN COEX1A.473
OBTZER = .TRUE. COEX1A.474
DO 40 I = 1,IX COEX1A.475
IF (IZERO(I).EQ.1) THEN COEX1A.476
ICOMP(JWORD) = IBCLR(ICOMP(JWORD),JBIT) COEX1A.477
IMAP(I) = 0 COEX1A.478
ELSE COEX1A.479
ICOMP(JWORD) = IBSET(ICOMP(JWORD),JBIT) COEX1A.480
END IF COEX1A.481
IF (JBIT.EQ.0) THEN COEX1A.482
JBIT = 63 COEX1A.483
JWORD = JWORD + 1 COEX1A.484
ELSE COEX1A.485
JBIT = JBIT - 1 COEX1A.486
END IF COEX1A.487
40 CONTINUE COEX1A.488
ELSE COEX1A.489
OBTZER = .FALSE. COEX1A.490
END IF COEX1A.491
C COEX1A.492
C IF BIT MAP INCLUDED FILL TO END OF CURRENT 32 BIT WORD COEX1A.493
C AND SET POINTERS TO START OF NEXT 32 BIT BOUNDARY. COEX1A.494
C COEX1A.495
IF (OBTZER .OR. OBTMIS) THEN COEX1A.496
DO 50 I = 0,JBIT COEX1A.497
ICOMP(JWORD) = IBSET(ICOMP(JWORD),I) COEX1A.498
50 CONTINUE COEX1A.499
IF (JBIT.NE.63) THEN COEX1A.500
IF (JBIT.GE.31) THEN COEX1A.501
JBIT = 31 COEX1A.502
ELSE COEX1A.503
JWORD = JWORD + 1 COEX1A.504
JBIT = 63 COEX1A.505
ENDIF COEX1A.506
ENDIF COEX1A.507
ELSE COEX1A.508
JBIT = 63 COEX1A.509
END IF COEX1A.510
C COEX1A.511
C IF BIT MAPPING ZEROS - COMPRESS OUT UNWANTED ZERO VALUES COEX1A.512
C (OTHERWISE PACK ALL NON RMDI DATA (JJ) ) TJ050593.96
C COEX1A.514
IF (OBTZER) THEN COEX1A.515
CDIR$ IVDEP COEX1A.516
DO I= 1,NPOINT COEX1A.517
ITEMP(I)=ITEMP(IGATH2(I)) COEX1A.518
END DO COEX1A.519
ELSE COEX1A.520
NPOINT=JJ COEX1A.521
END IF COEX1A.522
C COEX1A.523
C MOVE INTO OUTPUT ARRAY USING MINIMUM NUMBER OF BITS REQUIRED COEX1A.524
C COEX1A.525
DO 70 I = 1,NPOINT COEX1A.526
CALL INSTIN
(ICOMP(JWORD),2,JBIT,IBIT,ITEMP(I),0) COEX1A.527
JBIT = JBIT - IBIT COEX1A.528
IF (JBIT.LT.0) THEN COEX1A.529
JWORD = JWORD + 1 COEX1A.530
JBIT = JBIT + 64 COEX1A.531
END IF COEX1A.532
70 CONTINUE COEX1A.533
ELSEIF(IBIT.EQ.0) THEN COEX1A.534
C COEX1A.535
C IF BIT MAP INCLUDED FILL TO END OF CURRENT 32 BIT WORD COEX1A.536
C COEX1A.537
IF (OBTMIS) THEN COEX1A.538
DO 80 I = 0,JBIT COEX1A.539
ICOMP(JWORD) = IBSET(ICOMP(JWORD),I) COEX1A.540
80 CONTINUE COEX1A.541
END IF COEX1A.542
END IF COEX1A.543
C COEX1A.544
C CALCULATE LENGTH IN 32 BIT WORDS COEX1A.545
C COEX1A.546
NOP = JWORD*64 - JBIT - 1 COEX1A.547
NOP = (NOP+31)/32 COEX1A.548
C COEX1A.549
C SET FLAGS TO INDICATE BIT MAPS COEX1A.550
C COEX1A.551
IF (OBTZER) IBIT = IBIT + 128 COEX1A.552
IF (OBTMIS) IBIT = IBIT + 32 COEX1A.553
COEX1A.554
RETURN COEX1A.555
! Fujitsu vectorization directive GRB0F405.201
!OCL NOVREC GRB0F405.202
COEX1A.556
END COEX1A.557
SUBROUTINE XPND(IX,ICOMP,FIELD,APREC,IBIT,BASE,NOP,RMDI) 1,1TJ050593.97
CFPP$ NOCONCUR R COEX1A.559
CFPP$ EXPAND (EXTRIN) COEX1A.560
INTEGER IX,NOP,ICOMP(NOP+1),IBIT COEX1A.561
REAL FIELD(IX),APREC,BASE,RMDI TJ050593.98
INTEGER I,JWORD,JBIT,JJ,NPOINT,NZERO COEX1A.563
INTEGER IMAP(IX),IMIS(IX),IZERO(IX),IMIN(IX),ITEMP(IX) COEX1A.564
INTEGER IGATH(IX) COEX1A.565
REAL ATEMP(IX) COEX1A.566
LOGICAL OBTMIN,OBTMIS,OBTZER,OBTMAP,BTEST COEX1A.567
EXTERNAL INSTIN,EXTRIN COEX1A.568
C COEX1A.569
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.570
C COEX1A.571
C IX = LENTH OF FIELD COEX1A.572
C ICOMP = DATA FOR EXPANSION COEX1A.573
C FIELD = FIELD OF EXPANDED DATA COEX1A.574
C APREC = PRECISION COEX1A.575
C IBIT = NUMBER OF BITS INTO WHICH DATA PACKED COEX1A.576
C BASE = REFERENCE (MIN) VALUE FOR PACKED DATA COEX1A.577
C NOP = SIZE OF COMP COEX1A.578
C RMDI = MISSING DATA INDICATOR TJ050593.99
C COEX1A.580
C COEX1A.581
C INITIALISE VARIABLES/TEMP ARRAYS COEX1A.582
C COEX1A.583
OBTMAP = .FALSE. COEX1A.584
OBTMIS = .FALSE. COEX1A.585
OBTMIN = .FALSE. COEX1A.586
OBTZER = .FALSE. COEX1A.587
DO I = 1,IX COEX1A.588
IMAP(I) = 1 COEX1A.589
IMIS(I) = 0 COEX1A.590
IMIN(I) = 0 COEX1A.591
IZERO(I) = 0 COEX1A.592
END DO COEX1A.593
C COEX1A.594
C CHECK IF BITMAP USED FOR ZERO VALUES COEX1A.595
C COEX1A.596
IF (IBIT.GE.128) THEN COEX1A.597
OBTZER = .TRUE. COEX1A.598
OBTMAP = .TRUE. COEX1A.599
IBIT = IBIT-128 COEX1A.600
ELSE COEX1A.601
OBTMIN = .FALSE. COEX1A.602
ENDIF COEX1A.603
C COEX1A.604
C CHECK IF BITMAP USED FOR MINIMUM VALUES (NOT CURRENTLY USED) COEX1A.605
C COEX1A.606
IF (IBIT.GE.64) THEN COEX1A.607
OBTMIN = .TRUE. COEX1A.608
OBTMAP = .TRUE. COEX1A.609
IBIT = IBIT - 64 COEX1A.610
ELSE COEX1A.611
OBTMIN = .FALSE. COEX1A.612
ENDIF COEX1A.613
C COEX1A.614
C CHECK IF BITMAP USED FOR MISSING DATA COEX1A.615
C COEX1A.616
IF (IBIT.GE.32) THEN COEX1A.617
OBTMIS = .TRUE. COEX1A.618
OBTMAP = .TRUE. COEX1A.619
IBIT = IBIT - 32 COEX1A.620
ELSE COEX1A.621
OBTMIS = .FALSE. COEX1A.622
ENDIF COEX1A.623
C COEX1A.624
C SET START POSITION IN ICOMP COEX1A.625
C COEX1A.626
JWORD = 1 COEX1A.627
JBIT = 63 COEX1A.628
C COEX1A.629
C EXTRACT BITMAPS COEX1A.630
C COEX1A.631
IF (OBTMIS) THEN COEX1A.632
C COEX1A.633
C EXTRACT MISSING DATA BITMAP COEX1A.634
C COEX1A.635
DO 10 I=1,IX COEX1A.636
IF (BTEST(ICOMP(JWORD),JBIT)) THEN COEX1A.637
IMIS(I) = 1 COEX1A.638
IMAP(I) = 0 COEX1A.639
ELSE COEX1A.640
IMIS(I) = 0 COEX1A.641
ENDIF COEX1A.642
IF(JBIT.GT.0) THEN COEX1A.643
JBIT = JBIT - 1 COEX1A.644
ELSE COEX1A.645
JBIT = 63 COEX1A.646
JWORD = JWORD + 1 COEX1A.647
ENDIF COEX1A.648
10 CONTINUE COEX1A.649
ENDIF COEX1A.650
IF (OBTMIN) THEN COEX1A.651
C COEX1A.652
C EXTRACT MINIMUM VALUE BITMAP (NOT USED AT PRESENT) COEX1A.653
C COEX1A.654
DO 20 I=1,IX COEX1A.655
IF(BTEST(ICOMP(JWORD),JBIT)) THEN COEX1A.656
IMIN(I) = 1 COEX1A.657
IMAP(I) = 0 COEX1A.658
ELSE COEX1A.659
IMIN(I) = 0 COEX1A.660
ENDIF COEX1A.661
IF(JBIT.GT.0) THEN COEX1A.662
JBIT = JBIT - 1 COEX1A.663
ELSE COEX1A.664
JBIT = 63 COEX1A.665
JWORD = JWORD + 1 COEX1A.666
ENDIF COEX1A.667
20 CONTINUE COEX1A.668
ENDIF COEX1A.669
IF (OBTZER) THEN COEX1A.670
C COEX1A.671
C EXTRACT ZERO VALUE BITMAP COEX1A.672
C COEX1A.673
DO 30 I=1,IX COEX1A.674
IF(BTEST(ICOMP(JWORD),JBIT)) THEN COEX1A.675
IZERO(I)= 0 COEX1A.676
ELSE COEX1A.677
IZERO(I)= 1 COEX1A.678
IMAP(I) = 0 COEX1A.679
ENDIF COEX1A.680
IF(JBIT.GT.0) THEN COEX1A.681
JBIT = JBIT - 1 COEX1A.682
ELSE COEX1A.683
JBIT = 63 COEX1A.684
JWORD = JWORD + 1 COEX1A.685
ENDIF COEX1A.686
30 CONTINUE COEX1A.687
ENDIF COEX1A.688
C COEX1A.689
C IF BIT MAP USED FIND NUMBER OF POINTS STORED COEX1A.690
C AND RESET POINTERS TO BEGINNING OF 32 BIT BOUNDARY COEX1A.691
C COEX1A.692
IF (OBTMAP) THEN COEX1A.693
NPOINT = 0 COEX1A.694
DO 40 I=1,IX COEX1A.695
IF (IMAP(I).EQ.1) NPOINT = NPOINT + 1 COEX1A.696
40 CONTINUE COEX1A.697
IF (JBIT.NE.63) THEN COEX1A.698
IF (JBIT.GE.31) THEN COEX1A.699
JBIT = 31 COEX1A.700
ELSE COEX1A.701
JBIT = 63 COEX1A.702
JWORD = JWORD + 1 COEX1A.703
ENDIF COEX1A.704
ENDIF COEX1A.705
ELSE COEX1A.706
NPOINT = IX COEX1A.707
ENDIF COEX1A.708
IF (IBIT.GT.0) THEN COEX1A.709
C COEX1A.710
C UNPACK SCALED VALUES TO TEMP ARRAY COEX1A.711
C COEX1A.712
DO 50 I=1,NPOINT COEX1A.713
CALL EXTRIN
(ICOMP(JWORD),2,JBIT,IBIT,ITEMP(I),0) COEX1A.714
JBIT = JBIT - IBIT COEX1A.715
IF (JBIT.LT.0) THEN COEX1A.716
JWORD = JWORD + 1 COEX1A.717
JBIT = JBIT + 64 COEX1A.718
ENDIF COEX1A.719
50 CONTINUE COEX1A.720
C COEX1A.721
C ADD DIFFERENCES TO MINIMUM VALUE AND UNSCALE COEX1A.722
C COEX1A.723
DO 60 I=1,NPOINT COEX1A.724
ATEMP(I)=ITEMP(I)*APREC+BASE COEX1A.725
60 CONTINUE COEX1A.726
C COEX1A.727
C MOVE INTO UNPACKED ARRAY FIELD COEX1A.728
C COEX1A.729
C FIRST GET GATHER INDEX COEX1A.730
C GSS9F402.26
JJ =0 GSS9F402.27
DO I=1,IX GSS9F402.28
IF (IMAP(I).EQ.1) THEN GSS9F402.29
JJ =JJ+1 GSS9F402.30
IGATH(JJ)=I GSS9F402.31
END IF GSS9F402.32
END DO GSS9F402.33
C GSS9F402.34
DO I=1,IX COEX1A.733
FIELD(I)=0. COEX1A.734
END DO COEX1A.735
COEX1A.736
DO I=1,JJ COEX1A.737
FIELD(IGATH(I)) = ATEMP(I) COEX1A.738
END DO COEX1A.739
C COEX1A.740
C IF MINIMUMS BIT MAPPED FILL ZEROS IN FIELD WITH BASE COEX1A.741
C COEX1A.742
IF (OBTMIN) THEN COEX1A.743
DO 80 I=1,IX COEX1A.744
IF(IMIN(I).EQ.1) FIELD(I) = BASE COEX1A.745
80 CONTINUE COEX1A.746
ENDIF COEX1A.747
C COEX1A.748
C IF MISSING DATA BIT MAPPED FILL ZEROS IN FIELD WITH MDI COEX1A.749
C COEX1A.750
IF (OBTMIS) THEN COEX1A.751
DO 90 I=1,IX COEX1A.752
IF(IMIS(I).EQ.1) FIELD(I) = RMDI TJ050593.100
90 CONTINUE COEX1A.754
ENDIF COEX1A.755
COEX1A.756
ELSEIF (IBIT.EQ.0) THEN COEX1A.757
COEX1A.758
C COEX1A.759
C ALL POINTS IN ROW HAVE SAME VALUE E.G. POLE ROW COEX1A.760
C COEX1A.761
DO 100 I=1,IX COEX1A.762
FIELD(I)=BASE COEX1A.763
100 CONTINUE COEX1A.764
C COEX1A.765
C IF MISSING DATA BIT MAPPED FILL ZEROS IN FIELD WITH MDI COEX1A.766
C COEX1A.767
IF (OBTMIS) THEN COEX1A.768
DO 110 I=1,IX COEX1A.769
IF(IMIS(I).EQ.1) FIELD(I) = RMDI TJ050593.101
110 CONTINUE COEX1A.771
ENDIF COEX1A.772
ENDIF COEX1A.773
C COEX1A.774
RETURN COEX1A.775
COEX1A.776
END COEX1A.777
SUBROUTINE INSTIN(ICOMP,IWORD,ISTART,NBIT,INUM,ISIGN) 1COEX1A.778
CFPP$ NOCONCUR R COEX1A.779
INTEGER IWORD,ICOMP(IWORD),ISTART,NBIT,INUM,ISIGN COEX1A.780
INTEGER ISNUM,ISCOMP,NUM COEX1A.781
EXTERNAL MOVBIT COEX1A.782
C COEX1A.783
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.784
C COEX1A.785
C COEX1A.786
C SUBROUTINE TO INSERT AN INTEGER VALUE INTO NBITS OF ICOMP COEX1A.787
C ========================================================= COEX1A.788
C COEX1A.789
C ICOMP = FULLWORD ARRAY TO WHICH AN NBIT INTEGER IS TO BE ADDED COEX1A.790
C IWORD = DIMESION OF ICOMP (ENABLES INTEGER TO CROSS WORDS) COEX1A.791
C ISTART = POSITION IN ICOMP WHERE AN NBIT INTEGER TO BE ADDED COEX1A.792
C NBIT = NUMBER OF BITS INTO WHICH INTEGER IS TO BE STORED COEX1A.793
C INUM = FULLWORD INTEGER TO BE ADDED TO ICOMP COEX1A.794
C ISIGN = INDICATOR FOR STORAGE IN ICOMP: 0= +VE INTEGER COEX1A.795
C 1= SIGNED INTEGER COEX1A.796
C 2= SIGN BIT & +VE INT. COEX1A.797
C COEX1A.798
C START POSITIONS IN ICOMP COEX1A.799
C | | | | COEX1A.800
C | ^ | ^ | ^ | ^ ^ COEX1A.801
C 6666555555555544444444443333333333222222222211111111110000000000 COEX1A.802
C 3210987654321098765432109876543210987654321098765432109876543210 COEX1A.803
C COEX1A.804
C 0000000001111111111222222222233333333334444444444555555555566666 COEX1A.805
C 1234567890123456789012345678901234567890123456789012345678901234 COEX1A.806
C COEX1A.807
C COEX1A.808
ISNUM = 64 - NBIT + 1 COEX1A.809
ISCOMP = 64 - ISTART COEX1A.810
NUM = NBIT COEX1A.811
C COEX1A.812
C MOVE SIGN BIT COEX1A.813
C COEX1A.814
IF (ISIGN.EQ.2) THEN COEX1A.815
CALL MOVBIT(
INUM,1,1,ICOMP(1),ISCOMP) COEX1A.816
NUM = NBIT - 1 COEX1A.817
ISNUM = ISNUM + 1 COEX1A.818
ISCOMP = ISCOMP + 1 COEX1A.819
INUM = -INUM COEX1A.820
ENDIF COEX1A.821
C COEX1A.822
C MOVE INTEGER COEX1A.823
C COEX1A.824
CALL MOVBIT(
INUM,ISNUM,NUM,ICOMP(1),ISCOMP) COEX1A.825
C COEX1A.826
RETURN COEX1A.827
COEX1A.828
END COEX1A.829
SUBROUTINE EXTRIN(ICOMP,IWORD,ISTART,NBIT,INUM,ISIGN) 1COEX1A.830
CFPP$ NOCONCUR R COEX1A.831
INTEGER IWORD,ICOMP(IWORD),ISTART,NBIT,INUM,ISIGN COEX1A.832
INTEGER ISCOMP,ISNUM,NUM,IBIT COEX1A.833
LOGICAL BTEST COEX1A.834
INTEGER IBSET,IBCLR COEX1A.835
EXTERNAL MOVBIT COEX1A.836
C COEX1A.837
C CRAY VERSION C.1.1 16/11/90 P J SMITH COEX1A.838
C COEX1A.839
C COEX1A.840
C SUBROUTINE TO EXTRACT AN INTEGER VALUE FROM NBITS OF ICOMP COEX1A.841
C ========================================================== COEX1A.842
C COEX1A.843
C ICOMP = ARRAY FROM WHICH AN NBIT INTEGER IS TO BE RETRIEVED COEX1A.844
C IWORD = DIMESION OF ICOMP (ENABLES INTEGER TO CROSS WORDS) COEX1A.845
C ISTART = POSITION IN ICOMP WHERE AN NBIT INTEGER STARTS COEX1A.846
C NBIT = NUMBER OF BITS IN WHICH INTEGER HAS BEEN STORED COEX1A.847
C INUM = INTEGER EXTRACTED FROM ICOMP COEX1A.848
C ISIGN = INDICATOR FOR STORAGE IN ICOMP: 0= +VE INTEGER COEX1A.849
C 1= SIGNED INTEGER COEX1A.850
C 2= SIGN BIT & +VE INT. COEX1A.851
C COEX1A.852
C START POSITIONS IN ICOMP COEX1A.853
C | | | | COEX1A.854
C | ^ | ^ | ^ | ^ ^ COEX1A.855
C 6666555555555544444444443333333333222222222211111111110000000000 COEX1A.856
C 3210987654321098765432109876543210987654321098765432109876543210 COEX1A.857
C COEX1A.858
C 0000000001111111111222222222233333333334444444444555555555566666 COEX1A.859
C 1234567890123456789012345678901234567890123456789012345678901234 COEX1A.860
C COEX1A.861
C COEX1A.862
INUM = 0 COEX1A.863
ISCOMP = 64 - ISTART COEX1A.864
ISNUM = 64 - NBIT + 1 COEX1A.865
NUM = NBIT COEX1A.866
C COEX1A.867
C MOVE INTEGER COEX1A.868
C COEX1A.869
IF (ISIGN.EQ.0) THEN COEX1A.870
C COEX1A.871
C POSITIVE INTEGER COEX1A.872
C COEX1A.873
CALL MOVBIT(
ICOMP,ISCOMP,NUM,INUM,ISNUM) COEX1A.874
COEX1A.875
ELSEIF(ISIGN.EQ.1) THEN COEX1A.876
C COEX1A.877
C SIGNED INTEGER COEX1A.878
C COEX1A.879
C SIGN COEX1A.880
CALL MOVBIT(
ICOMP,ISCOMP,1,INUM,1) COEX1A.881
ISCOMP = ISCOMP + 1 COEX1A.882
ISNUM = ISNUM +1 COEX1A.883
NUM = NUM -1 COEX1A.884
C INTEGER COEX1A.885
CALL MOVBIT(
ICOMP,ISCOMP,NUM,INUM,ISNUM) COEX1A.886
C SET UNDIFINED IF INUM NEGATIVE COEX1A.887
IF (INUM.LT.0) THEN COEX1A.888
DO IBIT=NUM,63 COEX1A.889
INUM = IBSET(INUM,IBIT) COEX1A.890
END DO COEX1A.891
ENDIF COEX1A.892
COEX1A.893
ELSEIF(ISIGN.EQ.2) THEN COEX1A.894
C COEX1A.895
C SIGN BIT PLUS POSITIVE INTEGER COEX1A.896
C COEX1A.897
ISCOMP = ISCOMP + 1 COEX1A.898
ISNUM = ISNUM +1 COEX1A.899
NUM = NUM -1 COEX1A.900
C INTEGER COEX1A.901
CALL MOVBIT(
ICOMP,ISCOMP,NUM,INUM,ISNUM) COEX1A.902
C SIGN COEX1A.903
IF (BTEST(ICOMP(1),ISTART)) INUM = -INUM COEX1A.904
ENDIF COEX1A.905
COEX1A.906
RETURN COEX1A.907
COEX1A.908
END COEX1A.909
*IF DEF,T3E,AND,DEF,MPP GBCQF405.61
GBCQF405.62
SUBROUTINE MPP_COEX( 1,1GBCQF405.63
& FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI, GBCQF405.64
& LOCAL_LEVS,GLOBAL_LEVS,PE_FOR_LEVEL,LOCAL_LEVEL) GBCQF405.65
GBCQF405.66
IMPLICIT NONE GBCQF405.67
GBCQF405.68
INTEGER GBCQF405.69
& M ! IN: horizontal size of field GBCQF405.70
&, N ! IN: size of compressed data GBCQF405.71
&, IX ! IN: X dimension of field GBCQF405.72
&, IY ! IN: Y dimension of field GBCQF405.73
&, NUM ! OUT: Total number of compressed (32 bit) words output GBCQF405.74
&, ISC ! IN: Accuracy in power of 2 GBCQF405.75
GBCQF405.76
&, LOCAL_LEVS ! IN: number of levels held on this PE GBCQF405.77
&, GLOBAL_LEVS ! IN: total number of levels being COEX'd GBCQF405.78
GBCQF405.79
LOGICAL GBCQF405.80
& OCO ! IN: .TRUE.=compress .FALSE.=expand (unsupported) GBCQF405.81
GBCQF405.82
INTEGER GBCQF405.83
& ICOMP(N,LOCAL_LEVS) ! OUT: returned compressed data GBCQF405.84
GBCQF405.85
&, PE_FOR_LEVEL(GLOBAL_LEVS) ! IN: which PE holds each level GBCQF405.86
&, LOCAL_LEVEL(GLOBAL_LEVS) ! IN: local level number on that PE GBCQF405.87
GBCQF405.88
REAL GBCQF405.89
& FIELD(IX,IY,LOCAL_LEVS) ! IN: field to compress GBCQF405.90
&, RMDI ! IN: real missing data indicator GBCQF405.91
GBCQF405.92
! COMMON blocks and PARAMETERS GBCQF405.93
*CALL PARVARS
GBCQF405.94
GBCQF405.95
INTEGER min_rows_per_pe ! minimum number of rows per proc. GBCQF405.96
PARAMETER(min_rows_per_pe=1) GBCQF405.97
GBCQF405.98
! Local variables GBCQF405.99
GBCQF405.100
INTEGER K,JJ,IERR GBCQF405.101
GBCQF405.102
INTEGER GBCQF405.103
& total_rows ! total number of rows to be COEXd GBCQF405.104
&, rows_per_pe ! number of rows per PE when split up GBCQF405.105
&, n_full_pes ! number of PEs doing rows_per_pe rows GBCQF405.106
! ! all the others do rows_per_pe-1 rows GBCQF405.107
&, first_gr_index ! first global row index I will COEX GBCQF405.108
&, n_rows ! number of rows I will be COEXing GBCQF405.109
GBCQF405.110
! Functions GBCQF405.111
INTEGER CRI2IBM GBCQF405.112
GBCQF405.113
! Check not been called for expand (not yet supported) GBCQF405.114
GBCQF405.115
IF (.NOT. OCO) THEN GBCQF405.116
WRITE(6,*) 'MPP_COEX called for expand'
GBCQF405.117
WRITE(6,*) 'Not yet supported.' GBCQF405.118
GOTO 9999 GBCQF405.119
ENDIF GBCQF405.120
GBCQF405.121
! Calculate scale factor,columns and rows. GBCQF405.122
GBCQF405.123
IERR=CRI2IBM(2,1,ICOMP(1,1),32,ISC,1,64,32) GBCQF405.124
IERR=CRI2IBM(2,1,ICOMP(2,1),0,IX,1,64,16) GBCQF405.125
IERR=CRI2IBM(2,1,ICOMP(2,1),16,IY,1,64,16) GBCQF405.126
GBCQF405.127
! Copy these values to other levels GBCQF405.128
GBCQF405.129
DO K=2,LOCAL_LEVS GBCQF405.130
ICOMP(1,K)=ICOMP(1,1) GBCQF405.131
ICOMP(2,K)=ICOMP(2,1) GBCQF405.132
ENDDO GBCQF405.133
GBCQF405.134
! Calculate mapping GBCQF405.135
! Distribute each row of every level to a given processor GBCQF405.136
GBCQF405.137
total_rows=IY*GLOBAL_LEVS GBCQF405.138
GBCQF405.139
rows_per_pe=((total_rows-1)/nproc)+1 GBCQF405.140
GBCQF405.141
IF (rows_per_pe .GE. min_rows_per_pe) THEN GBCQF405.142
GBCQF405.143
! However, this can give a few too many rows, so only n_full_pes GBCQF405.144
! PEs will do rows_per_pe, and all the rest will do rows_per_pe-1 GBCQF405.145
GBCQF405.146
n_full_pes=nproc-(rows_per_pe*nproc - total_rows) GBCQF405.147
GBCQF405.148
! Calculate the first global_row_index that I am responsible for GBCQF405.149
! and the number of rows I will COEX GBCQF405.150
GBCQF405.151
IF (mype .LE. n_full_pes) THEN GBCQF405.152
first_gr_index=1+(mype*rows_per_pe) GBCQF405.153
ELSE GBCQF405.154
first_gr_index=1+(n_full_pes*rows_per_pe)+ GBCQF405.155
& (mype-n_full_pes)*(rows_per_pe-1) GBCQF405.156
ENDIF GBCQF405.157
GBCQF405.158
IF (mype .LT. n_full_pes) THEN GBCQF405.159
n_rows=rows_per_pe GBCQF405.160
ELSE GBCQF405.161
n_rows=rows_per_pe-1 GBCQF405.162
ENDIF GBCQF405.163
GBCQF405.164
ELSE ! we force some PE's to do min_rows_per_pe, others none GBCQF405.165
GBCQF405.166
IF ((mype*min_rows_per_pe) .GE. total_rows) THEN GBCQF405.167
n_rows=0 GBCQF405.168
first_gr_index=1 GBCQF405.169
ELSE GBCQF405.170
n_rows=MIN(min_rows_per_pe, GBCQF405.171
& (total_rows-(mype*min_rows_per_pe))) GBCQF405.172
first_gr_index=1+(mype*min_rows_per_pe) GBCQF405.173
ENDIF GBCQF405.174
GBCQF405.175
ENDIF GBCQF405.176
GBCQF405.177
CALL MPP_COEX2
( GBCQF405.178
& FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI, GBCQF405.179
& LOCAL_LEVS,GLOBAL_LEVS,PE_FOR_LEVEL,LOCAL_LEVEL, GBCQF405.180
& first_gr_index,n_rows) GBCQF405.181
GBCQF405.182
9999 CONTINUE GBCQF405.183
GBCQF405.184
RETURN GBCQF405.185
END GBCQF405.186
GBCQF405.187
SUBROUTINE MPP_COEX2( 1,1GBCQF405.188
& FIELD,M,ICOMP,N,IX,IY,NUM,ISC,OCO,RMDI, GBCQF405.189
& LOCAL_LEVS,GLOBAL_LEVS,PE_FOR_LEVEL,LOCAL_LEVEL, GBCQF405.190
& first_gr_index,n_rows) GBCQF405.191
GBCQF405.192
IMPLICIT NONE GBCQF405.193
GBCQF405.194
INTEGER GBCQF405.195
& M ! IN: horizontal size of field GBCQF405.196
&, N ! IN: size of compressed data GBCQF405.197
&, IX ! IN: X dimension of field GBCQF405.198
&, IY ! IN: Y dimension of field GBCQF405.199
&, NUM ! OUT: Total number of compressed (32 bit) words output GBCQF405.200
&, ISC ! IN: Accuracy in power of 2 GBCQF405.201
GBCQF405.202
&, LOCAL_LEVS ! IN: number of levels held on this PE GBCQF405.203
&, GLOBAL_LEVS ! IN: total number of levels being COEX'd GBCQF405.204
GBCQF405.205
&, first_gr_index ! IN: first global row index I will COEX GBCQF405.206
&, n_rows ! IN: number of rows I will COEX GBCQF405.207
GBCQF405.208
LOGICAL GBCQF405.209
& OCO ! IN: .TRUE.=compress .FALSE.=expand (unsupported) GBCQF405.210
GBCQF405.211
INTEGER GBCQF405.212
& ICOMP(N,LOCAL_LEVS) ! OUT: returned compressed data GBCQF405.213
GBCQF405.214
&, PE_FOR_LEVEL(GLOBAL_LEVS) ! IN: which PE holds each level GBCQF405.215
&, LOCAL_LEVEL(GLOBAL_LEVS) ! IN: local level number on that PE GBCQF405.216
GBCQF405.217
REAL GBCQF405.218
& FIELD(IX,IY,LOCAL_LEVS) ! IN: field to compress GBCQF405.219
&, RMDI ! IN: real missing data indicator GBCQF405.220
GBCQF405.221
! COMMON blocks and PARAMETERS GBCQF405.222
*CALL PARVARS
GBCQF405.223
GBCQF405.224
! Local variables GBCQF405.225
GBCQF405.226
INTEGER K,KK,JJ,II,IERR GBCQF405.227
GBCQF405.228
INTEGER GBCQF405.229
& gri ! global row index counter GBCQF405.230
&, local_row ! local row corresponding to gri GBCQF405.231
&, local_lev ! local level corresponding to gri GBCQF405.232
&, global_lev ! global level corresponding to gri GBCQF405.233
&, proc ! PE corresponding to gri GBCQF405.234
&, start_row ! first row in chunk to process GBCQF405.235
&, last_row ! last row this PE is responsible for GBCQF405.236
&, n_rows_in_chunk ! number of rows in chunk GBCQF405.237
GBCQF405.238
INTEGER GBCQF405.239
& IST,ICX,JCX,JBIT,JLEN,NOB,NOC GBCQF405.240
&, IC(IY),ICB(IY),ISTART(IY) GBCQF405.241
&, NOP(IY,LOCAL_LEVS),NOP_PROC(IY) GBCQF405.242
&, IBIT(IY,LOCAL_LEVS),IBIT_PROC GBCQF405.243
&, JCOMP(IX,IY,LOCAL_LEVS),JCOMP_PROC(IX,IY) GBCQF405.244
GBCQF405.245
REAL GBCQF405.246
& ACC,APREC GBCQF405.247
&, BASE(IY,LOCAL_LEVS),BASE_PROC GBCQF405.248
&, FIELD_PROC(IX,IY) GBCQF405.249
GBCQF405.250
! Shmem addressing stuff GBCQF405.251
INTEGER GBCQF405.252
& nvars GBCQF405.253
GBCQF405.254
PARAMETER(nvars=3) GBCQF405.255
GBCQF405.256
INTEGER GBCQF405.257
& address_list(nvars) GBCQF405.258
&, remote_address_list(nvars) GBCQF405.259
GBCQF405.260
COMMON /coex_shmem_align/ GBCQF405.261
& address_list GBCQF405.262
GBCQF405.263
INTEGER GBCQF405.264
& remote_FIELD(IX,IY,GLOBAL_LEVS) GBCQF405.265
&, remote_JCOMP(IX,IY,GLOBAL_LEVS) GBCQF405.266
&, remote_NOP(IY,GLOBAL_LEVS) GBCQF405.267
GBCQF405.268
POINTER GBCQF405.269
& (ptr_remote_FIELD,remote_FIELD), GBCQF405.270
& (ptr_remote_JCOMP,remote_JCOMP), GBCQF405.271
& (ptr_remote_NOP,remote_NOP) GBCQF405.272
GBCQF405.273
EQUIVALENCE GBCQF405.274
& (ptr_remote_FIELD,remote_address_list(1)), GBCQF405.275
& (ptr_remote_JCOMP,remote_address_list(2)), GBCQF405.276
& (ptr_remote_NOP,remote_address_list(3)) GBCQF405.277
GBCQF405.278
!---End of Shmem addressing stuff GBCQF405.279
GBCQF405.280
! Functions GBCQF405.281
INTEGER CRI2IBM GBCQF405.282
GBCQF405.283
! Check not been called for expand (not yet supported) GBCQF405.284
GBCQF405.285
IF (.NOT. OCO) THEN GBCQF405.286
WRITE(6,*) 'MPP_COEX called for expand'
GBCQF405.287
WRITE(6,*) 'Not yet supported.' GBCQF405.288
GOTO 9999 GBCQF405.289
ENDIF GBCQF405.290
GBCQF405.291
! Initialise output array (compressed data might not fill it) GBCQF405.292
GBCQF405.293
DO K=1,LOCAL_LEVS GBCQF405.294
DO JJ=1,IY GBCQF405.295
DO II=1,IX GBCQF405.296
JCOMP(II,JJ,K)=0.0 GBCQF405.297
ENDDO GBCQF405.298
ENDDO GBCQF405.299
ENDDO GBCQF405.300
GBCQF405.301
! Initialise SHMEM addressing stuff. GBCQF405.302
GBCQF405.303
address_list(1)=LOC(FIELD) GBCQF405.304
address_list(2)=LOC(JCOMP) GBCQF405.305
address_list(3)=LOC(NOP) GBCQF405.306
GBCQF405.307
CALL barrier(
) GBCQF405.308
GBCQF405.309
! Do the COEX GBCQF405.310
GBCQF405.311
ACC=ISC GBCQF405.312
APREC=2.0**ACC GBCQF405.313
GBCQF405.314
IF (n_rows .NE. 0) THEN GBCQF405.315
GBCQF405.316
gri=first_gr_index ! global row index GBCQF405.317
global_lev=((gri-1)/IY)+1 ! global level number GBCQF405.318
local_row=gri-(global_lev-1)*IY ! local row number GBCQF405.319
local_lev=LOCAL_LEVEL(global_lev) ! local level number GBCQF405.320
proc=PE_FOR_LEVEL(global_lev) ! PE this level is stored on GBCQF405.321
GBCQF405.322
! Get addresses of array on processor PROC GBCQF405.323
GBCQF405.324
CALL shmem_get(
remote_address_list,address_list, GBCQF405.325
& nvars,proc) GBCQF405.326
GBCQF405.327
! Implicit relationship by EQUIVALENCE GBCQF405.328
! ptr_remote_FIELD=remote_address_list(1) GBCQF405.329
! ptr_remote_JCOMP=remote_address_list(2) GBCQF405.330
! ptr_remote_NOP=remote_address_list(3) GBCQF405.331
GBCQF405.332
! Loop over the rows I am responsible for, getting the data into GBCQF405.333
! FIELD_PROC from the remote FIELD array, COEXing it, and then GBCQF405.334
! sending the compressed data back. GBCQF405.335
GBCQF405.336
start_row=gri GBCQF405.337
last_row=gri+n_rows-1 GBCQF405.338
GBCQF405.339
DO WHILE (start_row .LE. last_row) GBCQF405.340
GBCQF405.341
! Process a chunk of contiguous rows. GBCQF405.342
! Either do all the rows to the end of the level - or just up GBCQF405.343
! to last_row if this is less. GBCQF405.344
GBCQF405.345
n_rows_in_chunk= GBCQF405.346
& MIN(((global_lev*IY)-start_row+1), GBCQF405.347
& (last_row-start_row+1)) GBCQF405.348
GBCQF405.349
CALL shmem_get(
FIELD_PROC, GBCQF405.350
& remote_FIELD(1,local_row,local_lev), GBCQF405.351
& IX*n_rows_in_chunk,proc) GBCQF405.352
GBCQF405.353
! DO JJ=start_row,start_row+n_rows_in_chunk-1 GBCQF405.354
DO JJ=1,n_rows_in_chunk GBCQF405.355
GBCQF405.356
CALL CMPS
(IX,FIELD_PROC(1,JJ),JCOMP_PROC(2,JJ), GBCQF405.357
& NOP_PROC(JJ),APREC, GBCQF405.358
& IBIT_PROC,BASE_PROC,RMDI) GBCQF405.359
GBCQF405.360
! Encode base value, number of bits and length into JCOMP GBCQF405.361
GBCQF405.362
IERR=CRI2IBM(3,1,JCOMP_PROC(1,JJ),0, BASE_PROC,1,64,32) GBCQF405.363
IERR=CRI2IBM(2,1,JCOMP_PROC(1,JJ),32,IBIT_PROC,1,64,16) GBCQF405.364
IERR=CRI2IBM(2,1,JCOMP_PROC(1,JJ),48,NOP_PROC(JJ), 1,64,16) GBCQF405.365
GBCQF405.366
ENDDO ! JJ GBCQF405.367
GBCQF405.368
! Send the compressed data back to the PE from whence it came GBCQF405.369
GBCQF405.370
CALL shmem_put(
remote_JCOMP(1,local_row,local_lev), GBCQF405.371
& JCOMP_PROC, GBCQF405.372
& IX*n_rows_in_chunk,proc) GBCQF405.373
GBCQF405.374
CALL shmem_put(
remote_NOP(local_row,local_lev), GBCQF405.375
& NOP_PROC, GBCQF405.376
& n_rows_in_chunk,proc) GBCQF405.377
GBCQF405.378
! Increment the row counter, and move to the next level, unless GBCQF405.379
! we've now finished all our work GBCQF405.380
GBCQF405.381
start_row=start_row+n_rows_in_chunk GBCQF405.382
GBCQF405.383
IF (start_row .LE. last_row) THEN GBCQF405.384
local_row=1 GBCQF405.385
global_lev=global_lev+1 GBCQF405.386
local_lev=LOCAL_LEVEL(global_lev) GBCQF405.387
proc=PE_FOR_LEVEL(global_lev) GBCQF405.388
GBCQF405.389
CALL shmem_get(
remote_address_list,address_list, GBCQF405.390
& nvars,proc) GBCQF405.391
GBCQF405.392
! Implicit relation from EQUIVALENCE statement GBCQF405.393
! ptr_remote_FIELD=remote_address_list(1) GBCQF405.394
! ptr_remote_JCOMP=remote_address_list(2) GBCQF405.395
! ptr_remote_NOP=remote_address_list(3) GBCQF405.396
GBCQF405.397
ENDIF GBCQF405.398
GBCQF405.399
ENDDO ! DO WHILE over chunks GBCQF405.400
GBCQF405.401
ENDIF ! IF n_rows .NE. 0 GBCQF405.402
GBCQF405.403
CALL barrier(
) ! wait for all my compressed data to arrive back GBCQF405.404
GBCQF405.405
! Calculate positions in output array for packed row GBCQF405.406
! (First packed row starts at word 1, bit 31) GBCQF405.407
GBCQF405.408
DO KK=1,GLOBAL_LEVS GBCQF405.409
GBCQF405.410
IF (PE_FOR_LEVEL(KK) .EQ. mype) THEN GBCQF405.411
GBCQF405.412
K=LOCAL_LEVEL(KK) GBCQF405.413
GBCQF405.414
IC(1)=2 GBCQF405.415
ICB(1)=-1 GBCQF405.416
ISTART(1)=5 GBCQF405.417
GBCQF405.418
DO JJ=2,IY GBCQF405.419
IF (MOD(NOP(JJ-1,K),2) .EQ. 1) THEN GBCQF405.420
IC(JJ ) = IC(JJ-1) + NOP(JJ-1,K)/2 + 1 GBCQF405.421
ICB(JJ) = -ICB(JJ-1) GBCQF405.422
IF (ICB(JJ).GT.0) IC(JJ) = IC(JJ) + 1 GBCQF405.423
ELSE GBCQF405.424
IC(JJ) = IC(JJ-1) + (NOP(JJ-1,K)+1)/2 + 1 GBCQF405.425
ICB(JJ) = ICB(JJ-1) GBCQF405.426
ENDIF GBCQF405.427
ISTART(JJ) = 5 GBCQF405.428
IF(ICB(JJ).EQ.1) ISTART(JJ) = 1 GBCQF405.429
ENDDO GBCQF405.430
GBCQF405.431
! Move temporary array into output array GBCQF405.432
GBCQF405.433
DO JJ=1,IY GBCQF405.434
NOB = NOP(JJ,K)*4 + 8 GBCQF405.435
IST = ISTART(JJ) GBCQF405.436
ICX = IC(JJ) GBCQF405.437
CALL STRMOV(
JCOMP(1,JJ,K),1,NOB,ICOMP(ICX,K),IST) GBCQF405.438
ENDDO GBCQF405.439
GBCQF405.440
! Insert total length of this field GBCQF405.441
GBCQF405.442
NUM = IC(IY)*2 + NOP(IY,K) GBCQF405.443
IF (ICB(IY).LT.0) NUM = NUM + 1 GBCQF405.444
IERR=CRI2IBM(2,1,ICOMP(1,K),0,NUM,1,64,32) GBCQF405.445
GBCQF405.446
ENDIF ! If this level is on this PE GBCQF405.447
GBCQF405.448
ENDDO ! K loop over local levels GBCQF405.449
GBCQF405.450
9999 CONTINUE GBCQF405.451
GBCQF405.452
RETURN GBCQF405.453
END GBCQF405.454
GBCQF405.455
*ENDIF GBCQF405.456
*ENDIF COEX1A.910