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

Context Navigation

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

Last change on this file since 15450 was 15450, checked in by cetlod, 8 months ago
Some updates to make the PISCES/SED module usable. Totally orthogonal with no effect on other parts of the code
File size: 35.5 KB

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

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

Download in other formats: