New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trosk.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/trosk.F90 @ 15297

Last change on this file since 15297 was 15297, checked in by aumont, 3 years ago

major update of the sediment module

File size: 35.4 KB
Line 
1MODULE TROS
2!****************************************************************
3!* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY *
4!* DIFFERENTIAL EQUATIONS Y'=F(X,Y) BY ROSENBROCK METHOD.       *
5!* ------------------------------------------------------------ *
6!* ------------------------------------------------------------ *
7!* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77   *
8!*       [BIBLI 18].                                            *
9!*                                                              *
10!*                       F90 Release 1.0 By J-P Moreau, Paris   *
11!*                                (www.jpmoreau.fr)             *
12!****************************************************************
13  USE timing
14  USE in_out_manager, ONLY : ln_timing ! I/O manager
15  USE sed
16  USE sedfunc
17
18  IMPLICIT NONE
19  PRIVATE
20
21  PUBLIC ROSK
22
23  INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: NFCN, NJAC, NSTEP, NACCPT, NREJCT, NDEC, NSOL
24
25
26!define example #1
27  INTERFACE
28       SUBROUTINE JAC(NEQ,Y,DFY,LDFY,ACCMASK)
29         INTEGER, PARAMETER :: WP = KIND(1.0D0)
30         INTEGER, INTENT(IN) :: NEQ, LDFY
31         REAL(WP), DIMENSION(:,:,:), INTENT(OUT) :: DFY
32         INTEGER , DIMENSION(:), INTENT(IN) :: ACCMASK
33       END SUBROUTINE JAC
34  END INTERFACE
35
36
37  CONTAINS
38
39!**********************************************************************
40SUBROUTINE ROSK(ROSM,N,X,Y,XEND,H, RTOL,ATOL,ITOL,                  &
41           &    JAC, MLJAC, MUJAC, WORK,LWORK,IDID,ISTAT,RSTAT)
42! ---------------------------------------------------------------------
43!     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
44!     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
45!     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4
46!     (WITH STEP SIZE CONTROL).
47!     C.F. SECTION IV.7
48!
49!     AUTHORS: E. HAIRER AND G. WANNER
50!              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
51!              CH-1211 GENEVE 24, SWITZERLAND
52!              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
53!
54!     THIS CODE IS PART OF THE BOOK:
55!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
56!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
57!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
58!         SPRINGER-VERLAG (1990)
59!
60!     VERSION OF OCTOBER 12, 1990
61!
62!     INPUT PARAMETERS
63!     ----------------
64!     N           DIMENSION OF THE SYSTEM
65!
66!     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
67!                 VALUE OF F(X,Y):
68!                    SUBROUTINE FCN(N,X,Y,F)
69!                    REAL*8 X,Y(N),F(N)
70!                    F(1)=...   ETC.
71!
72!     X           INITIAL X-VALUE
73!
74!     Y(N)        INITIAL VALUES FOR Y
75!
76!     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
77!
78!     H           INITIAL STEP SIZE GUESS;
79!                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
80!                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
81!                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
82!                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
83!                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
84!                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
85!
86!     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
87!                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
88!
89!     ITOL        SWITCH FOR RTOL AND ATOL:
90!                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
91!                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
92!                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
93!                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
94!                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
95!                     RTOL(I)*ABS(Y(I))+ATOL(I).
96!
97!     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
98!                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
99!                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM:
100!                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
101!                    REAL*8 X,Y(N),DFY(LDFY,N)
102!                    DFY(1,1)= ...
103!                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
104!                 FURNISHED BY THE CALLING PROGRAM.
105!                 THE JACOBIAN IS TAKEN AS BANDED AND
106!                    THE PARTIAL DERIVATIVES ARE STORED
107!                    DIAGONAL-WISE AS
108!                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
109!
110!     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
111!                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
112!                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
113!                       THE MAIN DIAGONAL).
114!
115!     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
116!                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
117!                 NEED NOT BE DEFINED IF MLJAC=N.
118!
119!     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
120!                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
121!                 "LWORK" MUST BE AT LEAST
122!                             N*(LJAC+LE1+8)+5
123!                 WHERE
124!                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
125!                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
126!                 AND
127!                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
128!                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
129!
130!                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
131!                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
132!                 STORAGE REQUIREMENT IS
133!                             LWORK = 2*N*N+8*N+5.
134!
135!     LWORK       DECLARED LENGHT OF ARRAY "WORK".
136!
137! ----------------------------------------------------------------------
138!
139!     SOPHISTICATED SETTING OF PARAMETERS
140!     -----------------------------------
141!              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
142!              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5)
143!              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
144!              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
145!
146!    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
147!
148!    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
149!
150!    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
151!              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
152!                 WORK(3) <= HNEW/HOLD <= WORK(4)
153!              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
154!
155!    WORK(5)   AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS
156!              THE STEP SIZE IS MULTIPLIED BY WORK(5)
157!              DEFAULT VALUES: WORK(5)=0.1D0
158!
159!-----------------------------------------------------------------------
160!
161!     OUTPUT PARAMETERS
162!     -----------------
163!     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
164!                 (AFTER SUCCESSFUL RETURN X=XEND)
165!
166!     Y(N)        SOLUTION AT X
167!
168!     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
169!
170!     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
171!                   IDID=1  COMPUTATION SUCCESSFUL,
172!                   IDID=-1 COMPUTATION UNSUCCESSFUL.
173!
174! ---------------------------------------------------------
175! *** *** *** *** *** *** *** *** *** *** *** *** ***
176!          DECLARATIONS
177! *** *** *** *** *** *** *** *** *** *** *** *** ***
178      INTEGER, INTENT(in) :: ROSM, N, ITOL, MLJAC, MUJAC, LWORK
179      REAL(wp), DIMENSION(1), INTENT(in) :: ATOL, RTOL
180      INTEGER, INTENT(inout) :: IDID
181      INTEGER , DIMENSION(jpoce,3), INTENT(out) :: ISTAT
182      REAL(wp), DIMENSION(jpoce,2), INTENT(out) :: RSTAT
183
184      INTEGER :: NMAX, LDJAC, LDE, IEYNEW, IEDY1, IEDY, IEAK1
185      INTEGER :: IEAK2, IEAK3, IEAK4, IEFX, IEJAC, IEE
186      INTEGER :: ISTORE
187      REAL(wp) :: UROUND, HMAX, FAC1, FAC2, FACREJ, XEND, X
188      REAL(wp), DIMENSION(jpoce) :: H
189      REAL(wp), DIMENSION(jpoce, N) :: Y
190      REAL(wp), DIMENSION(LWORK) :: WORK
191      LOGICAL ARRET
192      EXTERNAL JAC
193! --------------------------------------------------------------------
194! --- COMMON STAT CAN BE USED FOR STATISTICS
195! ---    NFCN      NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
196!                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
197! ---    NJAC      NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
198!                  OR NUMERICALLY)
199! ---    NSTEP     NUMBER OF COMPUTED STEPS
200! ---    NACCPT    NUMBER OF ACCEPTED STEPS
201! ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
202!                  HAS BEEN ACCEPTED)
203! ---    NDEC      NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
204! ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
205! --------------------------------------------------------------------
206! *** *** *** *** *** *** ***
207!    SETTING THE PARAMETERS
208! *** *** *** *** *** *** ***
209
210      IF ( ln_timing ) CALL timing_start('rosk')
211
212      ALLOCATE (NFCN(jpoce), NJAC(jpoce), NSTEP(jpoce), NACCPT(jpoce), NREJCT(jpoce), NDEC(jpoce), NSOL(jpoce))
213
214      NFCN=0
215      NJAC=0
216      NSTEP=0
217      NACCPT=0
218      NREJCT=0
219      NDEC=0
220      NSOL=0
221      ARRET=.FALSE.
222! -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
223      NMAX = 100000
224! -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
225      IF(WORK(1) == 0.0)THEN
226         UROUND = 1.E-16
227      ELSE
228         UROUND = WORK(1)
229         IF(UROUND <= 1.E-14 .OR. UROUND >= 1.0)THEN
230            WRITE(NUMSED,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1)
231            ARRET=.TRUE.
232         END IF
233      END IF
234! -------- MAXIMAL STEP SIZE
235      IF(WORK(2) == 0.0)THEN
236         HMAX = XEND-X
237      ELSE
238         HMAX = WORK(2)
239      END IF
240! -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
241      IF(WORK(3) == 0.0)THEN
242         FAC1 = 5.0_wp
243      ELSE
244         FAC1 = 1.0/WORK(3)
245      END IF
246      IF(WORK(4) == 0.0)THEN
247         FAC2 = 1.0_wp / 6.0_wp
248      ELSE
249         FAC2 = 1.0_wp / WORK(4)
250      END IF
251! -------  FACREJ    FOR THE HUMP
252      IF(WORK(5) == 0.0)THEN
253         FACREJ = 0.1_wp
254      ELSE
255         FACREJ = WORK(5)
256      END IF
257! *** *** *** *** *** *** *** *** *** *** *** *** ***
258!         COMPUTATION OF ARRAY ENTRIES
259! *** *** *** *** *** *** *** *** *** *** *** *** ***
260! -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
261! -- JACOBIAN
262      LDJAC=MLJAC+MUJAC+1
263      LDE=2*MLJAC+MUJAC+1
264! ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
265      IEYNEW=6
266      IEDY1=IEYNEW+N
267      IEDY=IEDY1+N
268      IEAK1=IEDY+N
269      IEAK2=IEAK1+N
270      IEAK3=IEAK2+N
271      IEAK4=IEAK3+N
272      IEFX =IEAK4+N
273      IEJAC=IEFX +N
274      IEE  =IEJAC+N*LDJAC
275! ------ TOTAL STORAGE REQUIREMENT -----------
276      ISTORE=IEE+N*LDE-1
277      IF(ISTORE > LWORK)THEN
278         WRITE(NUMSED,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
279         ARRET=.TRUE.
280      END IF
281! ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
282      IF (ARRET) THEN
283         IDID=-1
284         RETURN
285      END IF
286
287! -------- CALL TO CORE INTEGRATOR ------------
288      IF (ROSM == 4) THEN
289         CALL RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,        &
290            MLJAC,MUJAC,IDID,                 &
291            NMAX,UROUND,FAC1,FAC2,FACREJ,     &
292            LDJAC,LDE,RSTAT )
293      ELSE IF (ROSM == 3) THEN
294         CALL RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,        &
295            MLJAC,MUJAC,IDID,                 &
296            NMAX,UROUND,FAC1,FAC2,FACREJ,     &
297            LDJAC,LDE,RSTAT )
298      ENDIF
299! ----------- RETURN -----------
300
301      ISTAT(:,1) = NFCN(:)
302      ISTAT(:,2) = NJAC(:)
303      ISTAT(:,3) = NSTEP(:)
304
305      DEALLOCATE (NFCN, NJAC, NSTEP, NACCPT, NREJCT, NDEC, NSOL )
306
307      IF ( ln_timing ) CALL timing_stop('rosk')
308
309      RETURN
310
311      END SUBROUTINE ROSK
312
313      SUBROUTINE RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,         &
314        MLJAC,MUJAC,IDID, NMAX,UROUND,FAC1,FAC2,FACREJ,       &
315        LFJAC,LE,RSTAT)
316! ----------------------------------------------------------
317!     CORE INTEGRATOR FOR ROS4
318!     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
319! ----------------------------------------------------------
320! ----------------------------------------------------------
321!         DECLARATIONS
322! ----------------------------------------------------------
323      IMPLICIT REAL(wp) (A-H,O-Z)
324      IMPLICIT INTEGER (I-N)
325
326      REAL(wp) :: ATOL(1),RTOL(1)
327      REAL(wp), DIMENSION(jpoce,N) :: Y, YNEW, DY1, DY, AK1, AK2, AK3, AK4, FX
328      REAL(wp), DIMENSION(jpoce,LFJAC,N) :: FJAC
329      REAL(wp), DIMENSION(jpoce, LE, N)  :: E
330      REAL(wp), DIMENSION(jpoce) :: H, HNEW, HMAXN, XI, FAC
331      REAL(wp), DIMENSION(jpoce) :: HC21, HC31, HC32, HC41, HC42, HC43
332      REAL(wp), DIMENSION(jpoce) :: XXOLD, HOPT, ERR
333      INTEGER, DIMENSION(jpoce,N) :: IP
334      LOGICAL, DIMENSION(jpoce) :: REJECT,RJECT2
335      INTEGER, DIMENSION(jpoce) :: ACCMASK, ENDMASK, ERRMASK
336      REAL(wp), DIMENSION(jpoce,2) :: RSTAT
337
338      IF ( ln_timing ) CALL timing_start('ro3cor')
339
340! ---- PREPARE BANDWIDTHS -----
341       MLE=MLJAC
342       MUE=MUJAC
343       MBJAC=MLJAC+MUJAC+1
344       MDIAG=MLE+MUE+1
345! *** *** *** *** *** *** ***
346!  INITIALISATIONS
347! *** *** *** *** *** *** ***
348      POSNEG=SIGN(1.0,XEND-X)
349      CALL RODAS3 (A21,A31,A32,A41,A42,A43,C21,C31,C32,C41,C42,C43,  &
350                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
351
352! --- INITIAL PREPARATIONS
353      DO ji = 1, jpoce
354         HMAXN(ji)=MIN(ABS(HMAX),ABS(XEND-X))
355         H(ji)=MIN(MAX(1.E-10,ABS(H(ji))),HMAXN(ji))
356         H(ji)=SIGN(H(ji),POSNEG)
357         REJECT(ji)=.FALSE.
358         XXOLD(ji)=X
359         XI(ji) = X
360      END DO
361      IRTRN = 1
362      ERRMASK(:) = 0
363      ENDMASK(:) = 0
364      ACCMASK(:) = ENDMASK(:)
365
366      IF (IRTRN < 0) GOTO 79
367! --- BASIC INTEGRATION STEP
368   1  CONTINUE
369      DO JI = 1, jpoce
370         IF (NSTEP(ji) > NMAX .OR. XI(ji)+0.1*H(ji) == XI(ji) .OR. ABS(H(ji)) <= UROUND) ERRMASK(ji) = 1
371         IF ((XI(ji)-XEND)*POSNEG+UROUND > 0.0) THEN
372             H(ji)=HOPT(ji)
373             ENDMASK(JI) = 1
374         END IF
375         IF ( ENDMASK(ji) == 0 ) THEN
376            HOPT(JI)=H(JI)
377            IF ((XI(ji)+H(ji)-XEND)*POSNEG > 0.0) H(ji)=XEND-XI(ji)
378         ENDIF
379      END DO
380
381      ACCMASK(:) = ENDMASK(:)
382
383      IF ( COUNT( ENDMASK(:) == 1 ) == jpoce ) THEN
384         IF ( ln_timing ) CALL timing_stop('ro3cor')
385         RETURN
386      ENDIF
387      IF ( COUNT( ERRMASK(:) == 1 ) > 0 ) GOTO 79
388
389      CALL sed_func(N,Y,DY1,ACCMASK)
390
391      NFCN(:)=NFCN(:)+1
392
393! *** *** *** *** *** *** ***
394!  COMPUTATION OF THE JACOBIAN
395! *** *** *** *** *** *** ***
396      NJAC(:)=NJAC(:)+1
397      CALL JAC(N,Y,FJAC,LFJAC,ACCMASK)
398   2  CONTINUE
399
400! *** *** *** *** *** *** ***
401!  COMPUTE THE STAGES
402! *** *** *** *** *** *** ***
403      DO ji = 1, jpoce
404         IF (ACCMASK(ji) == 0) THEN
405            NDEC(ji)=NDEC(ji)+1
406            HC21(ji)=C21/H(ji)
407            HC31(ji)=C31/H(ji)
408            HC32(ji)=C32/H(ji)
409            HC41(ji)=C41/H(ji)
410            HC42(ji)=C42/H(ji)
411            HC43(ji)=C43/H(ji)
412            FAC(ji)=1.0/(H(ji)*GAMMA)
413! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
414            DO 601 J=1,N
415            I1=MAX(1,MUJAC+2-J)
416            I2=MIN(MBJAC,N+MUJAC+1-J)
417            DO 600 I=I1,I2
418  600       E(ji,I+MLE,J)=-FJAC(ji,I,J)
419  601       E(ji,MDIAG,J)=E(ji,MDIAG,J)+FAC(ji)
420         ENDIF
421      END DO
422      CALL DECB(N,LE,E,MLE,MUE,IP,INFO,ACCMASK)
423
424! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
425! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
426! ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
427! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
428      DO ji = 1, jpoce
429         IF (ACCMASK(ji) == 0) THEN
430            AK1(ji,1:N)=DY1(ji,1:N)
431         ENDIF
432      END DO
433      CALL SOLB(N,LE,E,MLE,MUE,AK1,IP,ACCMASK)
434      DO ji = 1, jpoce
435         IF (ACCMASK(ji) == 0) THEN
436            AK2(ji,1:N)=DY1(ji,1:N)+HC21(ji)*AK1(ji,1:N)
437         ENDIF
438      END DO
439      CALL SOLB(N,LE,E,MLE,MUE,AK2,IP,ACCMASK)
440      DO ji = 1, jpoce
441         IF (ACCMASK(ji) == 0) THEN
442            YNEW(ji,1:N)=Y(ji,1:N)+A31*AK1(ji,1:N)+A32*AK2(ji,1:N)
443         ENDIF
444      END DO
445      CALL sed_func(N,YNEW,DY,ACCMASK)
446      DO ji = 1, jpoce
447         IF (ACCMASK(ji) == 0) THEN
448            AK3(ji,1:N)=DY(ji,1:N)+HC31(ji)*AK1(ji,1:N)+HC32(ji)*AK2(ji,1:N)
449         ENDIF
450      END DO
451      CALL SOLB(N,LE,E,MLE,MUE,AK3,IP,ACCMASK)
452      DO ji = 1, jpoce
453         IF (ACCMASK(ji) == 0) THEN
454            YNEW(ji,1:N)=Y(ji,1:N)+A41*AK1(ji,1:N)+A42*AK2(ji,1:N)+A43*AK3(ji,1:N)
455         ENDIF
456      END DO
457      CALL sed_func(N,YNEW,DY,ACCMASK)
458      DO ji = 1, jpoce
459         IF (ACCMASK(ji) == 0) THEN
460            DO I = 1, N
461               AK4(ji,I)=DY(ji,I)+HC41(ji)*AK1(ji,I)+HC42(ji)*AK2(ji,I)+HC43(ji)*AK3(ji,I)
462            END DO
463         ENDIF
464      END DO
465      CALL SOLB(N,LE,E,MLE,MUE,AK4,IP,ACCMASK)
466
467      DO ji = 1, jpoce
468         IF (ACCMASK(ji) == 0) THEN
469            NSOL(ji) = NSOL(ji)+4
470            NFCN(ji) = NFCN(ji)+2
471! *** *** *** *** *** *** ***
472!  ERROR ESTIMATION
473! *** *** *** *** *** *** ***
474            NSTEP(ji)=NSTEP(ji)+1
475! ------------ NEW SOLUTION ---------------
476            DO I = 1, N
477               YNEW(ji,I)=Y(ji,I)+B1*AK1(ji,I)+B2*AK2(ji,I)+B3*AK3(ji,I)+B4*AK4(ji,I)
478            END DO
479! ------------ COMPUTE ERROR ESTIMATION ----------------
480            ERR(JI) = 0.0_wp
481            DO I = 1, N
482               S = E1*AK1(ji,I)+E2*AK2(ji,I)+E3*AK3(ji,I)+E4*AK4(ji,I)
483               IF (ITOL == 0) THEN
484                  SK = ATOL(1)+RTOL(1)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I)))
485               ELSE
486                  SK = ATOL(I)+RTOL(I)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I)))
487               END IF
488               ERR(ji) = ERR(ji)+(S/SK)**2
489            END DO
490            ERR(ji) = SQRT(ERR(ji)/N)
491! --- COMPUTATION OF HNEW
492! --- WE REQUIRE .2<=HNEW/H<=6.
493            FAC(ji)  = MAX(FAC2,MIN(FAC1,(ERR(ji))**.333D0/.9D0))
494            HNEW(ji) = H(ji)/FAC(ji)
495! *** *** *** *** *** *** ***
496!  IS THE ERROR SMALL ENOUGH ?
497! *** *** *** *** *** *** ***
498            RJECT2(ji) = .TRUE.
499            IF (ERR(ji) <= 1.0) THEN
500! --- STEP IS ACCEPTED
501               NACCPT(ji) = NACCPT(ji)+1
502               Y(ji,1:N)  = YNEW(ji,1:N)
503               XXOLD(ji)  = XI(ji)
504               XI(ji) = XI(ji)+H(ji)
505               RSTAT(ji,2) = H(ji)
506               IF (IRTRN < 0) GOTO 79
507               IF (ABS(HNEW(ji)) > HMAXN(ji)) HNEW(ji)=POSNEG*HMAXN(ji)
508               IF (REJECT(ji)) HNEW(ji)=POSNEG*MIN(ABS(HNEW(ji)),ABS(H(ji)))
509               REJECT(ji) = .FALSE.
510               RJECT2(ji) = .FALSE.
511               IF (NACCPT(ji) == 1) RSTAT(ji,1) = (H(ji)+HNEW(ji))/2.0
512               H(ji) = HNEW(ji)
513               ACCMASK(ji) = 1
514            ELSE
515! --- STEP IS REJECTED
516               IF (RJECT2(ji)) HNEW(ji)   = H(ji)*FACREJ
517               IF (REJECT(ji)) RJECT2(ji) = .TRUE.
518               REJECT(ji) = .TRUE.
519               H(ji)=HNEW(ji)
520               IF (NACCPT(ji) >= 1) NREJCT(ji) = NREJCT(ji)+1
521               ACCMASK(ji) = 0
522            END IF
523         ENDIF
524      END DO
525      IF (COUNT( ACCMASK(:) == 0 ) > 0 ) GOTO 2
526      GOTO 1
527! --- EXIT
528 79   CONTINUE
529 979  FORMAT(' EXIT OF ROS3 AT X=',D16.7,'   H=',D16.7)
530      IDID=-1
531
532      IF ( ln_timing ) CALL timing_stop('ro3cor')
533
534      RETURN
535
536      END SUBROUTINE RO3COR
537
538      SUBROUTINE RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,         &
539        MLJAC,MUJAC,IDID, NMAX,UROUND,FAC1,FAC2,FACREJ,       &
540        LFJAC,LE,RSTAT)
541! ----------------------------------------------------------
542!     CORE INTEGRATOR FOR ROS4
543!     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
544! ----------------------------------------------------------
545! ----------------------------------------------------------
546!         DECLARATIONS
547! ----------------------------------------------------------
548      IMPLICIT REAL(wp) (A-H,O-Z)
549      IMPLICIT INTEGER (I-N)
550
551      REAL(wp) :: ATOL(1),RTOL(1)
552      REAL(wp), DIMENSION(jpoce,N) :: Y, YNEW, DY1, DY, AK1, AK2, AK3, AK4, FX
553      REAL(wp), DIMENSION(jpoce,N) :: AK5, AK6
554      REAL(wp), DIMENSION(jpoce,LFJAC,N) :: FJAC
555      REAL(wp), DIMENSION(jpoce, LE, N)  :: E
556      REAL(wp), DIMENSION(jpoce) :: H, HNEW, HMAXN, XI, FAC
557      REAL(wp), DIMENSION(jpoce) :: HC21, HC31, HC32, HC41, HC42, HC43
558      REAL(wp), DIMENSION(jpoce) :: HC51, HC52, HC53, HC54, HC61, HC62
559      REAL(wp), DIMENSION(jpoce) :: HC63, HC64, HC65
560      REAL(wp), DIMENSION(jpoce) :: XXOLD, HOPT, ERR
561      INTEGER, DIMENSION(jpoce,N) :: IP
562      LOGICAL, DIMENSION(jpoce) :: REJECT,RJECT2
563      INTEGER, DIMENSION(jpoce) :: ACCMASK, ENDMASK, ERRMASK
564      REAL(wp), DIMENSION(jpoce,2) :: RSTAT
565! ---- PREPARE BANDWIDTHS -----
566       MLE   = MLJAC
567       MUE   = MUJAC
568       MBJAC = MLJAC+MUJAC+1
569       MDIAG = MLE+MUE+1
570! *** *** *** *** *** *** ***
571!  INITIALISATIONS
572! *** *** *** *** *** *** ***
573      POSNEG = SIGN(1.0,XEND-X)
574      CALL RODAS4(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61,A62,A63,  &
575                A64,A65,C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,  &
576                C64,C65,B1,B2,B3,B4,B5,B6,E1,E2,E3,E4,E5,E6,GAMMA)
577
578! --- INITIAL PREPARATIONS
579      DO ji = 1, jpoce
580         HMAXN(ji) = MIN(ABS(HMAX),ABS(XEND-X))
581         H(ji) = MIN(MAX(1.E-10,ABS(H(ji))),HMAXN(ji))
582         H(ji) = SIGN(H(ji),POSNEG)
583         REJECT(ji) = .FALSE.
584         XXOLD(ji) = X
585         XI(ji) = X
586      END DO
587      IRTRN = 1
588      ERRMASK(:) = 0
589      ENDMASK(:) = 0
590      ACCMASK(:) = ENDMASK(:)
591
592      IF (IRTRN < 0) GOTO 79
593! --- BASIC INTEGRATION STEP
594   1  CONTINUE
595      DO JI = 1, jpoce
596         IF (NSTEP(ji) > NMAX .OR. XI(ji)+0.1*H(ji) == XI(ji) .OR. ABS(H(ji)) <= UROUND) THEN
597            ERRMASK(ji) = 1
598         ENDIF
599         IF ((XI(ji)-XEND)*POSNEG+UROUND > 0.0) THEN
600             H(ji) = HOPT(ji)
601             ENDMASK(JI) = 1
602         END IF
603         IF ( ENDMASK(ji) == 0 ) THEN
604            HOPT(JI) = H(JI)
605            IF ((XI(ji)+H(ji)-XEND)*POSNEG > 0.0) H(ji)=XEND-XI(ji)
606         ENDIF
607      END DO
608
609      ACCMASK(:) = ENDMASK(:)
610
611      IF ( COUNT( ENDMASK(:) == 1 ) == jpoce ) RETURN
612      IF ( COUNT( ERRMASK(:) == 1 ) > 0 ) GOTO 79
613
614      CALL sed_func(N,Y,DY1,ACCMASK)
615
616      NFCN(:) = NFCN(:) + 1
617! *** *** *** *** *** *** ***
618!  COMPUTATION OF THE JACOBIAN
619! *** *** *** *** *** *** ***
620      NJAC(:) = NJAC(:)+1
621      CALL JAC(N,Y,FJAC,LFJAC,jpoce,ACCMASK)
622   2  CONTINUE
623! *** *** *** *** *** *** ***
624!  COMPUTE THE STAGES
625! *** *** *** *** *** *** ***
626      DO ji = 1, jpoce
627         IF (ACCMASK(ji) == 0) THEN
628            NDEC(ji) = NDEC(ji)+1
629            HC21(ji) = C21/H(ji)
630            HC31(ji) = C31/H(ji)
631            HC32(ji) = C32/H(ji)
632            HC41(ji) = C41/H(ji)
633            HC42(ji) = C42/H(ji)
634            HC43(ji) = C43/H(ji)
635            HC51(ji) = C51/H(ji)
636            HC52(ji) = C52/H(ji)
637            HC53(ji) = C53/H(ji)
638            HC54(ji) = C54/H(ji)
639            HC61(ji) = C61/H(ji)
640            HC62(ji) = C62/H(ji)
641            HC63(ji) = C63/H(ji)
642            HC64(ji) = C64/H(ji)
643            HC65(ji) = C65/H(ji)
644
645            FAC(ji) = 1.0/(H(ji)*GAMMA)
646! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
647            DO 601 J=1,N
648            I1=MAX0(1,MUJAC+2-J)
649            I2=MIN0(MBJAC,N+MUJAC+1-J)
650            DO 600 I=I1,I2
651  600       E(ji,I+MLE,J)=-FJAC(ji,I,J)
652  601       E(ji,MDIAG,J)=E(ji,MDIAG,J)+FAC(ji)
653         ENDIF
654      END DO
655      CALL DECB(N,LE,E,MLE,MUE,IP,INFO,ACCMASK)
656
657! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
658! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
659! ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
660! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
661      DO ji = 1, jpoce
662         IF (ACCMASK(ji) == 0) THEN
663            AK1(ji,1:N) = DY1(ji,1:N)
664         ENDIF
665      END DO
666      CALL SOLB(N,LE,E,MLE,MUE,AK1,IP,ACCMASK)
667      DO ji = 1, jpoce
668         IF (ACCMASK(ji) == 0) THEN
669            YNEW(ji,1:N)=Y(ji,1:N)+A21*AK1(ji,1:N)
670         ENDIF
671      END DO
672      CALL sed_func(N,YNEW,DY,ACCMASK)
673      DO ji = 1, jpoce
674         IF (ACCMASK(ji) == 0) THEN
675            AK2(ji,1:N)=DY(ji,1:N)+HC21(ji)*AK1(ji,1:N)
676         ENDIF
677      END DO
678      CALL SOLB(N,LE,E,MLE,MUE,AK2,IP,ACCMASK)
679      DO ji = 1, jpoce
680         IF (ACCMASK(ji) == 0) THEN
681            YNEW(ji,1:N)=Y(ji,1:N)+A31*AK1(ji,1:N)+A32*AK2(ji,1:N)
682         ENDIF
683      END DO
684      CALL sed_func(N,YNEW,DY,ACCMASK)
685      DO ji = 1, jpoce
686         IF (ACCMASK(ji) == 0) THEN
687            AK3(ji,1:N)=DY(ji,1:N)+HC31(ji)*AK1(ji,1:N)+HC32(ji)*AK2(ji,1:N)
688         ENDIF
689      END DO
690      CALL SOLB(N,LE,E,MLE,MUE,AK3,IP,ACCMASK)
691      DO ji = 1, jpoce
692         IF (ACCMASK(ji) == 0) THEN
693            DO I = 1, N
694               YNEW(ji,I)=Y(ji,I)+A41*AK1(ji,I)+A42*AK2(ji,I)+A43*AK3(ji,I)
695            END DO
696         ENDIF
697      END DO
698      CALL sed_func(N,YNEW,DY,ACCMASK)
699      DO ji = 1, jpoce
700         IF (ACCMASK(ji) == 0) THEN
701            DO I = 1, N
702               AK4(ji,I)=DY(ji,I)+HC41(ji)*AK1(ji,I)+HC42(ji)*AK2(ji,I)+HC43(ji)*AK3(ji,I)
703            END DO
704         ENDIF
705      END DO
706      CALL SOLB(N,LE,E,MLE,MUE,AK4,IP,ACCMASK)
707      DO ji = 1, jpoce
708         IF (ACCMASK(ji) == 0) THEN
709            DO I = 1, N
710               YNEW(ji,I)=Y(ji,I)+A51*AK1(ji,I)+A52*AK2(ji,I)+A53*AK3(ji,I)+A54*AK4(ji,I)
711            END DO
712         ENDIF
713      END DO
714      CALL sed_func(N,YNEW,DY,ACCMASK)
715      DO ji = 1, jpoce
716         IF (ACCMASK(ji) == 0) THEN
717            DO I = 1, N
718               AK5(ji,I)=DY(ji,I)+HC51(ji)*AK1(ji,I)+HC52(ji)*AK2(ji,I)+HC53(ji)*AK3(ji,I)+HC54(ji)*AK4(ji,I)
719            END DO
720         ENDIF
721      END DO
722      CALL SOLB(N,LE,E,MLE,MUE,AK5,IP,ACCMASK)
723      DO ji = 1, jpoce
724         IF (ACCMASK(ji) == 0) THEN
725            DO I = 1, N
726               YNEW(ji,I)=Y(ji,I)+A61*AK1(ji,I)+A62*AK2(ji,I)+A63*AK3(ji,I)+A64*AK4(ji,I)+A65*AK5(ji,I)
727            END DO
728         ENDIF
729      END DO
730      CALL sed_func(N,YNEW,DY,ACCMASK)
731      DO ji = 1, jpoce
732         IF (ACCMASK(ji) == 0) THEN
733            DO I = 1, N
734               AK6(ji,I)=DY(ji,I)+HC61(ji)*AK1(ji,I)+HC62(ji)*AK2(ji,I)+HC63(ji)*AK3(ji,I)+HC64(ji)*AK4(ji,I)   &
735               &         + HC65(ji)*AK5(ji,I)
736            END DO
737         ENDIF
738      END DO
739      CALL SOLB(N,LE,E,MLE,MUE,AK6,IP,ACCMASK)
740
741      DO ji = 1, jpoce
742         IF (ACCMASK(ji) == 0) THEN
743            NSOL(ji) = NSOL(ji) + 6
744            NFCN(ji) = NFCN(ji) + 5
745! *** *** *** *** *** *** ***
746!  ERROR ESTIMATION
747! *** *** *** *** *** *** ***
748            NSTEP(ji) = NSTEP(ji)+1
749! ------------ NEW SOLUTION ---------------
750            DO 240 I = 1, N
751  240       YNEW(ji,I)=Y(ji,I)+B1*AK1(ji,I)+B2*AK2(ji,I)+B3*AK3(ji,I)+B4*AK4(ji,I)+B5*AK5(ji,I)+B6*AK6(ji,I)
752! ------------ COMPUTE ERROR ESTIMATION ----------------
753            ERR(JI) = 0.0_wp
754            DO 300 I = 1, N
755            S = E1*AK1(ji,I)+E2*AK2(ji,I)+E3*AK3(ji,I)+E4*AK4(ji,I)+E5*AK5(ji,I)+E6*AK6(ji,I)
756            IF (ITOL == 0) THEN
757               SK = ATOL(1)+RTOL(1)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I)))
758            ELSE
759               SK = ATOL(I)+RTOL(I)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I)))
760            END IF
761  300       ERR(ji) = ERR(ji)+(S/SK)**2
762            ERR(ji) = SQRT(ERR(ji)/N)
763! --- COMPUTATION OF HNEW
764! --- WE REQUIRE .2<=HNEW/H<=6.
765            FAC(ji)  = MAX(FAC2,MIN(FAC1,(ERR(ji))**0.25/0.9))
766            HNEW(ji) = H(ji)/FAC(ji)
767! *** *** *** *** *** *** ***
768!  IS THE ERROR SMALL ENOUGH ?
769! *** *** *** *** *** *** ***
770            RJECT2(ji) = .TRUE.
771            IF (ERR(ji) <= 1.0) THEN
772! --- STEP IS ACCEPTED
773               NACCPT(ji) = NACCPT(ji) + 1
774               Y(ji,1:N)  = YNEW(ji,1:N)
775               XXOLD(ji)  = XI(ji)
776               XI(ji) = XI(ji)+H(ji)
777               RSTAT(ji,2) = H(ji)
778               IF (IRTRN < 0) GOTO 79
779               IF (ABS(HNEW(ji)) > HMAXN(ji)) HNEW(ji)=POSNEG*HMAXN(ji)
780               IF (REJECT(ji)) HNEW(ji)=POSNEG*MIN(ABS(HNEW(ji)),ABS(H(ji)))
781               REJECT(ji) = .FALSE.
782               RJECT2(ji) = .FALSE.
783               IF (NACCPT(ji) == 1) RSTAT(ji,1) = (H(ji)+HNEW(ji))/2.0
784               H(ji) = HNEW(ji)
785               ACCMASK(ji) = 1
786            ELSE
787! --- STEP IS REJECTED
788               IF (RJECT2(ji)) HNEW(ji)=H(ji)*FACREJ
789               IF (REJECT(ji)) RJECT2(ji)=.TRUE.
790               REJECT(ji) = .TRUE.
791               H(ji) = HNEW(ji)
792               IF (NACCPT(ji) >= 1) NREJCT(ji)=NREJCT(ji)+1
793               ACCMASK(ji) = 0
794            END IF
795         ENDIF
796      END DO
797      IF (COUNT( ACCMASK(:) == 0 ) > 0 ) GOTO 2
798      GOTO 1
799! --- EXIT
800 79   CONTINUE
801 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7)
802      IDID=-1
803      RETURN
804
805      END SUBROUTINE RO4COR
806
807      SUBROUTINE RODAS3 (A21,A31,A32,A41,A42,A43,C21,C31,C32,C41,C42,C43,  &
808                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
809
810      REAL(wp), INTENT(out) :: A21, A31, A32, A41, A42, A43, C21, C31
811      REAL(wp), INTENT(out) :: C32, C41, C42, C43, B1, B2, B3, B4, E1
812      REAL(wp), INTENT(out) :: E2, E3, E4, GAMMA
813
814         A21= 0.0
815         A31= 2.0
816         A32= 0.0
817         A41= 2.0
818         A42= 0.0
819         A43= 1.0
820         C21= 4.0
821         C31= 1.0
822         C32=-1.0
823         C41= 1.0
824         C42=-1.0
825         C43=-8.0/3.0
826         B1= 2.0
827         B2= 0.0
828         B3= 1.0
829         B4= 1.0
830         E1= 0.0
831         E2= 0.0
832         E3= 0.0
833         E4= 1.0
834         GAMMA= 0.5
835      RETURN
836      END SUBROUTINE RODAS3
837
838      SUBROUTINE RODAS4(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61,A62,A63,  &
839                A64,A65,C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,  &
840                C64,C65,B1,B2,B3,B4,B5,B6,E1,E2,E3,E4,E5,E6,GAMMA)
841
842      REAL(wp), INTENT(out) :: A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61
843      REAL(wp), INTENT(out) :: A62,A63,A64,A65,C21,C31,C32,C41,C42,C43,C51
844      REAL(wp), INTENT(out) :: C52,C53,C54,C61,C62,C63,C64,C65,B1,B2,B3,B4,B5
845      REAL(wp), INTENT(out) :: B6,E1,E2,E3,E4,E5,E6,GAMMA
846
847      A21 = 0.1544000000000000E+01
848      A31 = 0.9466785280815826
849      A32 = 0.2557011698983284
850      A41 = 0.3314825187068521E+01
851      A42 = 0.2896124015972201E+01
852      A43 = 0.9986419139977817
853      A51 = 0.1221224509226641E+01
854      A52 = 0.6019134481288629E+01
855      A53 = 0.1253708332932087E+02
856      A54 =-0.6878860361058950
857      A61 = A51
858      A62 = A52
859      A63 = A53
860      A64 = A54
861      A65 = 1.0
862      C21 =-0.5668800000000000E+01
863      C31 =-0.2430093356833875E+01
864      C32 =-0.2063599157091915
865      C41 =-0.1073529058151375
866      C42 =-0.9594562251023355E+01
867      C43 =-0.2047028614809616E+02
868      C51 = 0.7496443313967647E+01
869      C52 =-0.1024680431464352E+02
870      C53 =-0.3399990352819905E+02
871      C54 = 0.1170890893206160E+02
872      C61 = 0.8083246795921522E+01
873      C62 =-0.7981132988064893E+01
874      C63 =-0.3152159432874371E+02
875      C64 = 0.1631930543123136E+02
876      C65 =-0.6058818238834054E+01
877      B1 = A51
878      B2 = A52
879      B3 = A53
880      B4 = A54
881      B5 = 1.0
882      B6 = 1.0
883      E1 = 0.0
884      E2 = 0.0
885      E3 = 0.0
886      E4 = 0.0
887      E5 = 0.0
888      E6 = 1.0
889      GAMMA= 0.25
890      RETURN
891      END SUBROUTINE RODAS4
892!
893      SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER, ACCMASK)
894      IMPLICIT INTEGER (I-N)
895      REAL(wp), DIMENSION(jpoce,NDIM,N) :: A
896      INTEGER, DIMENSION(jpoce, N) :: IP
897      INTEGER, DIMENSION(jpoce) :: ACCMASK
898      REAL(wp) :: T
899      INTEGER :: JI
900!-----------------------------------------------------------------------
901!  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
902!  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
903!  INPUT..
904!     N       ORDER OF THE ORIGINAL MATRIX A.
905!     NDIM    DECLARED DIMENSION OF ARRAY  A.
906!     A       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
907!                OF THE MATRIX ARE STORED IN THE COLUMNS OF  A  AND
908!                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
909!                ML+1 THROUGH 2*ML+MU+1 OF  A.
910!     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
911!     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
912!  OUTPUT..
913!     A       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
914!                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
915!     IP      INDEX VECTOR OF PIVOT INDICES.
916!     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
917!     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
918!                SINGULAR AT STAGE K.
919!  USE  SOLB  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
920!  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
921!  IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
922!
923!  REFERENCE..
924!     THIS IS A MODIFICATION OF
925!     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
926!     C.A.C.M. 15 (1972), P. 274.
927!-----------------------------------------------------------------------
928    IF ( ln_timing ) CALL timing_start('decb')
929
930    DO JI = 1, jpoce
931      IF (ACCMASK(ji) == 0) THEN
932      IER = 0
933      IP(ji,N) = 1
934      MD = ML + MU + 1
935      MD1 = MD + 1
936      JU = 0
937      IF (ML == 0) GO TO 70
938      IF (N == 1) GO TO 70
939      IF (N < MU+2) GO TO 7
940      DO 5 J = MU+2,N
941      DO 5 I = 1,ML
942  5   A(ji,I,J) = 0.0
943  7   NM1 = N - 1
944      DO 60 K = 1, NM1
945        KP1 = K + 1
946        M = MD
947        MDL = MIN(ML,N-K) + MD
948        DO 10 I = MD1, MDL
949          IF (ABS(A(ji,I,K)) > ABS(A(ji,M,K))) M = I
950 10     CONTINUE
951        IP(ji,K) = M + K - MD
952        T = A(ji,M,K)
953        IF (M == MD) GO TO 20
954        IP(ji,N) = -IP(ji,N)
955        A(ji,M,K) = A(ji,MD,K)
956        A(ji,MD,K) = T
957 20     CONTINUE
958        IF (T == 0.0) THEN
959           IER = K
960           IP(ji,N) = 0
961        ENDIF
962        T = 1.0/T
963        DO 30 I = MD1, MDL
964 30       A(ji,I,K) = -A(ji,I,K)*T
965        JU = MIN(MAX(JU,MU+IP(ji,K)),N)
966        MM = MD
967        IF (JU < KP1) GO TO 55
968        DO 50 J = KP1, JU
969          M = M - 1
970          MM = MM - 1
971          T = A(ji,M,J)
972          IF (M .EQ. MM) GO TO 35
973          A(ji,M,J) = A(ji,MM,J)
974          A(ji,MM,J) = T
975 35       CONTINUE
976          IF (T == 0.0) GO TO 45
977          JK = J - K
978          DO 40 I = MD1, MDL
979            IJK = I - JK
980 40         A(ji,IJK,J) = A(ji,IJK,J) + A(ji,I,K)*T
981 45       CONTINUE
982 50       CONTINUE
983 55     CONTINUE
984 60     CONTINUE
985 70   K = N
986      IF (A(ji,MD,N) == 0.0) THEN
987         IER = K
988         IP(ji,N) = 0
989      ENDIF
990      ENDIF
991      END DO
992
993      IF ( ln_timing ) CALL timing_stop('decb')
994
995      RETURN
996!----------------------- END OF SUBROUTINE DECB ------------------------
997      END SUBROUTINE DECB
998
999      SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP, ACCMASK)
1000      IMPLICIT INTEGER (I-N)
1001      REAL(wp) :: T
1002      REAL(wp), DIMENSION(jpoce,NDIM,N) :: A
1003      REAL(wp), DIMENSION(jpoce,N) :: B
1004      INTEGER, DIMENSION(jpoce,N)  :: IP 
1005      INTEGER :: IER
1006      INTEGER, DIMENSION(jpoce) :: ACCMASK
1007!-----------------------------------------------------------------------
1008!  SOLUTION OF LINEAR SYSTEM, A*X = B .
1009!  INPUT..
1010!    N      ORDER OF MATRIX A.
1011!    NDIM   DECLARED DIMENSION OF ARRAY  A .
1012!    A      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
1013!    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1014!    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1015!    B      RIGHT HAND SIDE VECTOR.
1016!    IP     PIVOT VECTOR OBTAINED FROM DECB.
1017!  DO NOT USE IF DECB HAS SET IER  <> 0.
1018!  OUTPUT..
1019!    B      SOLUTION VECTOR, X .
1020!-----------------------------------------------------------------------
1021
1022      IF ( ln_timing ) CALL timing_start('solb')
1023
1024      DO JI = 1, jpoce
1025      IF (ACCMASK(ji) == 0) THEN
1026      MD = ML + MU + 1
1027      MD1 = MD + 1
1028      MDM = MD - 1
1029      NM1 = N - 1
1030      IF (ML == 0) GO TO 25
1031      IF (N == 1) GO TO 50
1032      DO 20 K = 1, NM1
1033        M = IP(ji,K)
1034        T = B(ji,M)
1035        B(ji,M) = B(ji,K)
1036        B(ji,K) = T
1037        MDL = MIN(ML,N-K) + MD
1038        DO 10 I = MD1, MDL
1039          IMD = I + K - MD
1040 10       B(ji,IMD) = B(ji,IMD) + A(ji,I,K)*T
1041 20     CONTINUE
1042 25   CONTINUE
1043      DO 40 KB = 1, NM1
1044        K = N + 1 - KB
1045        B(ji,K) = B(ji,K)/A(ji,MD,K)
1046        T = -B(ji,K)
1047        KMD = MD - K
1048        LM = MAX(1,KMD+1)
1049        DO 30 I = LM, MDM
1050          IMD = I - KMD
1051 30       B(ji,IMD) = B(ji,IMD) + A(ji,I,K)*T
1052 40     CONTINUE
1053 50   B(ji,1) = B(ji,1)/A(ji,MD,1)
1054      ENDIF
1055      END DO
1056
1057      IF ( ln_timing ) CALL timing_stop('solb')
1058
1059      RETURN
1060!----------------------- END OF SUBROUTINE SOLB ------------------------
1061      END SUBROUTINE SOLB
1062
1063END MODULE TROS
Note: See TracBrowser for help on using the repository browser.