*IF DEF,C93_1A,OR,DEF,C93_2A GNF0F402.6
C ******************************COPYRIGHT****************************** GTS2F400.1819
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.1820
C GTS2F400.1821
C Use, duplication or disclosure of this code is subject to the GTS2F400.1822
C restrictions as set forth in the contract. GTS2F400.1823
C GTS2F400.1824
C Meteorological Office GTS2F400.1825
C London Road GTS2F400.1826
C BRACKNELL GTS2F400.1827
C Berkshire UK GTS2F400.1828
C RG12 2SZ GTS2F400.1829
C GTS2F400.1830
C If no contract has been raised with this copy of the code, the use, GTS2F400.1831
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.1832
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.1833
C Modelling at the above address. GTS2F400.1834
C ******************************COPYRIGHT****************************** GTS2F400.1835
C GTS2F400.1836
CLL SUBROUTINE DEL SQUARED FFT P ------------------------------------- DEL2P1A.3
CLL DEL2P1A.4
CLL PURPOSE: Solves DEL SQUARED Q = RHS on surface of sphere. DEL2P1A.5
CLL Uses Fast Fourier Fransforms to decompose right-hand- DEL2P1A.6
CLL side in x-direction and then solves in y-direction DEL2P1A.7
CLL for the fourier coefficients of the solution. DEL2P1A.8
CLL DEL2P1A.9
CLL N.B. 1. This version is for a problem where the solution Q and DEL2P1A.10
CLL the right-hand-side RHS are held at pressure points. DEL2P1A.11
CLL 2. A solution to this equation exists uniquely, upto DEL2P1A.12
CLL an arbitrary constant if and only if the sum of the DEL2P1A.13
CLL right-hand-side values over the sphere is equal to zero. DEL2P1A.14
CLL This is known as the Compatability Condition. DEL2P1A.15
CLL This routine assumes that this condition is satisfied. DEL2P1A.16
CLL This routine chooses the arbitrary constant as follows; DEL2P1A.17
CLL If A(y) is the coefficient of wave number 0 then the DEL2P1A.18
CLL constant is chosen by setting the mean of A to zero. DEL2P1A.19
CLL DEL2P1A.20
CLL NOT SUITABLE FOR SINGLE CLOUMN USE. DEL2P1A.21
CLL VERSION FOR CRAY Y-MP DEL2P1A.22
CLL DEL2P1A.23
CLL M.Mawson <- programmer of some or all of previous code or changes DEL2P1A.24
CLL DEL2P1A.25
CLL Model Modification history from model version 3.0: DEL2P1A.26
CLL version Date DEL2P1A.27
CLL 3.1 6/1/92 FPP directives to fix bug (R.S.Bell) SB210293.1
CLL 3.1 24/02/93 Tidy code to remove QA Fortran messages. MM240293.29
CLL vn4.2 07/11/96:Allow this routine to be used in C93_2A (Farnon) GNF0F402.5
CLL DEL2P1A.28
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4 DEL2P1A.29
CLL STANDARD B. DEL2P1A.30
CLL DEL2P1A.31
CLL SYSTEM COMPONENTS COVERED: DEL2P1A.32
CLL DEL2P1A.33
CLL SYSTEM TASK: DEL2P1A.34
CLL DEL2P1A.35
CLL DOCUMENTATION: U.M.D.P. No 14 by M.H. Mawson DEL2P1A.36
CLL DEL2P1A.37
CLLEND------------------------------------------------------------- DEL2P1A.38
DEL2P1A.39
C DEL2P1A.40
C*L ARGUMENTS:--------------------------------------------------- DEL2P1A.41
SUBROUTINE DEL_SQUARED_FFT_P ,2DEL2P1A.42
1 (Q,RHS,SEC_P_LATITUDE,COS_U_LATITUDE, DEL2P1A.43
2 TRIGS,IFAX,LATITUDE_STEP_INVERSE, DEL2P1A.44
3 EARTH_RADIUS_INVERSE,P_FIELD, DEL2P1A.45
4 ROW_LENGTH) DEL2P1A.46
DEL2P1A.47
IMPLICIT NONE DEL2P1A.48
DEL2P1A.49
INTEGER DEL2P1A.50
& P_FIELD !IN Horizontal dimension of pressure field DEL2P1A.51
&, ROW_LENGTH !IN Number of points on a row. DEL2P1A.52
&, IFAX(10) !IN Holds factors of row_length used in FFT's DEL2P1A.53
DEL2P1A.54
REAL DEL2P1A.55
& RHS(P_FIELD) !IN Holds right-hand-side. DEL2P1A.56
&,SEC_P_LATITUDE(P_FIELD) !IN 1./cos(lat) at pressure points DEL2P1A.57
&,COS_U_LATITUDE(P_FIELD-ROW_LENGTH) !IN cos(lat) at velocity point DEL2P1A.58
&,LATITUDE_STEP_INVERSE !IN. 1/(delta phi) DEL2P1A.59
&,EARTH_RADIUS_INVERSE !IN 1./(radius of earth). DEL2P1A.60
DEL2P1A.61
REAL DEL2P1A.62
& TRIGS(ROW_LENGTH) !IN Holds trigonometric terms used in FFT's DEL2P1A.63
DEL2P1A.64
REAL DEL2P1A.65
& Q(P_FIELD) !OUT Holds solution. DEL2P1A.66
DEL2P1A.67
C*--------------------------------------------------------------------- DEL2P1A.68
DEL2P1A.69
C*L DEFINE ARRAYS AND VARIABLES USED IN THIS ROUTINE----------------- DEL2P1A.70
C DEFINE LOCAL ARRAYS: 5 ARE REQUIRED DEL2P1A.71
DEL2P1A.72
REAL DEL2P1A.73
& RHS_DATA(ROW_LENGTH+2,P_FIELD/ROW_LENGTH) ! Fourier modes of DEL2P1A.74
&, Q_DATA(ROW_LENGTH+2,P_FIELD/ROW_LENGTH) ! Q and RHS. DEL2P1A.75
&, A_DIAG(P_FIELD/ROW_LENGTH,ROW_LENGTH+2) ! Matrix diagonal DEL2P1A.76
&, A_SUB_DIAG(P_FIELD/ROW_LENGTH,ROW_LENGTH+2) ! Sub diagonal DEL2P1A.77
&, A_SUP_DIAG(P_FIELD/ROW_LENGTH,ROW_LENGTH+2) ! Super diagonal DEL2P1A.78
DEL2P1A.79
C*--------------------------------------------------------------------- DEL2P1A.80
C DEFINE LOCAL VARIABLES DEL2P1A.81
DEL2P1A.82
C COUNT VARIABLES FOR DO LOOPS ETC. DEL2P1A.83
INTEGER DEL2P1A.84
& I,J,K,IK DEL2P1A.85
&, ROWS ! Number of rows in field. DEL2P1A.86
&, LOT ! Number of data vectors passed to FFT's. DEL2P1A.87
&, JUMP ! Number of storage locations between the DEL2P1A.88
& ! start of consecutive data vectors. DEL2P1A.89
&, INCREMENT ! Number of storage locations between each DEL2P1A.90
& ! element of the same data vector, 1, if DEL2P1A.91
& ! consecutive. DEL2P1A.92
&, FFT_ISIGN ! Parameter determining whether spectral MM240293.30
& ! to grid-point (1) or grid-point to MM240293.31
& ! spectral (-1) FFT's are required. MM240293.32
DEL2P1A.97
REAL DEL2P1A.98
& SCALAR ! Generic real work variable. DEL2P1A.99
&,FACTOR ! Holds factor in matrix gaussian DEL2P1A.100
& ! elimination. DEL2P1A.101
&,WAVE_NUMBER ! Holds wave number for which fourier DEL2P1A.102
& ! coefficients are being calculated. DEL2P1A.103
DEL2P1A.104
C*L EXTERNAL SUBROUTINE CALLS:------------------------------------ DEL2P1A.105
EXTERNAL FOURIER DEL2P1A.106
C*--------------------------------------------------------------------- DEL2P1A.107
DEL2P1A.108
CL MAXIMUM VECTOR LENGTH ASSUMED IS ROW_LENGTH+2 DEL2P1A.109
CL--------------------------------------------------------------------- DEL2P1A.110
CL INTERNAL STRUCTURE. DEL2P1A.111
CL--------------------------------------------------------------------- DEL2P1A.112
CL DEL2P1A.113
DEL2P1A.114
CL--------------------------------------------------------------------- DEL2P1A.115
CL SECTION 1. Put all RHS data in bigger array before calling DEL2P1A.116
CL fourier decomposition. Two extra addresses per row DEL2P1A.117
CL required. DEL2P1A.118
CL--------------------------------------------------------------------- DEL2P1A.119
DEL2P1A.120
ROWS = P_FIELD/ROW_LENGTH DEL2P1A.121
DEL2P1A.122
CFPP$ NODEPCHK DEL2P1A.123
! Fujitsu vectorization directive GRB0F405.595
!OCL NOVREC GRB0F405.596
DO 100 J=1,ROWS DEL2P1A.124
IK = (J-1)*ROW_LENGTH DEL2P1A.125
DO 110 I=1,ROW_LENGTH DEL2P1A.126
RHS_DATA(I,J) = RHS(IK+I) DEL2P1A.127
110 CONTINUE DEL2P1A.128
100 CONTINUE DEL2P1A.129
DEL2P1A.130
CL--------------------------------------------------------------------- DEL2P1A.131
CL SECTION 2. Call FOURIER to get fourier decomposition of data. DEL2P1A.132
CL--------------------------------------------------------------------- DEL2P1A.133
DEL2P1A.134
INCREMENT = 1 DEL2P1A.135
JUMP = ROW_LENGTH+2 DEL2P1A.136
FFT_ISIGN = -1 MM240293.33
LOT = ROWS DEL2P1A.138
CALL FOURIER
(RHS_DATA,TRIGS,IFAX,INCREMENT,JUMP,ROW_LENGTH, DEL2P1A.139
& LOT,FFT_ISIGN) MM240293.34
DEL2P1A.141
CL--------------------------------------------------------------------- DEL2P1A.142
CL SECTION 3. Solve equation. DEL2P1A.143
CL--------------------------------------------------------------------- DEL2P1A.144
DEL2P1A.145
C --------------------------------------------------------------------- DEL2P1A.146
CL SECTION 3.1 Solve for real constant term. DEL2P1A.147
C --------------------------------------------------------------------- DEL2P1A.148
DEL2P1A.149
CL i) Set up tri-diagonal matrix left-hand-side. DEL2P1A.150
DEL2P1A.151
SCALAR = EARTH_RADIUS_INVERSE*EARTH_RADIUS_INVERSE* DEL2P1A.152
& LATITUDE_STEP_INVERSE*LATITUDE_STEP_INVERSE DEL2P1A.153
J=1 DEL2P1A.154
WAVE_NUMBER = 0. DEL2P1A.155
A_SUB_DIAG(1,J) = 0. DEL2P1A.156
A_SUP_DIAG(ROWS,J) = 0. DEL2P1A.157
DEL2P1A.158
C set up A matrix. DEL2P1A.159
C Northern boundary condition. DEL2P1A.160
I=1 DEL2P1A.161
K = (I-1)*ROW_LENGTH+1 DEL2P1A.162
A_DIAG(I,J) = -2.*SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.163
& COS_U_LATITUDE(K) DEL2P1A.164
A_SUP_DIAG(I,J) = 2.*SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.165
& COS_U_LATITUDE(K) DEL2P1A.166
DEL2P1A.167
C Inner points. DEL2P1A.168
DO 310 I=2,ROWS-1 DEL2P1A.169
K = (I-1)*ROW_LENGTH+1 DEL2P1A.170
A_DIAG(I,J) = -SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.171
& (COS_U_LATITUDE(K)+COS_U_LATITUDE(K-ROW_LENGTH)) DEL2P1A.172
A_SUP_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.173
& COS_U_LATITUDE(K) DEL2P1A.174
A_SUB_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.175
& COS_U_LATITUDE(K-ROW_LENGTH) DEL2P1A.176
310 CONTINUE DEL2P1A.177
DEL2P1A.178
C Southern boundary condition. DEL2P1A.179
I=ROWS DEL2P1A.180
K = (I-1)*ROW_LENGTH+1 DEL2P1A.181
A_DIAG(I,J) = -2.*SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.182
& COS_U_LATITUDE(K-ROW_LENGTH) DEL2P1A.183
A_SUB_DIAG(I,J) = 2.*SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.184
& COS_U_LATITUDE(K-ROW_LENGTH) DEL2P1A.185
DEL2P1A.186
CL ii) Solve matrix system. DEL2P1A.187
DEL2P1A.188
A_DIAG(1,J) = 1./A_DIAG(1,J) DEL2P1A.189
DO 312 I=2,ROWS DEL2P1A.190
FACTOR = A_SUB_DIAG(I,J) * A_DIAG(I-1,J) DEL2P1A.191
A_DIAG(I,J) = 1./(A_DIAG(I,J) - FACTOR*A_SUP_DIAG(I-1,J)) DEL2P1A.192
RHS_DATA(J,I) = RHS_DATA(J,I) - FACTOR*RHS_DATA(J,I-1) DEL2P1A.193
312 CONTINUE DEL2P1A.194
DEL2P1A.195
C Back substitute to get solution. DEL2P1A.196
DEL2P1A.197
Q_DATA(J,ROWS) = A_DIAG(ROWS,J)*RHS_DATA(J,ROWS) DEL2P1A.198
DO 314 I= ROWS-1,1,-1 DEL2P1A.199
Q_DATA(J,I) = A_DIAG(I,J)*(RHS_DATA(J,I)- DEL2P1A.200
& A_SUP_DIAG(I,J)*Q_DATA(J,I+1)) DEL2P1A.201
314 CONTINUE DEL2P1A.202
DEL2P1A.203
C Set constant imaginery mode to zero. DEL2P1A.204
C Remove arbitrary constant from real constant mode. DEL2P1A.205
C Remove mean value as guess to arbitrary constant. DEL2P1A.206
DEL2P1A.207
DO I=1,ROWS DEL2P1A.208
Q_DATA(2,I) = 0. DEL2P1A.209
END DO DEL2P1A.210
SCALAR = 0. DEL2P1A.211
DO I=1,ROWS DEL2P1A.212
SCALAR = SCALAR + Q_DATA(1,I) DEL2P1A.213
END DO DEL2P1A.214
SCALAR = SCALAR / ROWS DEL2P1A.215
CFPP$ NOINNER SB210293.2
DO I=1,ROWS DEL2P1A.216
Q_DATA(1,I) = Q_DATA(1,I) - SCALAR DEL2P1A.217
END DO DEL2P1A.218
DEL2P1A.219
C --------------------------------------------------------------------- DEL2P1A.220
CL SECTION 3.2 Solve for wave-number > 0 modes. DEL2P1A.221
CL Solution at poles is zero and this is thus the DEL2P1A.222
CL boundary condition for the solver. DEL2P1A.223
C --------------------------------------------------------------------- DEL2P1A.224
DEL2P1A.225
CL Set solution at poles = 0 for wave numbers > 0 DEL2P1A.226
DEL2P1A.227
CFPP$ NOINNER SB210293.3
DO I=3,ROW_LENGTH+2 DEL2P1A.228
Q_DATA(I,1) = 0. DEL2P1A.229
Q_DATA(I,ROWS) = 0. DEL2P1A.230
END DO DEL2P1A.231
DEL2P1A.232
CL i) Set up tri-diagonal matrix system left-hand-side. DEL2P1A.233
C Real coefficients only. DEL2P1A.234
DEL2P1A.235
SCALAR = EARTH_RADIUS_INVERSE*EARTH_RADIUS_INVERSE* DEL2P1A.236
& LATITUDE_STEP_INVERSE*LATITUDE_STEP_INVERSE DEL2P1A.237
DO J=3,ROW_LENGTH+2,2 DEL2P1A.238
WAVE_NUMBER = (J+1)/2-1 DEL2P1A.239
A_SUB_DIAG(2,J) = 0. DEL2P1A.240
A_SUP_DIAG(ROWS-1,J) = 0. DEL2P1A.241
DEL2P1A.242
C Set up A matrix. DEL2P1A.243
C Northern boundary condition. DEL2P1A.244
I=2 DEL2P1A.245
K = (I-1)*ROW_LENGTH+1 DEL2P1A.246
DEL2P1A.247
A_SUP_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR*COS_U_LATITUDE(K) DEL2P1A.248
A_DIAG(I,J) = -SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.249
& (COS_U_LATITUDE(K)+COS_U_LATITUDE(K-ROW_LENGTH)) DEL2P1A.250
& - WAVE_NUMBER*WAVE_NUMBER*EARTH_RADIUS_INVERSE DEL2P1A.251
& *EARTH_RADIUS_INVERSE*SEC_P_LATITUDE(K) DEL2P1A.252
& *SEC_P_LATITUDE(K) DEL2P1A.253
DEL2P1A.254
C Inner points. DEL2P1A.255
DO 320 I=3,ROWS-2 DEL2P1A.256
K = (I-1)*ROW_LENGTH+1 DEL2P1A.257
A_SUP_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.258
& COS_U_LATITUDE(K) DEL2P1A.259
A_SUB_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.260
& COS_U_LATITUDE(K-ROW_LENGTH) DEL2P1A.261
A_DIAG(I,J) = -SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.262
& (COS_U_LATITUDE(K)+COS_U_LATITUDE(K-ROW_LENGTH)) DEL2P1A.263
& - WAVE_NUMBER*WAVE_NUMBER*EARTH_RADIUS_INVERSE DEL2P1A.264
& *EARTH_RADIUS_INVERSE*SEC_P_LATITUDE(K) DEL2P1A.265
& *SEC_P_LATITUDE(K) DEL2P1A.266
320 CONTINUE DEL2P1A.267
DEL2P1A.268
C Southern boundary condition. DEL2P1A.269
I=ROWS-1 DEL2P1A.270
K = (I-1)*ROW_LENGTH+1 DEL2P1A.271
A_SUB_DIAG(I,J) = SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.272
& COS_U_LATITUDE(K-ROW_LENGTH) DEL2P1A.273
A_DIAG(I,J) = -SEC_P_LATITUDE(K)*SCALAR* DEL2P1A.274
& (COS_U_LATITUDE(K)+COS_U_LATITUDE(K-ROW_LENGTH)) DEL2P1A.275
& - WAVE_NUMBER*WAVE_NUMBER*EARTH_RADIUS_INVERSE DEL2P1A.276
& *EARTH_RADIUS_INVERSE*SEC_P_LATITUDE(K) DEL2P1A.277
& *SEC_P_LATITUDE(K) DEL2P1A.278
END DO DEL2P1A.279
DEL2P1A.280
C Matrix for imaginery coefficients is copy of real one for same wave DEL2P1A.281
C number. DEL2P1A.282
DO 322 J=4,ROW_LENGTH+2,2 DEL2P1A.283
DO I=2,ROWS-1 DEL2P1A.284
A_DIAG(I,J) = A_DIAG(I,J-1) DEL2P1A.285
A_SUB_DIAG(I,J) = A_SUB_DIAG(I,J-1) DEL2P1A.286
A_SUP_DIAG(I,J) = A_SUP_DIAG(I,J-1) DEL2P1A.287
END DO DEL2P1A.288
322 CONTINUE DEL2P1A.289
DEL2P1A.290
CL ii) Solve matrix system for all right-hand-sides. DEL2P1A.291
DEL2P1A.292
CFPP$ NOINNER SB210293.4
DO 324 J=3,ROW_LENGTH+2 DEL2P1A.293
A_DIAG(2,J) = 1./A_DIAG(2,J) DEL2P1A.294
324 CONTINUE DEL2P1A.295
DO 326 I=3,ROWS-1 DEL2P1A.296
DO J=3,ROW_LENGTH+2 DEL2P1A.297
FACTOR = A_SUB_DIAG(I,J) * A_DIAG(I-1,J) DEL2P1A.298
A_DIAG(I,J) = 1./(A_DIAG(I,J) - FACTOR*A_SUP_DIAG(I-1,J)) DEL2P1A.299
RHS_DATA(J,I) = RHS_DATA(J,I) - FACTOR*RHS_DATA(J,I-1) DEL2P1A.300
END DO DEL2P1A.301
326 CONTINUE DEL2P1A.302
DEL2P1A.303
C Back substitute to get solution. DEL2P1A.304
DEL2P1A.305
DO 328 J=3,ROW_LENGTH+2 DEL2P1A.306
Q_DATA(J,ROWS-1) = A_DIAG(ROWS-1,J)*RHS_DATA(J,ROWS-1) DEL2P1A.307
328 CONTINUE DEL2P1A.308
DO 329 I= ROWS-2,2,-1 DEL2P1A.309
DO J=3,ROW_LENGTH+2 DEL2P1A.310
Q_DATA(J,I) = A_DIAG(I,J)*(RHS_DATA(J,I)- DEL2P1A.311
& A_SUP_DIAG(I,J)*Q_DATA(J,I+1)) DEL2P1A.312
END DO DEL2P1A.313
329 CONTINUE DEL2P1A.314
DEL2P1A.315
CL--------------------------------------------------------------------- DEL2P1A.316
CL SECTION 4. Call FOURIER to create grid-point fields from DEL2P1A.317
CL fourier solution modes. DEL2P1A.318
CL--------------------------------------------------------------------- DEL2P1A.319
DEL2P1A.320
INCREMENT = 1 DEL2P1A.321
JUMP = ROW_LENGTH+2 DEL2P1A.322
FFT_ISIGN = 1 MM240293.35
CALL FOURIER(Q_DATA,TRIGS,IFAX,INCREMENT,JUMP,ROW_LENGTH, DEL2P1A.324
& LOT,FFT_ISIGN) MM240293.36
DEL2P1A.326
CL--------------------------------------------------------------------- DEL2P1A.327
CL SECTION 5. Copy solution into output array. DEL2P1A.328
CL--------------------------------------------------------------------- DEL2P1A.329
DEL2P1A.330
DO 500 J=1,ROWS DEL2P1A.331
IK = (J-1)*ROW_LENGTH DEL2P1A.332
DO 510 I=1,ROW_LENGTH DEL2P1A.333
Q(IK+I) = Q_DATA(I,J) DEL2P1A.334
510 CONTINUE DEL2P1A.335
500 CONTINUE DEL2P1A.336
DEL2P1A.337
CL END OF ROUTINE DEL_SQUARED_FFT_P DEL2P1A.338
DEL2P1A.339
RETURN DEL2P1A.340
END DEL2P1A.341
*ENDIF DEL2P1A.342