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

Context Navigation

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, 2 years ago
major update of the sediment module
File size: 35.4 KB

Line
1	MODULE 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	!**********************************************************************
40	SUBROUTINE 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 = 2NN+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)+A31AK1(ji,1:N)+A32AK2(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)+A41AK1(ji,1:N)+A42AK2(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)+B1AK1(ji,I)+B2AK2(ji,I)+B3AK3(ji,I)+B4AK4(ji,I)
478	END DO
479	! ------------ COMPUTE ERROR ESTIMATION ----------------
480	ERR(JI) = 0.0_wp
481	DO I = 1, N
482	S = E1AK1(ji,I)+E2AK2(ji,I)+E3AK3(ji,I)+E4AK4(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)+A31AK1(ji,1:N)+A32AK2(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)+A41AK1(ji,I)+A42AK2(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)+A51AK1(ji,I)+A52AK2(ji,I)+A53AK3(ji,I)+A54AK4(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)+A61AK1(ji,I)+A62AK2(ji,I)+A63AK3(ji,I)+A64AK4(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)+B1AK1(ji,I)+B2AK2(ji,I)+B3AK3(ji,I)+B4AK4(ji,I)+B5AK5(ji,I)+B6AK6(ji,I)
752	! ------------ COMPUTE ERROR ESTIMATION ----------------
753	ERR(JI) = 0.0_wp
754	DO 300 I = 1, N
755	S = E1AK1(ji,I)+E2AK2(ji,I)+E3AK3(ji,I)+E4AK4(ji,I)+E5AK5(ji,I)+E6AK6(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
1063	END MODULE TROS

Note: See TracBrowser for help on using the repository browser.

Download in other formats: