SUBROUTINE CG_RELAX( 1,16CGRELAX.32
*CALL ARGSIZE 
CGRELAX.33
*CALL ARGOCALL CGRELAX.34
*CALL ARGOINDX ORH5F402.2
& P,PB,PTD,PTDB ORH3F405.88
&,FDIR,HDIR CGRELAX.37
&,ZTD CGRELAX.38
& ) CGRELAX.39
! for correct operation fdir and hdir should be given as arguments and CGRELAX.40
! thus passed into the dump, and set to zero in the dump for a new run. CGRELAX.41
CGRELAX.42
C----------------------------------------------------------------------- CGRELAX.43
CLL MODEL Date MODIFICATION HISTORY FROM MODEL VERSION 3.4 CGRELAX.44
CLL VERSION CGRELAX.45
CLL 4.1 14.5.96 Complete solver re-write to use CG (KAH) CGRELAX.46
CLL CGRELAX.47
CLL CGRELAX.48
C----------------------------------------------------------------------- CGRELAX.49
C CG_RELAX takes as input the vorticity driving function computed in CGRELAX.50
C "clinic" (ZTD) and, using conjugate gradient iteration CGRELAX.51
C solves the Laplacian equation for the external mode of CGRELAX.52
C velocity in terms of a mass transport stream function (P). CGRELAX.53
C CGRELAX.54
C For information on the method used in this code see: CGRELAX.55
C (for example) R. Barrett et al "Templates for the Solution CGRELAX.56
C of Linear Systems: Building Blocks for Iterative Methods" CGRELAX.57
C Society for Industrial and Applied Mathematics (1994). CGRELAX.58
C CGRELAX.59
C----------------------------------------------------------------------- CGRELAX.60
IMPLICIT NONE CGRELAX.61
CGRELAX.62
C Define global data: CGRELAX.63
*CALL TYPSIZE CGRELAX.64
*CALL TYPOINDX PXORDER.9
*CALL TYPOCALL CGRELAX.65
*CALL UMSCALAR CGRELAX.66
*CALL CNTLOCN CGRELAX.67
*CALL OARRYSIZ CGRELAX.68
*CALL OTIMER CGRELAX.69
CGRELAX.70
C Define argument list: CGRELAX.71
REAL ORH3F405.90
& P(IMT_STREAM,0:JMT_STREAM+1), !stream function ORH1F402.50
& PB(IMT_STREAM,0:JMT_STREAM+1), !stream function at prev tstep ORH1F402.51
& PTD(IMT_STREAM,JMT_STREAM), !chng in strm fn across timestep CGRELAX.77
& PTDB(IMT_STREAM,JMT_STREAM), !chng in strm fn across prvs tstep CGRELAX.78
& ZTD(IMT_STREAM,JMT_STREAM) !Vorticity change across timesteps CGRELAX.79
CGRELAX.80
C Define local variables: CGRELAX.81
REAL FX, ! Temporary value CGRELAX.82
& FXA, ! Temporary value CGRELAX.83
& FXB, ! Temporary value CGRELAX.84
& FXC, ! Temporary value CGRELAX.85
& FXD, ! Temporary value CGRELAX.86
& CRTP, ! Altered val of CRIT (convergence criteria) CGRELAX.87
& RESMAX ! Max of residuals at all grid points CGRELAX.88
CGRELAX.89
INTEGER ORH3F405.89
& I, ! Grid point index (Zonal) CGRELAX.91
& J, ! Grid point index (Meridonal) CGRELAX.92
& L, ! Loop control for ocean segment CGRELAX.93
& IE, ! End index for I loop CGRELAX.94
& IS ! Start index for I loop CGRELAX.95
CGRELAX.96
REAL CGRELAX.97
& RESIS, ! Residual of relaxation of island CGRELAX.98
& RESIS_ROW(JMT), ! Residual of relaxation of island ORH5F402.196
& TEST1, ! Product of reciprocals of 4 surrounding depths CGRELAX.99
& TEST2 ! Sum of reciprocals of 4 surrounding depths CGRELAX.100
CGRELAX.101
REAL CGRELAX.102
& CFN(IMT,JMT), ! coefficients in Laplacian star stencil CGRELAX.103
& CFS(IMT,JMT), ! (north, south, east, west etc) CGRELAX.104
& CFE(IMT,JMT), CGRELAX.105
& CFW(IMT,JMT), CGRELAX.106
& RES(IMT,JMT), ! residual CGRELAX.107
& COF(IMT,JMT), ! normalising coefficient CGRELAX.108
& COFIS(NISLE), CGRELAX.109
& COFIS_ROW(JMT), ORH5F402.197
& CPF(IMT,JMT), ! normalisation factor for Laplacian stencil CGRELAX.110
& H(IMT,JMT) ! total depth at UV points CGRELAX.111
CGRELAX.112
INTEGER CGRELAX.113
& ISLE, ! Do loop index indicates current island CGRELAX.114
& JE, ! End index for J loop CGRELAX.115
& JS, ! Start index for J loop CGRELAX.116
& N, ! Do loop index for island segments CGRELAX.117
& NSEG, ! End index for island segment DO loops CGRELAX.118
& ISMASK(IMT,JMT) ! grid point type indicator CGRELAX.119
! (0=interior, 1=perimeter, 2=land) CGRELAX.120
& ,INFO ! MPP call info ORH5F402.4
CGRELAX.121
C Definitions for Conjugate gradient: CGRELAX.122
REAL CGRELAX.123
& GDIR(IMT,JMT), ! CG gradient directions etc CGRELAX.124
& FDIR(IMT,JMT), CGRELAX.125
& HDIR(IMT,JMT) CGRELAX.126
REAL ASUM(JMT) ! used for temporary sums CGRELAX.127
REAL beta, ! Conjugate Gradient scalars CGRELAX.128
& asor, CGRELAX.129
& oapr, CGRELAX.130
& apr2, CGRELAX.131
& apr4 CGRELAX.132
CGRELAX.133
C Flags used: (Pulled in from an include deck or whatever) CGRELAX.134
C L_OTIMER - Ocean Timer CGRELAX.135
C L_OCYCLIC - Ocean Cyclic CGRELAX.136
C L_OISLANDS - Ocean Islands CGRELAX.137
C L_OROTATE - Ocean Rotation CGRELAX.138
C L_OSYMM - Ocean Symmetry CGRELAX.139
C MIX - 1/0 its a mixing timestep CGRELAX.140
CGRELAX.141
CGRELAX.142
C======================================================================= CGRELAX.143
IF (L_OTIMER) CALL TIMER('CGRELAX ',3) CGRELAX.144
CGRELAX.145
C Begin introductory section to prepare for the relaxation: CGRELAX.146
CGRELAX.147
C Initialize the work area: CGRELAX.148
DO J=J_1,J_JMT ORH5F402.5
DO I=1,IMT CGRELAX.150
CFN(I,J)=0.0 CGRELAX.151
CFS(I,J)=0.0 CGRELAX.152
CFE(I,J)=0.0 CGRELAX.153
CFW(I,J)=0.0 CGRELAX.154
CPF(I,J)=0.0 CGRELAX.155
ASUM(J)=0.0 ORH5F402.198
IF( L_OISLANDS ) THEN CGRELAX.156
COF(I,J)=1.0 CGRELAX.157
ISMASK(I,J)=0 CGRELAX.158
ENDIF CGRELAX.159
ENDDO CGRELAX.160
ENDDO CGRELAX.161
CGRELAX.162
IF (L_OISLANDS) THEN CGRELAX.163
C Form island mask: (ISMASK=0 over interior ocean points CGRELAX.164
C ISMASK=1 over perimeter ocean points CGRELAX.165
C ISMASK=2 over land points) CGRELAX.166
DO I=2,IMU CGRELAX.167
DO J=J_2,J_JMTM1 ORH5F402.6
TEST1=HR(I,J)*HR(I-1,J)*HR(I,J-1)*HR(I-1,J-1) CGRELAX.169
TEST2=HR(I,J)+HR(I-1,J)+HR(I,J-1)+HR(I-1,J-1) CGRELAX.170
IF(TEST1.EQ.0.0) ISMASK(I,J)=1 CGRELAX.171
IF(TEST2.EQ.0.0) ISMASK(I,J)=2 CGRELAX.172
ENDDO CGRELAX.173
ENDDO CGRELAX.174
CGRELAX.175
ENDIF CGRELAX.176
CGRELAX.177
C Calculate depth field from its reciprocal: CGRELAX.178
C Compute reciprocals with an epsilon added to avoid zero divide: CGRELAX.179
C (never mind the fact that it might be wrong!) CGRELAX.180
DO J=J_2,J_JMTM1 ORH5F402.7
DO I=1,IMT CGRELAX.182
H(I,J)=1.0/(HR(I,J)+1.E-20) CGRELAX.183
ENDDO CGRELAX.184
ENDDO CGRELAX.185
CGRELAX.186
IF (L_OISLANDS) THEN ! reset H over land to zero: CGRELAX.187
DO J=J_1,J_JMT ORH5F402.8
DO I=1,IMT CGRELAX.189
IF( HR(I,J) .EQ. 0.0) H(I,J)=0.0 CGRELAX.190
ENDDO CGRELAX.191
ENDDO CGRELAX.192
ENDIF CGRELAX.193
CGRELAX.194
IF( L_OCYCLIC )THEN ! set cyclic boundary conditions on H: CGRELAX.195
DO J=J_1,J_JMT ORH5F402.9
H(1,J)=H(IMUM1,J) CGRELAX.197
ENDDO CGRELAX.198
ENDIF CGRELAX.199
CGRELAX.200
*IF DEF,MPP ORH5F402.10
! Populate halo regions for H ORH5F402.11
CALL SWAPBOUNDS(H,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.12
*ENDIF ORH5F402.13
C Generate arrays of coefficients for relaxation: CGRELAX.201
C CGRELAX.202
C 1st, compute coefficients of the Laplacian star: CGRELAX.203
C (hold non-interior points to zero using start and end indices) CGRELAX.204
DO 130 J=J_3,J_JMTM1 ORH5F402.14
DO 130 L=1,LSEG CGRELAX.206
IS=ISZ(J,L) CGRELAX.207
IF(IS.EQ.0) GO TO 130 CGRELAX.208
IE=IEZ(J,L) CGRELAX.209
FXA=2.0*CSTR(J)*CSTR(J) CGRELAX.210
FXB=2.0*CS(J )*CSTR(J)*DYTR(J)*DYUR(J ) CGRELAX.211
FXC=2.0*CS(J-1)*CSTR(J)*DYTR(J)*DYUR(J-1) CGRELAX.212
FX=1.0 CGRELAX.213
DO I=IS,IE CGRELAX.214
CFN(I,J)=FXB/(H(I-1,J)+H(I,J)) CGRELAX.215
CFS(I,J)=FXC/(H(I-1,J-1)+H(I,J-1)) CGRELAX.216
CFE(I,J)=FXA*DXUR(I)*DXTR(I)/(H(I,J)+H(I,J-1)) CGRELAX.217
CFW(I,J)=FXA*DXUR(I-1)*DXTR(I)/(H(I-1,J)+H(I-1,J-1)) CGRELAX.218
CPF(I,J)=FX/(CFN(I,J)+CFS(I,J)+CFE(I,J)+CFW(I,J)) CGRELAX.219
ENDDO CGRELAX.220
CGRELAX.221
C 2nd, augment coefficients for implicit treatment of Coriolis term: CGRELAX.222
IF(ACOR .NE. 0.0) THEN CGRELAX.223
CGRELAX.224
IF (L_OROTATE) THEN CGRELAX.225
FX=-0.5*C2DTSF*ACOR*CSTR(J)*DYTR(J) CGRELAX.226
DO I=IS,IE CGRELAX.227
CFN(I,J)=CFN(I,J) CGRELAX.228
& +(HR(I,J)*CORIOLIS(I,J)-HR(I-1,J)*CORIOLIS(I-1,J)) CGRELAX.229
& *FX*DXTR(I) CGRELAX.230
CFS(I,J)=CFS(I,J) CGRELAX.231
& -(HR(I,J-1)*CORIOLIS(I,J-1)-HR(I-1,J-1)*CORIOLIS(I-1,J-1)) CGRELAX.232
& *FX*DXTR(I) CGRELAX.233
CFE(I,J)=CFE(I,J) CGRELAX.234
& -(HR(I,J)*CORIOLIS(I,J)-HR(I,J-1)*CORIOLIS(I,J-1)) CGRELAX.235
& *FX*DXTR(I) CGRELAX.236
CFW(I,J)=CFW(I,J) CGRELAX.237
& +(HR(I-1,J)*CORIOLIS(I-1,J)-HR(I-1,J-1)*CORIOLIS(I-1,J-1)) CGRELAX.238
& *FX*DXTR(I) CGRELAX.239
ENDDO ! over I CGRELAX.240
ELSE CGRELAX.241
FX=-C2DTSF*ACOR*CSTR(J)*DYTR(J)*OMEGA CGRELAX.242
DO I=IS,IE CGRELAX.243
CGRELAX.244
CFN(I,J)=CFN(I,J)+(HR(I,J )-HR(I-1,J ))*SINE(J ) CGRELAX.245
& *FX*DXTR(I) CGRELAX.246
CFS(I,J)=CFS(I,J)-(HR(I,J-1)-HR(I-1,J-1))*SINE(J-1) CGRELAX.247
& *FX*DXTR(I) CGRELAX.248
CFE(I,J)=CFE(I,J)-(HR(I ,J)*SINE(J)-HR(I ,J-1)*SINE(J-1)) CGRELAX.249
& *FX*DXTR(I) CGRELAX.250
CFW(I,J)=CFW(I,J)+(HR(I-1,J)*SINE(J)-HR(I-1,J-1)*SINE(J-1)) CGRELAX.251
& *FX*DXTR(I) CGRELAX.252
ENDDO ! over I CGRELAX.253
ENDIF ! L_OROTATE = false CGRELAX.254
CGRELAX.255
ENDIF CGRELAX.256
CGRELAX.257
C 3rd, normalize coefficients and forcing term for efficiency: CGRELAX.258
DO I=IS,IE CGRELAX.259
CFN(I,J)=CFN(I,J)*CPF(I,J) CGRELAX.260
CFS(I,J)=CFS(I,J)*CPF(I,J) CGRELAX.261
CFE(I,J)=CFE(I,J)*CPF(I,J) CGRELAX.262
CFW(I,J)=CFW(I,J)*CPF(I,J) CGRELAX.263
ENDDO CGRELAX.264
130 CONTINUE ! end of J and L loops CGRELAX.265
IF (L_OISLANDS) THEN CGRELAX.266
CGRELAX.267
C 4th, compute coefficients on island perimeter points: CGRELAX.268
DO 180 ISLE=1,NISLE CGRELAX.269
DO J=J_1,J_JMT ORH5F402.199
COFIS_ROW(J)=0.0 ORH5F402.200
ENDDO ORH5F402.201
COFIS(ISLE)=0.0 CGRELAX.270
NSEG=ISEG(ISLE) CGRELAX.271
DO 175 N=1,NSEG CGRELAX.272
IS=ISIS(ISLE,N) CGRELAX.273
IE=IEIS(ISLE,N) CGRELAX.274
JS=JSIS(ISLE,N) CGRELAX.275
JE=JEIS(ISLE,N) CGRELAX.276
*IF DEF,MPP ORH5F402.15
c Adjust J indices to be relative to local values for each PE ORH5F402.16
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.17
c Make sure the J loop doesn't get executed ORH5F402.18
JS = 1 ORH5F402.19
JE = 0 ORH5F402.20
ELSE ORH5F402.21
IF (JST.GT.JS) THEN ORH5F402.22
JS = JST ORH5F402.23
ENDIF ORH5F402.24
IF (JFIN.LT.JE) THEN ORH5F402.25
JE = JFIN ORH5F402.26
ENDIF ORH5F402.27
c Adjust to local versions of indices. ORH5F402.28
JS = JS - J_OFFSET ORH5F402.29
JE = JE - J_OFFSET ORH5F402.30
ENDIF ORH5F402.31
*ENDIF ORH5F402.32
DO 175 J=JS,JE CGRELAX.277
FXA=2.0*CSTR(J)*CSTR(J) CGRELAX.278
FXB=2.0*CS(J )*DYUR(J )*DYTR(J)*CSTR(J) CGRELAX.279
FXC=2.0*CS(J-1)*DYUR(J-1)*DYTR(J)*CSTR(J) CGRELAX.280
IF (L_OROTATE) THEN CGRELAX.281
FXD=-0.5*C2DTSF*ACOR*CSTR(J)*DYTR(J) CGRELAX.282
ELSE CGRELAX.283
FXD=-C2DTSF*ACOR*CSTR(J)*DYTR(J)*OMEGA CGRELAX.284
ENDIF CGRELAX.285
CGRELAX.286
DO 175 I=IS,IE CGRELAX.287
IF(ISMASK(I,J).NE.1) GO TO 174 CGRELAX.288
CGRELAX.289
IF (L_OROTATE) THEN CGRELAX.290
CGRELAX.291
IF(HR(I-1,J ).NE.0.0 .OR. HR(I ,J ).NE.0.0) CGRELAX.292
& CFN(I,J)=FXB/(H(I-1,J)+H(I,J)) CGRELAX.293
CGRELAX.294
IF(HR(I-1,J-1).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.295
& CFS(I,J)=FXC/(H(I-1,J-1)+H(I,J-1)) CGRELAX.296
CGRELAX.297
IF(HR(I ,J ).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.298
& CFE(I,J)=FXA*DXTR(I)*DXUR(I)/(H(I,J)+H(I,J-1)) CGRELAX.299
CGRELAX.300
IF(HR(I-1,J ).NE.0.0 .OR. HR(I-1,J-1).NE.0.0) CGRELAX.301
& CFW(I,J)=FXA*DXTR(I)*DXUR(I-1)/(H(I-1,J)+H(I-1,J-1)) CGRELAX.302
ELSE CGRELAX.303
CGRELAX.304
IF(HR(I-1,J ).NE.0.0 .OR. HR(I ,J ).NE.0.0) CGRELAX.305
& CFN(I,J)=FXB/(H(I-1,J)+H(I,J)) CGRELAX.306
CGRELAX.307
IF(HR(I-1,J-1).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.308
& CFS(I,J)=FXC/(H(I-1,J-1)+H(I,J-1)) CGRELAX.309
CGRELAX.310
IF(HR(I ,J ).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.311
& CFE(I,J)=FXA*DXTR(I)*DXUR(I)/(H(I,J)+H(I,J-1)) CGRELAX.312
CGRELAX.313
IF(HR(I-1,J ).NE.0.0 .OR. HR(I-1,J-1).NE.0.0) CGRELAX.314
& CFW(I,J)=FXA*DXTR(I)*DXUR(I-1)/(H(I-1,J)+H(I-1,J-1)) CGRELAX.315
ENDIF CGRELAX.316
ORH5F402.202
CPF(I,J) = 1.0/(CFN(i,j)+CFS(I,J)+CFE(I,J)+CFW(I,J)) CGRELAX.318
CGRELAX.319
IF (L_OROTATE) THEN CGRELAX.320
CGRELAX.321
IF(HR(I-1,J ).NE.0.0 .OR. HR(I ,J ).NE.0.0) CGRELAX.322
& CFN(I,J)=FXB/(H(I-1,J)+H(I,J)) CGRELAX.323
& +FXD*DXTR(I)*(HR(I,J)*CORIOLIS(I,J) CGRELAX.324
& -HR(I-1,J)*CORIOLIS(I-1,J)) CGRELAX.325
CGRELAX.326
IF(HR(I-1,J-1).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.327
& CFS(I,J)=FXC/(H(I-1,J-1)+H(I,J-1)) CGRELAX.328
& -FXD*DXTR(I)*(HR(I,J-1)*CORIOLIS(I,J-1) CGRELAX.329
& -HR(I-1,J-1)*CORIOLIS(I-1,J-1)) CGRELAX.330
CGRELAX.331
IF(HR(I ,J ).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.332
& CFE(I,J)=FXA*DXTR(I)*DXUR(I)/(H(I,J)+H(I,J-1)) CGRELAX.333
& -FXD*DXTR(I)*(HR(I,J)*CORIOLIS(I,J) CGRELAX.334
& -HR(I,J-1)*CORIOLIS(I,J-1)) CGRELAX.335
CGRELAX.336
IF(HR(I-1,J ).NE.0.0 .OR. HR(I-1,J-1).NE.0.0) CGRELAX.337
& CFW(I,J)=FXA*DXTR(I)*DXUR(I-1)/(H(I-1,J)+H(I-1,J-1)) CGRELAX.338
& +FXD*DXTR(I)*(HR(I-1,J)*CORIOLIS(I-1,J) CGRELAX.339
& -HR(I-1,J-1)*CORIOLIS(I-1,J-1)) CGRELAX.340
ELSE CGRELAX.341
CGRELAX.342
IF(HR(I-1,J ).NE.0.0 .OR. HR(I ,J ).NE.0.0) CGRELAX.343
& CFN(I,J)=FXB/(H(I-1,J)+H(I,J)) CGRELAX.344
& +FXD*DXTR(I)*(HR(I,J)-HR(I-1,J))*SINE(J) CGRELAX.345
CGRELAX.346
IF(HR(I-1,J-1).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.347
& CFS(I,J)=FXC/(H(I-1,J-1)+H(I,J-1)) CGRELAX.348
& -FXD*DXTR(I)*(HR(I,J-1)-HR(I-1,J-1))*SINE(J-1) CGRELAX.349
CGRELAX.350
IF(HR(I ,J ).NE.0.0 .OR. HR(I ,J-1).NE.0.0) CGRELAX.351
& CFE(I,J)=FXA*DXTR(I)*DXUR(I)/(H(I,J)+H(I,J-1)) CGRELAX.352
& -FXD*DXTR(I)*(HR(I,J)*SINE(J)-HR(I,J-1)*SINE(J-1)) CGRELAX.353
CGRELAX.354
IF(HR(I-1,J ).NE.0.0 .OR. HR(I-1,J-1).NE.0.0) CGRELAX.355
& CFW(I,J)=FXA*DXTR(I)*DXUR(I-1)/(H(I-1,J)+H(I-1,J-1)) CGRELAX.356
& +FXD*DXTR(I)* CGRELAX.357
& (HR(I-1,J)*SINE(J)-HR(I-1,J-1)*SINE(J-1)) CGRELAX.358
ENDIF CGRELAX.359
CGRELAX.360
CGRELAX.361
CFN(I,J)=CFN(I,J)*CPF(I,J) CGRELAX.362
CFS(I,J)=CFS(I,J)*CPF(I,J) CGRELAX.363
CFE(I,J)=CFE(I,J)*CPF(I,J) CGRELAX.364
CFW(I,J)=CFW(I,J)*CPF(I,J) CGRELAX.365
COF(I,J)=CST(J)*DXT(I)*DYT(J)/CPF(I,J) CGRELAX.366
COFIS_ROW(J)=COFIS_ROW(J)+COF(I,J) ORH5F402.203
174 CONTINUE ! jumped to by a goto above CGRELAX.368
175 CONTINUE ! end of N, J and I loops CGRELAX.369
ORH5F402.204
! Now we must add all COFIS values for this island ORH5F402.205
! in a manner which will allow bit comparison. ORH5F402.206
CALL CGSUMMATION ORH5F402.207
& (COFIS_ROW,JMT_GLOBAL,J_1,J_JMT,JMT, ORH5F402.208
& JST,JFIN,O_MYPE,O_NPROC,COFIS(ISLE)) ORH5F402.209
COFIS(ISLE)=1.0/COFIS(ISLE) CGRELAX.370
ORH5F402.210
180 CONTINUE ! end of isle loop CGRELAX.371
ENDIF ! L_OISLANDS = true CGRELAX.372
CGRELAX.373
DO ISLE=1,NISLE CGRELAX.374
NSEG=ISEG(ISLE) CGRELAX.375
DO N=1,NSEG CGRELAX.376
IS=ISIS(ISLE,N) CGRELAX.377
IE=IEIS(ISLE,N) CGRELAX.378
JS=JSIS(ISLE,N) CGRELAX.379
JE=JEIS(ISLE,N) CGRELAX.380
*IF DEF,MPP ORH5F402.33
c Adjust J indices to be relative to local values for each PE ORH5F402.34
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.35
c Make sure the J loop doesn't get executed ORH5F402.36
JS = 1 ORH5F402.37
JE = 0 ORH5F402.38
ELSE ORH5F402.39
IF (JST.GT.JS) THEN ORH5F402.40
JS = JST ORH5F402.41
ENDIF ORH5F402.42
IF (JFIN.LT.JE) THEN ORH5F402.43
JE = JFIN ORH5F402.44
ENDIF ORH5F402.45
c Adjust to local versions of indices. ORH5F402.46
JS = JS - J_OFFSET ORH5F402.47
JE = JE - J_OFFSET ORH5F402.48
ENDIF ORH5F402.49
*ENDIF ORH5F402.50
DO J=JS,JE CGRELAX.381
DO I=IS,IE CGRELAX.382
IF(ISMASK(I,J).EQ.1) THEN CGRELAX.383
CGRELAX.384
COF(I,J) = COF(I,J)*COFIS(ISLE) CGRELAX.385
ENDIF ORH5F402.211
ENDDO CGRELAX.387
ENDDO CGRELAX.388
ENDDO CGRELAX.389
ENDDO CGRELAX.390
CGRELAX.391
IF (L_OCYCLIC) THEN CGRELAX.392
DO J=J_1,J_JMT ORH5F402.51
COF(1,j) = 0.0 CGRELAX.394
COF(IMT,J) = 0.0 CGRELAX.395
ENDDO CGRELAX.396
ENDIF CGRELAX.397
CGRELAX.398
DO J=J_1,J_JMT ORH5F402.52
DO I = 1, IMT CGRELAX.400
ZTD(I,J)=ZTD(I,J)*CPF(I,J) CGRELAX.401
ENDDO CGRELAX.402
ENDDO CGRELAX.403
C Compute a first guess for the relaxation by extrapolating the two CGRELAX.404
C previous solutions forward in time: CGRELAX.405
C 1st, save values in ptd in array res: CGRELAX.406
C 2nd, perform time extrapolation (accounting for mixing timestep): CGRELAX.407
C 3rd, determine max ptd over all points and use this CGRELAX.408
C to compute criterion for convergence: CGRELAX.409
C 4th, now copy original values from ptd (now in res) into PTDB: CGRELAX.410
C 5th, set res and gdir to zero prior to using them in cg context: CGRELAX.411
CRTP = 0.0 CGRELAX.412
FXA=1.0 CGRELAX.413
FXB=2.0 CGRELAX.414
IF(MIX .NE. 0) FXA=0.5 CGRELAX.415
DO J=J_1,J_JMT ORH5F402.53
DO I=1,IMT CGRELAX.417
res(I,J)=PTD(I,J) CGRELAX.418
PTD(I,J)=FXA * (FXB*PTD(I,J)-PTDB(I,J)) CGRELAX.419
CRTP = MAX(ABS(PTD(I,J)),CRTP) CGRELAX.420
PTDB(I,J)=res(I,J) CGRELAX.421
RES(I,J)=0.0 CGRELAX.422
GDIR(I,J)=0.0 CGRELAX.423
ENDDO CGRELAX.424
ENDDO CGRELAX.425
*IF DEF,MPP ORH5F402.166
CALL GC_RMAX (1,O_NPROC,INFO,CRTP) ORH5F402.167
*ENDIF ORH5F402.168
CRTP=CRTP*CRIT CGRELAX.426
IF (CRTP.LE.0.0) CRTP=1.0 ! covers new runs when ptd is zero CGRELAX.427
CGRELAX.428
CGRELAX.429
C This is the `stop when ||r|| < stop_tol x ||b||' criterion CGRELAX.430
C where stop_tol is passed in from CRIT, and the inf-norms are used CGRELAX.431
CGRELAX.432
beta = 0.0 CGRELAX.433
MSCAN=0 CGRELAX.434
CGRELAX.435
*IF DEF,MPP ORH5F402.172
CALL SWAPBOUNDS(PTD,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.173
*ENDIF ORH5F402.174
C compute initial residuals: CGRELAX.436
DO J=J_2,J_JMTM1 ORH5F402.54
do i=2,imtm1 CGRELAX.438
res(i,j) = cfn(i,j)*ptd(i,j+1) + CGRELAX.439
& cfs(i,j)*ptd(i,j-1) + CGRELAX.440
& cfe(i,j)*ptd(i+1,j) + CGRELAX.441
& cfw(i,j)*ptd(i-1,j) - CGRELAX.442
& ptd(i,j) - ztd(i,j) CGRELAX.443
if( L_OISLANDS ) res(i,j) = res(i,j)*cof(i,j) CGRELAX.444
enddo CGRELAX.445
enddo CGRELAX.446
CGRELAX.447
CGRELAX.448
C Find initial island residual using a line integral: CGRELAX.449
if( L_OISLANDS) then CGRELAX.450
DO ISLE=1,NISLE CGRELAX.451
DO J=J_1,J_JMT ORH5F402.212
RESIS_ROW(J)=0.0 ORH5F402.213
ENDDO ORH5F402.214
RESIS=0.0 CGRELAX.452
NSEG=ISEG(ISLE) CGRELAX.453
DO N=1,NSEG CGRELAX.454
IS=ISIS(ISLE,N) CGRELAX.455
IE=IEIS(ISLE,N) CGRELAX.456
JS=JSIS(ISLE,N) CGRELAX.457
JE=JEIS(ISLE,N) CGRELAX.458
*IF DEF,MPP ORH5F402.55
c Adjust J indices to be relative to local values for each PE ORH5F402.56
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.57
c Make sure the J loop doesn't get executed ORH5F402.58
JS = 1 ORH5F402.59
JE = 0 ORH5F402.60
ELSE ORH5F402.61
IF (JST.GT.JS) THEN ORH5F402.62
JS = JST ORH5F402.63
ENDIF ORH5F402.64
IF (JFIN.LT.JE) THEN ORH5F402.65
JE = JFIN ORH5F402.66
ENDIF ORH5F402.67
c Adjust to local versions of indices. ORH5F402.68
JS = JS - J_OFFSET ORH5F402.69
JE = JE - J_OFFSET ORH5F402.70
ENDIF ORH5F402.71
*ENDIF ORH5F402.72
DO J=JS,JE CGRELAX.459
DO I=IS,IE CGRELAX.460
IF(ISMASK(I,J).EQ.1) THEN CGRELAX.461
resis_ROW(J) = resis_ROW(J) + res(i,j) ORH5F402.215
ENDIF CGRELAX.463
ENDDO ! over I CGRELAX.464
ENDDO ! over J CGRELAX.465
ENDDO CGRELAX.466
ORH5F402.216
CALL CGSUMMATION ORH5F402.217
& (RESIS_ROW,JMT_GLOBAL,J_1,J_JMT,JMT, ORH5F402.218
& JST,JFIN,O_MYPE,O_NPROC,RESIS) ORH5F402.219
ORH5F402.220
DO N=1,NSEG CGRELAX.467
IS=ISIS(ISLE,N) CGRELAX.468
IE=IEIS(ISLE,N) CGRELAX.469
JS=JSIS(ISLE,N) CGRELAX.470
JE=JEIS(ISLE,N) CGRELAX.471
*IF DEF,MPP ORH5F402.73
c Adjust J indices to be relative to local values for each PE ORH5F402.74
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.75
c Make sure the J loop doesn't get executed ORH5F402.76
JS = 1 ORH5F402.77
JE = 0 ORH5F402.78
ELSE ORH5F402.79
IF (JST.GT.JS) THEN ORH5F402.80
JS = JST ORH5F402.81
ENDIF ORH5F402.82
IF (JFIN.LT.JE) THEN ORH5F402.83
JE = JFIN ORH5F402.84
ENDIF ORH5F402.85
c Adjust to local versions of indices. ORH5F402.86
JS = JS - J_OFFSET ORH5F402.87
JE = JE - J_OFFSET ORH5F402.88
ENDIF ORH5F402.89
*ENDIF ORH5F402.90
DO J=JS,JE CGRELAX.472
DO I=IS,IE CGRELAX.473
IF(ISMASK(I,J).EQ.1) THEN CGRELAX.474
res(i,j) = resis CGRELAX.475
ENDIF CGRELAX.476
ENDDO ! over I CGRELAX.477
ENDDO ! over J CGRELAX.478
ENDDO CGRELAX.479
ENDDO CGRELAX.480
endif CGRELAX.481
if( L_OCYCLIC) then ! set cyclic bc CGRELAX.482
DO J=J_2,J_JMTM1 ORH5F402.91
res(1,j) = res(imtm1,j) CGRELAX.484
res(imt,j) = res(2,j) CGRELAX.485
enddo CGRELAX.486
endif CGRELAX.487
*IF DEF,MPP ORH5F402.175
CALL SWAPBOUNDS(RES,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.176
*ENDIF ORH5F402.177
if( L_OSYMM) then ! set symmetry bc CGRELAX.488
do i=1,imt CGRELAX.489
*IF DEF,MPP ORH5F402.138
IF (JST.LE.JMT_GLOBAL.AND.JFIN.GE.JMT_GLOBAL) THEN ORH5F402.139
*ENDIF ORH5F402.140
res(i,j_jmt) = - res(i,j_jmtm1) ORH5F402.141
*IF DEF,MPP ORH5F402.142
ENDIF ORH5F402.143
*ENDIF ORH5F402.144
enddo CGRELAX.491
endif CGRELAX.492
CGRELAX.493
CGRELAX.494
C End introductory section. CGRELAX.495
C======================================================================= CGRELAX.496
CGRELAX.497
C Begin the relaxation: CGRELAX.498
300 CONTINUE CGRELAX.499
MSCAN=MSCAN+1 ! MSCAN will therefore be 1 on first scan. CGRELAX.500
CGRELAX.501
*IF DEF,MPP ORH5F402.178
CALL SWAPBOUNDS(RES,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.179
*ENDIF ORH5F402.180
DO J=J_2,J_JMTM1 ORH5F402.92
do i=2,imtm1 CGRELAX.503
gdir(i,j) = cfn(i,j) * res(i,j+1) + CGRELAX.504
& cfs(i,j) * res(i,j-1) + CGRELAX.505
& cfe(i,j) * res(i+1,j) + CGRELAX.506
& cfw(i,j) * res(i-1,j) - CGRELAX.507
& res(i,j) CGRELAX.508
if( L_OISLANDS) then CGRELAX.509
gdir(i,j) = gdir(i,j)*cof(i,j) CGRELAX.510
endif CGRELAX.511
enddo CGRELAX.512
enddo CGRELAX.513
CGRELAX.514
CGRELAX.515
C Find island residuals using a line integral: CGRELAX.516
C Also update gdir on perimeter & interior of island: CGRELAX.517
if( L_OISLANDS) then CGRELAX.518
CGRELAX.519
DO ISLE=1,NISLE CGRELAX.520
DO J = J_1,J_JMT ORH5F402.221
RESIS_ROW(J) = 0.0 ORH5F402.222
ENDDO ORH5F402.223
RESIS=0.0 CGRELAX.521
NSEG=ISEG(ISLE) CGRELAX.522
DO N=1,NSEG CGRELAX.523
IS=ISIS(ISLE,N) CGRELAX.524
IE=IEIS(ISLE,N) CGRELAX.525
JS=JSIS(ISLE,N) CGRELAX.526
JE=JEIS(ISLE,N) CGRELAX.527
*IF DEF,MPP ORH5F402.93
c Adjust J indices to be relative to local values for each PE ORH5F402.94
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.95
c Make sure the J loop doesn't get executed ORH5F402.96
JS = 1 ORH5F402.97
JE = 0 ORH5F402.98
ELSE ORH5F402.99
IF (JST.GT.JS) THEN ORH5F402.100
JS = JST ORH5F402.101
ENDIF ORH5F402.102
IF (JFIN.LT.JE) THEN ORH5F402.103
JE = JFIN ORH5F402.104
ENDIF ORH5F402.105
c Adjust to local versions of indices. ORH5F402.106
JS = JS - J_OFFSET ORH5F402.107
JE = JE - J_OFFSET ORH5F402.108
ENDIF ORH5F402.109
*ENDIF ORH5F402.110
DO J=JS,JE CGRELAX.528
DO I=IS,IE CGRELAX.529
IF(ISMASK(I,J).EQ.1) THEN CGRELAX.530
resis_ROW(J) = resis_ROW(J) + gdir(i,j) ORH5F402.224
ENDIF CGRELAX.532
ENDDO ! over I CGRELAX.533
ENDDO ! over J CGRELAX.534
ENDDO CGRELAX.535
ORH5F402.225
CALL CGSUMMATION ORH5F402.226
& (RESIS_ROW,JMT_GLOBAL,J_1,J_JMT,JMT, ORH5F402.227
& JST,JFIN,O_MYPE,O_NPROC,RESIS) ORH5F402.228
ORH5F402.229
DO N=1,NSEG CGRELAX.536
IS=ISIS(ISLE,N) CGRELAX.537
IE=IEIS(ISLE,N) CGRELAX.538
JS=JSIS(ISLE,N) CGRELAX.539
JE=JEIS(ISLE,N) CGRELAX.540
*IF DEF,MPP ORH5F402.111
c Adjust J indices to be relative to local values for each PE ORH5F402.112
IF ((JST.GT.JE).OR.(JFIN.LT.JS)) THEN ORH5F402.113
c Make sure the J loop doesn't get executed ORH5F402.114
JS = 1 ORH5F402.115
JE = 0 ORH5F402.116
ELSE ORH5F402.117
IF (JST.GT.JS) THEN ORH5F402.118
JS = JST ORH5F402.119
ENDIF ORH5F402.120
IF (JFIN.LT.JE) THEN ORH5F402.121
JE = JFIN ORH5F402.122
ENDIF ORH5F402.123
c Adjust to local versions of indices. ORH5F402.124
JS = JS - J_OFFSET ORH5F402.125
JE = JE - J_OFFSET ORH5F402.126
ENDIF ORH5F402.127
*ENDIF ORH5F402.128
DO J=JS,JE CGRELAX.541
DO I=IS,IE CGRELAX.542
IF(ISMASK(I,J).EQ.1) THEN CGRELAX.543
gdir(i,j) = resis CGRELAX.544
ENDIF CGRELAX.545
ENDDO ! over I CGRELAX.546
ENDDO ! over J CGRELAX.547
ENDDO CGRELAX.548
ENDDO CGRELAX.549
endif CGRELAX.550
if( L_OCYCLIC) then CGRELAX.551
DO J=J_2,J_JMTM1 ORH5F402.129
gdir(1,j) = gdir(imtm1,j) CGRELAX.553
gdir(imt,j) = gdir(2,j) CGRELAX.554
enddo CGRELAX.555
endif CGRELAX.556
if( L_OSYMM) then ! set symmetry bc CGRELAX.557
*IF DEF,MPP ORH5F402.181
CALL SWAPBOUNDS(GDIR,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.182
*ENDIF ORH5F402.183
do i=1,imt CGRELAX.558
*IF DEF,MPP ORH5F402.145
IF (JST.LE.JMT_GLOBAL.AND.JFIN.GE.JMT_GLOBAL) THEN ORH5F402.146
*ENDIF ORH5F402.147
gdir(i,j_jmt) = - gdir(i,j_jmtm1) ORH5F402.148
*IF DEF,MPP ORH5F402.149
ENDIF ORH5F402.150
*ENDIF ORH5F402.151
enddo CGRELAX.560
*IF DEF,MPP ORH5F402.184
CALL SWAPBOUNDS(GDIR,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.185
*ENDIF ORH5F402.186
endif CGRELAX.561
C gdir now contains the preconditioned solution. CGRELAX.562
CGRELAX.563
C Compute dot product of res and gdir: CGRELAX.564
apr4 = 0.0 CGRELAX.565
DO J=J_2,J_JMTM1 ORH5F402.130
asum(j) = 0.0 CGRELAX.567
do i=2,imtm1 CGRELAX.568
asum(j) = asum(j) + gdir(i,j)*res(i,j) CGRELAX.569
enddo CGRELAX.570
enddo CGRELAX.572
ORH5F402.230
CALL CGSUMMATION ORH5F402.231
& (asum,JMT_GLOBAL,J_1,J_JMT,JMT, ORH5F402.232
& JST,JFIN,O_MYPE,O_NPROC,apr4) ORH5F402.233
CGRELAX.573
CGRELAX.574
if (mscan .gt. 1) beta = apr4/oapr CGRELAX.575
oapr = apr4 CGRELAX.576
CGRELAX.577
C update search directions: CGRELAX.578
DO J=J_2,J_JMTM1 ORH5F402.131
do i=1,imt CGRELAX.580
hdir(i,j) = res(i,j) + beta*hdir(i,j) CGRELAX.581
fdir(i,j) = gdir(i,j) + beta*fdir(i,j) CGRELAX.582
enddo CGRELAX.583
enddo CGRELAX.584
CGRELAX.585
C apr2 is dot product of fdir with itself: CGRELAX.586
apr2 = 0.0 CGRELAX.587
DO J=J_2,J_JMTM1 ORH5F402.132
asum(j) = 0.0 CGRELAX.589
do i=2,imtm1 CGRELAX.590
asum(j) = asum(j) + fdir(i,j)**2 CGRELAX.591
enddo CGRELAX.592
enddo CGRELAX.594
CGRELAX.595
CALL CGSUMMATION ORH5F402.234
& (asum,JMT_GLOBAL,J_1,J_JMT,JMT, ORH5F402.235
& JST,JFIN,O_MYPE,O_NPROC,apr2) ORH5F402.236
CGRELAX.596
C Calculate stepsize (asor) & make a correction to ptd res: CGRELAX.597
asor = - apr4/apr2 CGRELAX.598
CGRELAX.599
CGRELAX.600
DO J=J_2,J_JMTM1 ORH5F402.133
do i=1,imt CGRELAX.602
ptd(i,j) = ptd(i,j) + asor*hdir(i,j) CGRELAX.603
res(i,j) = res(i,j) + asor*fdir(i,j) CGRELAX.604
enddo CGRELAX.605
enddo CGRELAX.606
CGRELAX.607
if( L_OCYCLIC) then ! set cyclic condition CGRELAX.608
DO J=J_2,J_JMTM1 ORH5F402.134
res(1,j) = res(imtm1,j) CGRELAX.610
res(imt,j) = res(2,j) CGRELAX.611
enddo CGRELAX.612
endif CGRELAX.613
CGRELAX.614
C Find the maximum absolute residual to determine convergence: CGRELAX.615
resmax = 0.0 CGRELAX.616
DO J=J_3,jscan ORH5F402.135
do i=2,imtm1 CGRELAX.618
resmax = max(abs(res(i,j)), resmax) CGRELAX.619
enddo CGRELAX.620
enddo CGRELAX.621
*IF DEF,MPP ORH5F402.169
CALL GC_RMAX (1,O_NPROC,INFO,resmax) ORH5F402.170
*ENDIF ORH5F402.171
if (L_OSYMM) then CGRELAX.622
*IF DEF,MPP ORH5F402.187
CALL SWAPBOUNDS(RES,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.188
*ENDIF ORH5F402.189
do i=1,imt CGRELAX.623
*IF DEF,MPP ORH5F402.152
IF (JST.LE.JMT_GLOBAL.AND.JFIN.GE.JMT_GLOBAL) THEN ORH5F402.153
*ENDIF ORH5F402.154
res(i,j_jmt) = - res(i,j_jmtm1) ORH5F402.155
*IF DEF,MPP ORH5F402.156
ENDIF ORH5F402.157
*ENDIF ORH5F402.158
enddo CGRELAX.625
*IF DEF,MPP ORH5F402.190
CALL SWAPBOUNDS(RES,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.191
*ENDIF ORH5F402.192
endif CGRELAX.626
CGRELAX.627
CGRELAX.628
CGRELAX.629
C Do another iteration unless the max residual < crtp or CGRELAX.630
C the number of scans reaches mxscan: CGRELAX.631
if (resmax .ge. crtp .and. mscan .lt. mxscan) go to 300 CGRELAX.632
C end of the iteration loop CGRELAX.633
CGRELAX.634
if (L_OSYMM) then ! set symmetry on northern wall CGRELAX.635
*IF DEF,MPP ORH5F402.193
CALL SWAPBOUNDS(PTD,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH5F402.194
*ENDIF ORH5F402.195
do i=1,imt CGRELAX.636
*IF DEF,MPP ORH5F402.159
IF (JST.LE.JMT_GLOBAL.AND.JFIN.GE.JMT_GLOBAL) THEN ORH5F402.160
*ENDIF ORH5F402.161
ptd(i,j_jmt) = - ptd(i,j_jmtm1) ORH5F402.162
*IF DEF,MPP ORH5F402.163
ENDIF ORH5F402.164
*ENDIF ORH5F402.165
enddo CGRELAX.638
endif CGRELAX.639
CGRELAX.640
C--------------------------------------------------------------------- CGRELAX.641
C Update the stream function based upon the relaxation solution: CGRELAX.642
C (Can use res as a work buffer as we are no longer interested in CGRELAX.643
C its contents for this call.) CGRELAX.644
DO J=J_1,J_JMT ORH5F402.136
DO I=1,IMT CGRELAX.646
res(i,j) = pb(i,j) + ptd(i,j) CGRELAX.647
pb(i,j) = p(i,j) CGRELAX.648
p(i,j) = res(i,j) CGRELAX.649
ENDDO CGRELAX.650
ENDDO CGRELAX.651
C so pb is old value, p is current value, ready for next call or dump CGRELAX.652
CGRELAX.653
C On a mixing timestep, alter PTD to be consistent with CGRELAX.654
C normal, leap-frog stepping: CGRELAX.655
IF( MIX .NE. 0 )THEN CGRELAX.656
DO J=J_1,J_JMT ORH5F402.137
DO I=1,IMT CGRELAX.658
PTD(I,J)=2.0*PTD(I,J) CGRELAX.659
ENDDO CGRELAX.660
ENDDO CGRELAX.661
ENDIF CGRELAX.662
CGRELAX.663
CGRELAX.664
IF (L_OTIMER) CALL TIMER('CGRELAX ',4) CGRELAX.665
CGRELAX.666
RETURN CGRELAX.667
END CGRELAX.668
*ENDIF CGRELAX.669