File indexing completed on 2023-09-03 05:09:54 UTC
view on githubraw file Latest commit 74487008 on 2023-09-03 01:50:18 UTC
113b21cd45 Mart*0001 #include "PACKAGES_CONFIG.h"
527dc572d4 Jean*0002 #include "CPP_OPTIONS.h"
113b21cd45 Mart*0003 #ifdef ALLOW_AUTODIFF
0004 # include "AUTODIFF_OPTIONS.h"
0005 #endif
527dc572d4 Jean*0006
0007
0008
0009
0d59e56009 Jean*0010 SUBROUTINE SOLVE_PENTADIAGONAL(
527dc572d4 Jean*0011 I iMin,iMax, jMin,jMax,
0012 U a5d, b5d, c5d, d5d, e5d,
0013 U y5d,
0014 O errCode,
0015 I bi, bj, myThid )
0016
0017
0d59e56009 Jean*0018
527dc572d4 Jean*0019
0020
0021
c3aa667b91 Jean*0022 #ifdef SOLVE_DIAGONAL_LOWMEMORY
0023
0024
0025
0026
0027
0028
0029
0030 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
527dc572d4 Jean*0031
0032
0033
0034
0035 IMPLICIT NONE
0036
0037 #include "SIZE.h"
0038 #include "EEPARAMS.h"
0039
0040
c3aa667b91 Jean*0041
90bf640dcf Jean*0042
c3aa667b91 Jean*0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
527dc572d4 Jean*0053
c3aa667b91 Jean*0054 #ifdef SOLVE_DIAGONAL_LOWMEMORY
0055
0056
0057 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
527dc572d4 Jean*0058 INTEGER iMin,iMax,jMin,jMax
0d59e56009 Jean*0059 _RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0060 _RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0061 _RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0062 _RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0063 _RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
23a7f3050f Jean*0064 _RL y5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
527dc572d4 Jean*0065 INTEGER errCode
0066 INTEGER bi, bj, myThid
0067
7448700841 Mart*0068
0069
0070 #if ( defined INCLUDE_IMPLVERTADV_CODE ||
527dc572d4 Jean*0071
0072
0073 INTEGER i,j,k
351f151dc9 Patr*0074 #ifndef SOLVE_DIAGONAL_LOWMEMORY
0d59e56009 Jean*0075 _RL tmpVar, recVar
0076 _RL y5d_m1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
6b36f7748a Patr*0077 # ifdef SOLVE_DIAGONAL_KINNER
0078 _RL c5d_prime(Nr)
0079 _RL d5d_prime(Nr)
0080 _RL e5d_prime(Nr)
0081 _RL y5d_prime(Nr)
0082 _RL y5d_update(Nr)
7448700841 Mart*0083 # ifdef ALLOW_AUTODIFF_TAMC
0084
0085 INTEGER kkey
0086 # endif
6b36f7748a Patr*0087 # else
0d59e56009 Jean*0088 _RL c5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0089 _RL d5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0090 _RL e5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0091 _RL y5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
33a04da737 Mart*0092 # endif
351f151dc9 Patr*0093 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
527dc572d4 Jean*0094
0095
351f151dc9 Patr*0096 #ifdef SOLVE_DIAGONAL_LOWMEMORY
33a04da737 Mart*0097
c3aa667b91 Jean*0098 IF ( errCode .LT. 0 ) THEN
0099 errCode = 0
0100
0101 DO k=1,Nr
0102
0103 IF (k.EQ.1) THEN
0104 DO j=jMin,jMax
0105 DO i=iMin,iMax
0106 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
0107 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
0108 ELSE
0109 c5d(i,j,k) = 0. _d 0
0110 errCode = 1
0111 ENDIF
0112 ENDDO
90bf640dcf Jean*0113 ENDDO
527dc572d4 Jean*0114
c3aa667b91 Jean*0115 ELSEIF (k.EQ.2) THEN
0116 DO j=jMin,jMax
0117 DO i=iMin,iMax
0118
0119 b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)
0120 c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)
0121 d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)
0122 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
0123 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
0124 ELSE
0125 c5d(i,j,k) = 0. _d 0
0126 errCode = 1
0127 ENDIF
0128 ENDDO
0129 ENDDO
0130
0131 ELSE
0132
0133 DO j=jMin,jMax
0134 DO i=iMin,iMax
0135
0136 a5d(i,j,k) = a5d(i,j,k)*c5d(i,j,k-2)
0137 b5d(i,j,k) = b5d(i,j,k) - a5d(i,j,k)*d5d(i,j,k-2)
0138 c5d(i,j,k) = c5d(i,j,k) - a5d(i,j,k)*e5d(i,j,k-2)
0139
0140 b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)
0141 c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)
0142 d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)
0143 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
0144 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
0145 ELSE
0146 c5d(i,j,k) = 0. _d 0
0147 errCode = 1
0148 ENDIF
0149 ENDDO
0150 ENDDO
0151
0152 ENDIF
0153 ENDDO
0154
0155
0156 ENDIF
0157
0158 DO k=2,Nr
0159
0160 IF (k.EQ.2) THEN
90bf640dcf Jean*0161 DO j=jMin,jMax
0162 DO i=iMin,iMax
23a7f3050f Jean*0163 y5d(i,j,k) = y5d(i,j,k) - b5d(i,j,k)*y5d(i,j,k-1)
90bf640dcf Jean*0164 ENDDO
0165 ENDDO
0166 ELSE
0167 DO j=jMin,jMax
0168 DO i=iMin,iMax
23a7f3050f Jean*0169 y5d(i,j,k) = y5d(i,j,k) - b5d(i,j,k)*y5d(i,j,k-1)
0170 & - a5d(i,j,k)*y5d(i,j,k-2)
90bf640dcf Jean*0171 ENDDO
527dc572d4 Jean*0172 ENDDO
90bf640dcf Jean*0173
0174 ENDIF
527dc572d4 Jean*0175 ENDDO
0176
0177
90bf640dcf Jean*0178 DO k=Nr,1,-1
0179 IF (k.EQ.Nr) THEN
0180 DO j=jMin,jMax
0181 DO i=iMin,iMax
23a7f3050f Jean*0182 y5d(i,j,k) = y5d(i,j,k)*c5d(i,j,k)
90bf640dcf Jean*0183 ENDDO
0184 ENDDO
0185 ELSEIF (k.EQ.Nr-1) THEN
0186 DO j=jMin,jMax
0187 DO i=iMin,iMax
23a7f3050f Jean*0188 y5d(i,j,k) = ( y5d(i,j,k)
0189 & - d5d(i,j,k)*y5d(i,j,k+1)
0190 & )*c5d(i,j,k)
90bf640dcf Jean*0191 ENDDO
0192 ENDDO
0193 ELSE
0194 DO j=jMin,jMax
0195 DO i=iMin,iMax
23a7f3050f Jean*0196 y5d(i,j,k) = ( y5d(i,j,k)
0197 & - d5d(i,j,k)*y5d(i,j,k+1)
0198 & - e5d(i,j,k)*y5d(i,j,k+2)
0199 & )*c5d(i,j,k)
90bf640dcf Jean*0200 ENDDO
527dc572d4 Jean*0201 ENDDO
90bf640dcf Jean*0202
0203 ENDIF
527dc572d4 Jean*0204 ENDDO
0205
351f151dc9 Patr*0206 #else /* ndef SOLVE_DIAGONAL_LOWMEMORY */
d666da8b50 Gael*0207
c3aa667b91 Jean*0208 errCode = 0
0209
1dd6de01dd Jean*0210 #ifdef SOLVE_DIAGONAL_KINNER
7448700841 Mart*0211
0212 #ifdef ALLOW_AUTODIFF_TAMC
0213
0214 #endif
d666da8b50 Gael*0215
0216 DO k=1,Nr
0d59e56009 Jean*0217 DO j=1-OLy,sNy+OLy
0218 DO i=1-OLx,sNx+OLx
23a7f3050f Jean*0219 y5d_m1(i,j,k) = y5d(i,j,k)
6b36f7748a Patr*0220 ENDDO
0221 ENDDO
d666da8b50 Gael*0222 ENDDO
6b36f7748a Patr*0223
d666da8b50 Gael*0224
0d59e56009 Jean*0225 DO j=1-OLy,sNy+OLy
0226 DO i=1-OLx,sNx+OLx
d666da8b50 Gael*0227
6b36f7748a Patr*0228 DO k=1,Nr
0229 c5d_prime(k) = 0. _d 0
0230 d5d_prime(k) = 0. _d 0
0231 e5d_prime(k) = 0. _d 0
0232 y5d_prime(k) = 0. _d 0
0233 y5d_update(k) = 0. _d 0
0234 ENDDO
d666da8b50 Gael*0235
6b36f7748a Patr*0236 DO k=1,Nr
d666da8b50 Gael*0237
0238
7448700841 Mart*0239 #ifdef ALLOW_AUTODIFF_TAMC
0240 kkey = k + Nr*(i-1+OLx) + Nr*(sNx+2*OLx)*(j-1+OLy)
0241
0242
0243
0244 #endif
d666da8b50 Gael*0245 IF (k.EQ.1) THEN
0246
0247 c5d_prime(k) = c5d(i,j,k)
0248 d5d_prime(k) = d5d(i,j,k)
0249 e5d_prime(k) = e5d(i,j,k)
0250 y5d_prime(k) = y5d_m1(i,j,k)
0251 ELSEIF (k.EQ.2) THEN
0252
0253 c5d_prime(k) = c5d(i,j,k)
0254 & -b5d(i,j,k)*d5d_prime(k-1)
0255 d5d_prime(k) = d5d(i,j,k)
0256 & -b5d(i,j,k)*e5d_prime(k-1)
0257 e5d_prime(k) = e5d(i,j,k)
0258 y5d_prime(k) = y5d_m1(i,j,k)
0259 & -b5d(i,j,k)*y5d_prime(k-1)
0260 ELSE
0261
7448700841 Mart*0262 #ifdef ALLOW_AUTODIFF_TAMC
0263
0264 #endif
d666da8b50 Gael*0265 c5d_prime(k) = c5d(i,j,k)
0266 & -a5d(i,j,k)*e5d_prime(k-2)
0267 & -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*d5d_prime(k-1)
0268 d5d_prime(k) = d5d(i,j,k)
0269 & -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*e5d_prime(k-1)
0270 e5d_prime(k) = e5d(i,j,k)
0271 y5d_prime(k) = y5d_m1(i,j,k)
0272 & -a5d(i,j,k)*y5d_prime(k-2)
0273 & -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*y5d_prime(k-1)
0274 ENDIF
0275
0276
0d59e56009 Jean*0277 tmpVar = c5d_prime(k)
0278 IF ( tmpVar.NE.0. _d 0 ) THEN
0279 recVar = 1. _d 0 / tmpVar
0280 d5d_prime(k) = d5d_prime(k)*recVar
0281 e5d_prime(k) = e5d_prime(k)*recVar
0282 y5d_prime(k) = y5d_prime(k)*recVar
d666da8b50 Gael*0283 ELSE
0284 d5d_prime(k) = 0. _d 0
0285 e5d_prime(k) = 0. _d 0
0286 y5d_prime(k) = 0. _d 0
0287 errCode = 1
0288 ENDIF
0289
6b36f7748a Patr*0290
0291 ENDDO
d666da8b50 Gael*0292
0293
6b36f7748a Patr*0294 DO k=Nr,1,-1
0295 IF (k.EQ.Nr) THEN
d666da8b50 Gael*0296 y5d_update(k) = y5d_prime(k)
6b36f7748a Patr*0297 ELSEIF (k.EQ.Nr-1) THEN
d666da8b50 Gael*0298 y5d_update(k) = y5d_prime(k)
0299 & - y5d_update(k+1)*d5d_prime(k)
6b36f7748a Patr*0300 ELSE
d666da8b50 Gael*0301 y5d_update(k) = y5d_prime(k)
0302 & - y5d_update(k+1)*d5d_prime(k)
0303 & - y5d_update(k+2)*e5d_prime(k)
6b36f7748a Patr*0304 ENDIF
0305 ENDDO
d666da8b50 Gael*0306
0307
6b36f7748a Patr*0308 DO k=1,Nr
23a7f3050f Jean*0309 y5d(i,j,k) = y5d_update(k)
6b36f7748a Patr*0310 ENDDO
0311
1dd6de01dd Jean*0312
6b36f7748a Patr*0313 ENDDO
d666da8b50 Gael*0314 ENDDO
0315
6b36f7748a Patr*0316 #else /* ndef SOLVE_DIAGONAL_KINNER */
0317
113b21cd45 Mart*0318 #ifdef ALLOW_AUTODIFF_TAMC
0319
0320 #endif
1dd6de01dd Jean*0321
6b36f7748a Patr*0322 DO k=1,Nr
0d59e56009 Jean*0323 DO j=1-OLy,sNy+OLy
0324 DO i=1-OLx,sNx+OLx
6b36f7748a Patr*0325 c5d_prime(i,j,k) = 0. _d 0
0326 d5d_prime(i,j,k) = 0. _d 0
0327 e5d_prime(i,j,k) = 0. _d 0
0328 y5d_prime(i,j,k) = 0. _d 0
1dd6de01dd Jean*0329 y5d_m1(i,j,k) = y5d(i,j,k)
6b36f7748a Patr*0330 ENDDO
0331 ENDDO
d666da8b50 Gael*0332 ENDDO
6b36f7748a Patr*0333
0334 DO k=1,Nr
1dd6de01dd Jean*0335
113b21cd45 Mart*0336 #ifdef ALLOW_AUTODIFF_TAMC
0337
0338
0339
0340 #endif
e0fb0cb0b7 Mart*0341 IF ( k .EQ. 1 ) THEN
1dd6de01dd Jean*0342
e0fb0cb0b7 Mart*0343 DO j=1-OLy,sNy+OLy
0344 DO i=1-OLx,sNx+OLx
0345 c5d_prime(i,j,k) = c5d(i,j,k)
0346 d5d_prime(i,j,k) = d5d(i,j,k)
0347 e5d_prime(i,j,k) = e5d(i,j,k)
0348 y5d_prime(i,j,k) = y5d_m1(i,j,k)
0349 ENDDO
0350 ENDDO
0351 ELSEIF ( k .EQ. 2 ) THEN
1dd6de01dd Jean*0352
e0fb0cb0b7 Mart*0353 DO j=1-OLy,sNy+OLy
0354 DO i=1-OLx,sNx+OLx
0355 tmpVar = b5d(i,j,k)
0356 c5d_prime(i,j,k) = c5d(i,j,k) - tmpVar*d5d_prime(i,j,k-1)
0357 d5d_prime(i,j,k) = d5d(i,j,k) - tmpVar*e5d_prime(i,j,k-1)
0358 e5d_prime(i,j,k) = e5d(i,j,k)
0359 y5d_prime(i,j,k) = y5d_m1(i,j,k) - tmpVar*y5d_prime(i,j,k-1)
0360 ENDDO
0361 ENDDO
0362 ELSE
1dd6de01dd Jean*0363
e0fb0cb0b7 Mart*0364 DO j=1-OLy,sNy+OLy
0365 DO i=1-OLx,sNx+OLx
0366 tmpVar = b5d(i,j,k)-a5d(i,j,k)*d5d_prime(i,j,k-2)
0367 c5d_prime(i,j,k) = c5d(i,j,k) - tmpVar*d5d_prime(i,j,k-1)
0368 & -a5d(i,j,k)*e5d_prime(i,j,k-2)
0369 d5d_prime(i,j,k) = d5d(i,j,k) - tmpVar*e5d_prime(i,j,k-1)
0370 e5d_prime(i,j,k) = e5d(i,j,k)
0371 y5d_prime(i,j,k) = y5d_m1(i,j,k) - tmpVar*y5d_prime(i,j,k-1)
0372 & -a5d(i,j,k)*y5d_prime(i,j,k-2)
0373 ENDDO
0374 ENDDO
0375 ENDIF
6b36f7748a Patr*0376
1dd6de01dd Jean*0377
e0fb0cb0b7 Mart*0378 DO j=1-OLy,sNy+OLy
0379 DO i=1-OLx,sNx+OLx
0380 tmpVar = c5d_prime(i,j,k)
0381 IF ( tmpVar.NE.0. _d 0 ) THEN
0382 recVar = 1. _d 0 / tmpVar
0383 d5d_prime(i,j,k) = d5d_prime(i,j,k)*recVar
0384 e5d_prime(i,j,k) = e5d_prime(i,j,k)*recVar
0385 y5d_prime(i,j,k) = y5d_prime(i,j,k)*recVar
0386 ELSE
0387 d5d_prime(i,j,k) = 0. _d 0
0388 e5d_prime(i,j,k) = 0. _d 0
0389 y5d_prime(i,j,k) = 0. _d 0
0390 errCode = 1
0391 ENDIF
6b36f7748a Patr*0392 ENDDO
0393 ENDDO
0394
1dd6de01dd Jean*0395
d666da8b50 Gael*0396 ENDDO
0397
6b36f7748a Patr*0398
0399 DO k=Nr,1,-1
e0fb0cb0b7 Mart*0400 IF ( k .EQ. Nr ) THEN
0401 DO j=1-OLy,sNy+OLy
0402 DO i=1-OLx,sNx+OLx
1dd6de01dd Jean*0403 y5d(i,j,k) = y5d_prime(i,j,k)
e0fb0cb0b7 Mart*0404 ENDDO
0405 ENDDO
0406 ELSEIF ( k .EQ. Nr-1 ) THEN
0407 DO j=1-OLy,sNy+OLy
0408 DO i=1-OLx,sNx+OLx
1dd6de01dd Jean*0409 y5d(i,j,k) = y5d_prime(i,j,k)
0410 & - y5d(i,j,k+1)*d5d_prime(i,j,k)
e0fb0cb0b7 Mart*0411 ENDDO
0412 ENDDO
0413 ELSE
0414 DO j=1-OLy,sNy+OLy
0415 DO i=1-OLx,sNx+OLx
1dd6de01dd Jean*0416 y5d(i,j,k) = y5d_prime(i,j,k)
0417 & - y5d(i,j,k+1)*d5d_prime(i,j,k)
0418 & - y5d(i,j,k+2)*e5d_prime(i,j,k)
e0fb0cb0b7 Mart*0419 ENDDO
6b36f7748a Patr*0420 ENDDO
e0fb0cb0b7 Mart*0421 ENDIF
0422
6b36f7748a Patr*0423 ENDDO
0424
0425 #endif /* SOLVE_DIAGONAL_KINNER */
0426
351f151dc9 Patr*0427 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
d666da8b50 Gael*0428
7448700841 Mart*0429 #endif /* INCLUDE_IMPLVERTADV_CODE */
0430
527dc572d4 Jean*0431 RETURN
0432 END