Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

trosk.F90 @ 15548

Last change on this file since 15548 was 15548, checked in by gsamson, 3 years ago
update branch to the head of the trunk (r15547); ticket #2632
File size: 35.5 KB

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

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

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

Context Navigation

source: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/TOP/PISCES/SED/trosk.F90 @ 15548

Download in other formats: