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/trunk/src/TOP/PISCES/SED – NEMO

source: NEMO/trunk/src/TOP/PISCES/SED/trosk.F90 @ 15450

Last change on this file since 15450 was 15450, checked in by cetlod, 8 months ago

Some updates to make the PISCES/SED module usable. Totally orthogonal with no effect on other parts of the code

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