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 @ 15545

Last change on this file since 15545 was 15545, checked in by acc, 7 months ago

Add missing declaration to src/TOP/PISCES/SED/trosk.F90 which prevents compilation with debug settings

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