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