Back to home page

MITgcm

 
 

    


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 CBOP
                0008 C     !ROUTINE: SOLVE_PENTADIAGONAL
                0009 C     !INTERFACE:
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 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
0d59e56009 Jean*0018 C     | S/R SOLVE_PENTADIAGONAL
527dc572d4 Jean*0019 C     | o Solve a penta-diagonal system A*X=Y (dimension Nr)
                0020 C     *==========================================================*
                0021 C     | o Used to solve implicitly vertical advection & diffusion
c3aa667b91 Jean*0022 #ifdef SOLVE_DIAGONAL_LOWMEMORY
                0023 C     | o Allows 2 types of call:
                0024 C     | 1) with INPUT errCode < 0 (first call):
                0025 C     |   Solve system A*X=Y by LU decomposition and return
                0026 C     |     inverse matrix as output (in a5d,b5d,c5d,d5d,e5d)
                0027 C     | 2) with INPUT errCode = 0 (subsequent calls):
                0028 C     |   Solve system A*Xp=Yp by using inverse matrix coeff
                0029 C     |     from first call output.
                0030 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
527dc572d4 Jean*0031 C     *==========================================================*
                0032 C     \ev
                0033 
                0034 C     !USES:
                0035       IMPLICIT NONE
                0036 C     == Global data ==
                0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 
                0040 C     !INPUT/OUTPUT PARAMETERS:
c3aa667b91 Jean*0041 C     -- INPUT: --
90bf640dcf Jean*0042 C     iMin,iMax,jMin,jMax  :: computational domain
c3aa667b91 Jean*0043 C     a5d     :: 2nd  lower diagonal of the pentadiagonal matrix
                0044 C     b5d     :: 1rst lower diagonal of the pentadiagonal matrix
                0045 C     c5d     :: main diagonal       of the pentadiagonal matrix
                0046 C     d5d     :: 1rst upper diagonal of the pentadiagonal matrix
                0047 C     e5d     :: 2nd  upper diagonal of the pentadiagonal matrix
                0048 C     y5d     :: Y vector (R.H.S.);
                0049 C     bi,bj   :: tile indices
                0050 C     myThid  :: thread number
                0051 C     -- OUTPUT: --
                0052 C     y5d     :: X = solution of A*X=Y
527dc572d4 Jean*0053 C     errCode :: > 0 if singular matrix
c3aa667b91 Jean*0054 #ifdef SOLVE_DIAGONAL_LOWMEMORY
                0055 C     a5d,b5d,c5d,d5d,e5d :: inverse matrix coeff to enable to find Xp solution
                0056 C                            of A*Xp=Yp with subsequent call to this routine
                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 C/*   This flag combination is only false if INCLUDE_IMPLVERTADV_CODE is
                0069 C     not defined while ALLOW_AUTODIFF is. */
                0070 #if ( defined INCLUDE_IMPLVERTADV_CODE || !defined ALLOW_AUTODIFF )
527dc572d4 Jean*0071 C     !LOCAL VARIABLES:
                0072 C     == Local variables ==
                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 C     kkey :: tape key, depends on i,j,k
                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 CEOP
                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 C--   forward sweep (starting from top) part-1: matrix L.U decomposition
                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 C--        [k] <- [k] - b_k/c_k-1 * [k-1]
                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 C--   Middle of forward sweep
                0133          DO j=jMin,jMax
                0134           DO i=iMin,iMax
                0135 C--        [k] <- [k] - a_k/c_k-2 * [k-2]
                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 C--        [k] <- [k] - b_k/c_k-1 * [k-1]
                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 C-      end if k= .. ; end of k loop
                0152         ENDIF
                0153        ENDDO
                0154 
                0155 C--   end if errCode < 0
                0156       ENDIF
                0157 
                0158       DO k=2,Nr
                0159 C--   forward sweep (starting from top) part-2: process RHS
                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 C-      end if k= .. ; end of k loop
                0174        ENDIF
527dc572d4 Jean*0175       ENDDO
                0176 
                0177 C--   Backward sweep (starting from bottom)
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 C-      end if k= .. ; end of k loop
                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 CADJ INIT loctape_solvepenta = COMMON, Nr*(sNy+2*OLy)*(sNx+2*OLx)
                0214 #endif
d666da8b50 Gael*0215 C--   Temporary array
                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 C--   Main loop
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 C--   forward sweep (starting from top)
                0238 
7448700841 Mart*0239 #ifdef ALLOW_AUTODIFF_TAMC
                0240           kkey = k + Nr*(i-1+OLx) + Nr*(sNx+2*OLx)*(j-1+OLy)
                0241 CADJ STORE d5d_prime(k) = loctape_solvepenta, key = kkey
                0242 CADJ STORE e5d_prime(k) = loctape_solvepenta, key = kkey
                0243 CADJ STORE y5d_prime(k) = loctape_solvepenta, key = kkey
                0244 #endif
d666da8b50 Gael*0245           IF (k.EQ.1) THEN
                0246 c just copy terms
                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 c subtract one term
                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 c subtract two terms
7448700841 Mart*0262 #ifdef ALLOW_AUTODIFF_TAMC
                0263 CADJ STORE d5d_prime(k) = loctape_solvepenta, key = kkey
                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 c normalization
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 C-- k loop
                0291         ENDDO
d666da8b50 Gael*0292 
                0293 C--   Backward sweep (starting from bottom)
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 C--   Update array
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 C-    end i,j loop
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 CADJ INIT loctape_solvepenta = COMMON, Nr
                0320 #endif
1dd6de01dd Jean*0321 C--   Init. + copy to temp. array
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 C--   Forward sweep (starting from top)
113b21cd45 Mart*0336 #ifdef ALLOW_AUTODIFF_TAMC
                0337 CADJ STORE d5d_prime(:,:,k) = loctape_solvepenta, key = k
                0338 CADJ STORE e5d_prime(:,:,k) = loctape_solvepenta, key = k
                0339 CADJ STORE y5d_prime(:,:,k) = loctape_solvepenta, key = k
                0340 #endif
e0fb0cb0b7 Mart*0341        IF ( k .EQ. 1 ) THEN
1dd6de01dd Jean*0342 C-    just copy terms
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 C-    subtract one term
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 C-    subtract two terms
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 C-    normalization
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 C-    end k-loop
d666da8b50 Gael*0396       ENDDO
                0397 
6b36f7748a Patr*0398 C--   Backward sweep (starting from bottom)
                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 C-    k-loop
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