forked from KineticPreProcessor/KPP-Standalone
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gckpp_Integrator.F90
2706 lines (2382 loc) · 92.7 KB
/
gckpp_Integrator.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Numerical Integrator (Time-Stepping) File
!
! Generated by KPP-3.1.0 symbolic chemistry Kinetics PreProcessor
! (https:/github.com/KineticPreProcessor/KPP
! KPP is distributed under GPL, the general public licence
! (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2022, A. Sandu, Michigan Tech, Virginia Tech
! With important contributions from:
! M. Damian, Villanova University, Philadelphia, PA, USA
! R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! M. Long, Renaissance Fiber, LLC, North Carolina, USA
! H. Lin, Harvard University, Cambridge, MA, USA
! R. Yantosca, Harvard University, Cambridge, MA, USA
!
! File : gckpp_Integrator.F90
! Equation file : gckpp.kpp
! Output root filename : gckpp
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! INTEGRATE - Integrator routine
! Arguments :
! TIN - Start Time for Integration
! TOUT - End Time for Integration
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! Rosenbrock - Implementation of several Rosenbrock methods: !
! * Ros2 !
! * Ros3 !
! * Ros4 !
! * Rodas3 !
! * Rodas4 !
! By default the code employs the KPP sparse linear algebra routines !
! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
! !
! (C) Adrian Sandu, August 2004 !
! Virginia Polytechnic Institute and State University !
! Contact: [email protected] !
! Revised by Philipp Miehe and Adrian Sandu, May 2006 !
! !
! Revised by Mike Long and Haipeng Lin to add auto-reduce fun. !
! Harvard University, Atmospheric Chemistry Modeling Group !
! Contact: [email protected] April 2022 !
! !
! This implementation is part of KPP - the Kinetic PreProcessor !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
MODULE gckpp_Integrator
USE gckpp_Parameters
USE gckpp_Global
IMPLICIT NONE
PUBLIC
SAVE
!~~~> Flags to determine if we should call the UPDATE_* routines from within
!~~~> the integrator. If using KPP in an external model, you might want to
!~~~> disable these calls (via ICNTRL(15)) to avoid excess computations.
LOGICAL, PRIVATE :: Do_Update_RCONST
LOGICAL, PRIVATE :: Do_Update_PHOTO
LOGICAL, PRIVATE :: Do_Update_SUN
!~~~> Statistics on the work performed by the Rosenbrock method
INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
Ntexit=1, Nhexit=2, Nhnew = 3, &
NARthr=4
CONTAINS
SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, &
ISTATUS_U, RSTATUS_U, IERR_U )
USE gckpp_Util, ONLY : Integrator_Update_Options
IMPLICIT NONE
REAL(kind=dp), INTENT(IN) :: TIN ! Start Time
REAL(kind=dp), INTENT(IN) :: TOUT ! End Time
!~~~> Optional input parameters and statistics
INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
REAL(kind=dp), INTENT(IN), OPTIONAL :: RCNTRL_U(20)
INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
REAL(kind=dp), INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
INTEGER, INTENT(OUT), OPTIONAL :: IERR_U
REAL(kind=dp) :: RCNTRL(20), RSTATUS(20)
INTEGER :: ICNTRL(20), ISTATUS(20), IERR
INTEGER, SAVE :: Ntotal = 0
!~~~> Zero input and output arrays for safety's sake
ICNTRL = 0
RCNTRL = 0.0_dp
ISTATUS = 0
RSTATUS = 0.0_dp
!~~~> fine-tune the integrator:
ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int.
!~~~> if optional parameters are given, and if they are /= 0,
! then use them to overwrite default settings
IF ( PRESENT( ICNTRL_U ) ) THEN
WHERE( ICNTRL_U /= 0 ) ICNTRL = ICNTRL_U
ENDIF
IF ( PRESENT( RCNTRL_U ) ) THEN
WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U
ENDIF
!~~~> Determine the settings of the Do_Update_* flags, which determine
!~~~> whether or not we need to call Update_* routines in the integrator
!~~~> (or not, if we are calling them from a higher-level)
! ICNTRL(15) = -1 ! Do not call Update_* functions within the integrator
! = 0 ! Status quo
! = 1 ! Call Update_RCONST from within the integrator
! = 2 ! Call Update_PHOTO from within the integrator
! = 3 ! Call Update_RCONST and Update_PHOTO from w/in the int.
! = 4 ! Call Update_SUN from within the integrator
! = 5 ! Call Update_SUN and Update_RCONST from within the int.
! = 6 ! Call Update_SUN and Update_PHOTO from within the int.
! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int.
CALL Integrator_Update_Options( ICNTRL(15), &
Do_Update_RCONST, &
Do_Update_PHOTO, &
Do_Update_Sun )
!~~~> In order to remove the prior EQUIVALENCE statements (which
!~~~> are not thread-safe), we now have declared VAR and FIX as
!~~~> threadprivate pointer variables that can point to C.
VAR => C(1:NVAR )
FIX => C(NVAR+1:NSPEC)
!~~~> Call the integrator
CALL Rosenbrock( NVAR, VAR, TIN, TOUT, ATOL, RTOL, &
RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR )
!~~~> Free pointers
VAR => NULL()
FIX => NULL()
!~~~> Debug option: show number of steps
!Ntotal = Ntotal + ISTATUS(Nstp)
!PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=', VAR(ind_O3)
STEPMIN = RSTATUS(Nhexit)
!~~~> if optional parameters are given for output
!~~~> use them to store information in them
IF ( PRESENT( ISTATUS_U ) ) ISTATUS_U = ISTATUS
IF ( PRESENT( RSTATUS_U ) ) RSTATUS_U = RSTATUS
IF ( PRESENT( IERR_U ) ) IERR_U = IERR
END SUBROUTINE INTEGRATE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
AbsTol,RelTol, &
RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
!
! G = 1/(H*gamma(1)) - Jac(t0,Y0)
! T_i = t0 + Alpha(i)*H
! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
! gamma(i)*dF/dT(t0, Y0)
! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
!
! For details on Rosenbrock methods and their implementation consult:
! E. Hairer and G. Wanner
! "Solving ODEs II. Stiff and differential-algebraic problems".
! Springer series in computational mathematics, Springer-Verlag, 1996.
! The codes contained in the book inspired this implementation.
!
! (C) Adrian Sandu, August 2004
! Virginia Polytechnic Institute and State University
! Contact: [email protected]
! Revised by Philipp Miehe and Adrian Sandu, May 2006
! This implementation is part of KPP - the Kinetic PreProcessor
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~> INPUT ARGUMENTS:
!
!- Y(N) = vector of initial conditions (at T=Tstart)
!- [Tstart,Tend] = time range of integration
! (if Tstart>Tend the integration is performed backwards in time)
!- RelTol, AbsTol = user precribed accuracy
!- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
! returns Ydot = Y' = F(T,Y)
!- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
! returns Jcb = dFun/dY
!- ICNTRL(1:20) = integer inputs parameters
!- RCNTRL(1:20) = real inputs parameters
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~> OUTPUT ARGUMENTS:
!
!- Y(N) -> vector of final states (at T->Tend)
!- ISTATUS(1:20) -> integer output parameters
!- RSTATUS(1:20) -> real output parameters
!- IERR -> job status upon return
! success (positive value) or
! failure (negative value)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~> INPUT PARAMETERS:
!
! Note: For input parameters equal to zero the default values of the
! corresponding variables are used.
!
! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
!
! ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
! = 1: AbsTol, RelTol are scalars
!
! ICNTRL(3) -> selection of a particular Rosenbrock method
! = 0 : Rodas3 (default)
! = 1 : Ros2
! = 2 : Ros3
! = 3 : Ros4
! = 4 : Rodas3
! = 5 : Rodas4
!
! ICNTRL(4) -> maximum number of integration steps
! For ICNTRL(4)=0) the default value of 200000 is used
!
! ICNTRL(12) -> use auto-reduce solver? set threshold in RCNTRL(12)
! ICNTRL(13) -> ... append slow species when auto-reducing?
! ICNTRL(14) -> choose a target species instead for determining threshold?
! if yes, specify idx. then RCNTRL(12) is obsolete.
!
! ICNTRL(15) -> Toggles calling of Update_* functions w/in the integrator
! = -1 : Do not call Update_* functions within the integrator
! = 0 : Status quo
! = 1 : Call Update_RCONST from within the integrator
! = 2 : Call Update_PHOTO from within the integrator
! = 3 : Call Update_RCONST and Update_PHOTO from w/in the int.
! = 4 : Call Update_SUN from within the integrator
! = 5 : Call Update_SUN and Update_RCONST from within the int.
! = 6 : Call Update_SUN and Update_PHOTO from within the int.
! = 7 : Call Update_SUN, Update_PHOTO, Update_RCONST w/in the int.
!
! ICNTRL(16) ->
! = 0 : allow negative concentrations (default)
! = 1 : set negative concentrations to zero
!
! RCNTRL(1) -> Hmin, lower bound for the integration step size
! It is strongly recommended to keep Hmin = ZERO
! RCNTRL(2) -> Hmax, upper bound for the integration step size
! RCNTRL(3) -> Hstart, starting value for the integration step size
!
! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
! (default=0.1)
! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
! than the predicted value (default=0.9)
!
! RCNTRL(12) -> threshold for auto-reduction (req. ICNTRL(12)) (default=100)
! RCNTRL(14) -> AR threshold ratio (default=0.01)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
! OUTPUT ARGUMENTS:
! -----------------
!
! T -> T value for which the solution has been computed
! (after successful return T=Tend).
!
! Y(N) -> Numerical solution at T
!
! IDID -> Reports on successfulness upon return:
! = 1 for success
! < 0 for error (value equals error code)
!
! ISTATUS(1) -> No. of function calls
! ISTATUS(2) -> No. of jacobian calls
! ISTATUS(3) -> No. of steps
! ISTATUS(4) -> No. of accepted steps
! ISTATUS(5) -> No. of rejected steps (except at very beginning)
! ISTATUS(6) -> No. of LU decompositions
! ISTATUS(7) -> No. of forward/backward substitutions
! ISTATUS(8) -> No. of singular matrix decompositions
!
! RSTATUS(1) -> Texit, the time corresponding to the
! computed Y upon return
! RSTATUS(2) -> Hexit, last accepted step before exit
! RSTATUS(3) -> Hnew, last predicted step (not yet taken)
! For multiple restarts, use Hnew as Hstart
! in the subsequent run
! RSTATUS(4) -> ARthr, last auto-reduction threshold determined
! only if AR is on (ICNTRL(12)) and key spc (ICNTRL(14))
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
USE gckpp_Parameters
USE gckpp_LinearAlgebra
IMPLICIT NONE
!~~~> Arguments
INTEGER, INTENT(IN) :: N
REAL(kind=dp), INTENT(INOUT) :: Y(N)
REAL(kind=dp), INTENT(IN) :: Tstart,Tend
REAL(kind=dp), INTENT(IN) :: AbsTol(N),RelTol(N)
INTEGER, INTENT(IN) :: ICNTRL(20)
REAL(kind=dp), INTENT(IN) :: RCNTRL(20)
INTEGER, INTENT(INOUT) :: ISTATUS(20)
REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
INTEGER, INTENT(OUT) :: IERR
!~~~> Parameters of the Rosenbrock method, up to 6 stages
INTEGER :: ros_S, rosMethod
INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5, RG3=6
REAL(kind=dp) :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
ros_Alpha(6), ros_Gamma(6), ros_ELO
LOGICAL :: ros_NewF(6)
CHARACTER(LEN=12) :: ros_Name
!~~~> Local variables
REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
REAL(kind=dp) :: Hmin, Hmax, Hstart
REAL(kind=dp) :: Texit, Redux_Threshold
INTEGER :: i, UplimTol, Max_no_steps
LOGICAL :: Autonomous, VectorTol, Autoreduce, Autoreduce_Append
INTEGER :: AR_target_spc
REAL(kind=dp) :: AR_thr_ratio
!~~~> Parameters
REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
!~~~> Initialize statistics
ISTATUS(1:8) = 0
RSTATUS(1:4) = ZERO
!~~~> Autonomous or time dependent ODE. Default is time dependent.
Autonomous = .NOT.(ICNTRL(1) == 0)
!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
IF (ICNTRL(2) == 0) THEN
VectorTol = .TRUE.
UplimTol = N
ELSE
VectorTol = .FALSE.
UplimTol = 1
END IF
!~~~> Initialize the particular Rosenbrock method selected
SELECT CASE (ICNTRL(3))
CASE (1)
CALL Ros2
CASE (2)
CALL Ros3
CASE (3)
CALL Ros4
CASE (0,4)
CALL Rodas3
CASE (5)
CALL Rodas4
CASE (6)
CALL Rang3
CASE DEFAULT
PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3)
CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
RETURN
END SELECT
!~~~> The maximum number of steps admitted
IF (ICNTRL(4) == 0) THEN
Max_no_steps = 200000
ELSEIF (ICNTRL(4) > 0) THEN
Max_no_steps=ICNTRL(4)
ELSE
PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> Auto-reduction toggle
Autoreduce = .false.
IF (ICNTRL(12) == 1) Autoreduce = .true.
Autoreduce_Append = ICNTRL(13) == 1
!~~~> Target species (if zero, uses the regular threshold)
AR_target_spc = ICNTRL(14)
!~~~> Unit roundoff (1+Roundoff>1)
Roundoff = WLAMCH('E')
!~~~> Lower bound on the step size: (positive value)
IF (RCNTRL(1) == ZERO) THEN
Hmin = ZERO
ELSEIF (RCNTRL(1) > ZERO) THEN
Hmin = RCNTRL(1)
ELSE
PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> Upper bound on the step size: (positive value)
IF (RCNTRL(2) == ZERO) THEN
Hmax = ABS(Tend-Tstart)
ELSEIF (RCNTRL(2) > ZERO) THEN
Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
ELSE
PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> Starting step size: (positive value)
IF (RCNTRL(3) == ZERO) THEN
Hstart = MAX(Hmin,DeltaMin)
ELSEIF (RCNTRL(3) > ZERO) THEN
Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
ELSE
PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> Step size can be changed s.t. FacMin < Hnew/Hold < FacMax
IF (RCNTRL(4) == ZERO) THEN
FacMin = 0.2_dp
ELSEIF (RCNTRL(4) > ZERO) THEN
FacMin = RCNTRL(4)
ELSE
PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
RETURN
END IF
IF (RCNTRL(5) == ZERO) THEN
FacMax = 6.0_dp
ELSEIF (RCNTRL(5) > ZERO) THEN
FacMax = RCNTRL(5)
ELSE
PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> FacRej: Factor to decrease step after 2 succesive rejections
IF (RCNTRL(6) == ZERO) THEN
FacRej = 0.1_dp
ELSEIF (RCNTRL(6) > ZERO) THEN
FacRej = RCNTRL(6)
ELSE
PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> FacSafe: Safety Factor in the computation of new step size
IF (RCNTRL(7) == ZERO) THEN
FacSafe = 0.9_dp
ELSEIF (RCNTRL(7) > ZERO) THEN
FacSafe = RCNTRL(7)
ELSE
PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
RETURN
END IF
!~~~> Check if tolerances are reasonable
DO i=1,UplimTol
IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
.OR. (RelTol(i) >= 1.0_dp) ) THEN
PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
PRINT * , ' RelTol(',i,') = ',RelTol(i)
CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
RETURN
END IF
END DO
!~~~> Auto-reduction threshold
Redux_threshold = 1.d2
IF (RCNTRL(12) > ZERO) THEN
Redux_Threshold = RCNTRL(12)
ELSEIF (RCNTRL(12) < ZERO) THEN
Autoreduce = .false.
ENDIF
!~~~> Auto-reduction threshold ratio (only if ICNTRL(14) is not zero)
AR_thr_ratio = RCNTRL(14)
!~~~> CALL Auto-reducing Rosenbrock method
IF ( Autoreduce .and. .not. Autoreduce_Append ) THEN
! ros_yIntegrator is the aggressively micro-optimized revision by Haipeng Lin.
! ros_cIntegrator is the original auto-reduce implementation by Mike Long.
CALL ros_yIntegrator(Y, Tstart, Tend, Texit, &
AbsTol, RelTol, &
! Integration parameters
Autonomous, VectorTol, Max_no_steps, &
Roundoff, Hmin, Hmax, Hstart, &
FacMin, FacMax, FacRej, FacSafe, &
! Autoreduce threshold
redux_threshold, AR_target_spc, &
AR_thr_ratio, &
! Error indicator
IERR)
ENDIF
!~~~> CALL Auto-reducing Rosenbrock method capable of append
! (this version is less aggressively optimized.)
IF ( Autoreduce .and. Autoreduce_Append ) THEN
! ros_yIntegratorA is the append version of the AR integrator. It has less
! optimizations because of the need to update Prod/Loss.
CALL ros_yIntegratorA(Y, Tstart, Tend, Texit, &
AbsTol, RelTol, &
! Integration parameters
Autonomous, VectorTol, Max_no_steps, &
Roundoff, Hmin, Hmax, Hstart, &
FacMin, FacMax, FacRej, FacSafe, &
! Autoreduce threshold
redux_threshold, AR_target_spc, &
AR_thr_ratio, &
! Error indicator
IERR)
ENDIF
!~~~> CALL Normal Rosenbrock method
IF ( .not. Autoreduce .or. IERR .eq. -99 ) &
CALL ros_Integrator(Y, Tstart, Tend, Texit, &
AbsTol, RelTol, &
! Integration parameters
Autonomous, VectorTol, Max_no_steps, &
Roundoff, Hmin, Hmax, Hstart, &
FacMin, FacMax, FacRej, FacSafe, &
! Error indicator
IERR)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS ! SUBROUTINES internal to Rosenbrock
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Handles all error messages
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
REAL(kind=dp), INTENT(IN) :: T, H
INTEGER, INTENT(IN) :: Code
INTEGER, INTENT(OUT) :: IERR
IERR = Code
PRINT * , &
'Forced exit from Rosenbrock due to the following error:'
SELECT CASE (Code)
CASE (-1)
PRINT * , '--> Improper value for maximal no of steps'
CASE (-2)
PRINT * , '--> Selected Rosenbrock method not implemented'
CASE (-3)
PRINT * , '--> Hmin/Hmax/Hstart must be positive'
CASE (-4)
PRINT * , '--> FacMin/FacMax/FacRej must be positive'
CASE (-5)
PRINT * , '--> Improper tolerance values'
CASE (-6)
PRINT * , '--> No of steps exceeds maximum bound'
CASE (-7)
PRINT * , '--> Step size too small: T + 10*H = T', &
' or H < Roundoff'
CASE (-8)
PRINT * , '--> Matrix is repeatedly singular'
CASE DEFAULT
PRINT *, 'Unknown Error code: ', Code
END SELECT
PRINT *, "T=", T, "and H=", H
END SUBROUTINE ros_ErrorMsg
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, &
AbsTol, RelTol, &
!~~~> Integration parameters
Autonomous, VectorTol, Max_no_steps, &
Roundoff, Hmin, Hmax, Hstart, &
FacMin, FacMax, FacRej, FacSafe, &
!~~~> Error indicator
IERR )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Template for the implementation of a generic Rosenbrock method
! defined by ros_S (no of stages)
! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
IMPLICIT NONE
!~~~> Input: the initial condition at Tstart; Output: the solution at T
REAL(kind=dp), INTENT(INOUT) :: Y(N)
!~~~> Input: integration interval
REAL(kind=dp), INTENT(IN) :: Tstart,Tend
!~~~> Output: time at which the solution is returned (T=Tend if success)
REAL(kind=dp), INTENT(OUT) :: T
!~~~> Input: tolerances
REAL(kind=dp), INTENT(IN) :: AbsTol(N), RelTol(N)
!~~~> Input: integration parameters
LOGICAL, INTENT(IN) :: Autonomous, VectorTol
REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
INTEGER, INTENT(IN) :: Max_no_steps
REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
!~~~> Output: Error indicator
INTEGER, INTENT(OUT) :: IERR
! ~~~~ Local variables
REAL(kind=dp) :: Ynew(N), Fcn0(N), Fcn(N)
REAL(kind=dp) :: K(N*ros_S), dFdT(N)
#ifdef FULL_ALGEBRA
REAL(kind=dp) :: Jac0(N,N), Ghimj(N,N)
#else
REAL(kind=dp) :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
#endif
REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau
REAL(kind=dp) :: Err, Yerr(N)
INTEGER :: Pivot(N), Direction, ioffset, j, istage
LOGICAL :: RejectLastH, RejectMoreH, Singular
!~~~> Local parameters
REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
!~~~> Locally called functions
! REAL(kind=dp) WLAMCH
! EXTERNAL WLAMCH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~> Initial preparations
DO_SLV = .true.
DO_FUN = .true.
DO_JVS = .true.
T = Tstart
RSTATUS(Nhexit) = ZERO
H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
IF (Tend >= Tstart) THEN
Direction = +1
ELSE
Direction = -1
END IF
H = Direction*H
RejectLastH=.FALSE.
RejectMoreH=.FALSE.
!~~~> Time loop begins below
TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
CALL ros_ErrorMsg(-6,T,H,IERR)
RETURN
END IF
IF ( ((T+0.1_dp*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
CALL ros_ErrorMsg(-7,T,H,IERR)
RETURN
END IF
!~~~> Limit H if necessary to avoid going beyond Tend
H = MIN(H,ABS(Tend-T))
!~~~> Compute the function at current time
CALL FunTemplate( T, Y, Fcn0 )
ISTATUS(Nfun) = ISTATUS(Nfun) + 1
!~~~> Compute the function derivative with respect to T
IF (.NOT.Autonomous) THEN
CALL ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT )
END IF
!~~~> Compute the Jacobian at current time
CALL JacTemplate( T, Y, Jac0 )
ISTATUS(Njac) = ISTATUS(Njac) + 1
!~~~> Repeat step calculation until current step accepted
UntilAccepted: DO
CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
Jac0,Ghimj,Pivot,Singular)
IF (Singular) THEN ! More than 5 consecutive failed decompositions
CALL ros_ErrorMsg(-8,T,H,IERR)
RETURN
END IF
!~~~> Compute the stages
Stage: DO istage = 1, ros_S
! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
ioffset = N*(istage-1)
! For the 1st istage the function has been computed previously
IF ( istage == 1 ) THEN
!slim: CALL WCOPY(N,Fcn0,1,Fcn,1)
Fcn(1:N) = Fcn0(1:N)
! istage>1 and a new function evaluation is needed at the current istage
ELSEIF ( ros_NewF(istage) ) THEN
!slim: CALL WCOPY(N,Y,1,Ynew,1)
Ynew(1:N) = Y(1:N)
DO j = 1, istage-1
CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), &
K(N*(j-1)+1),1,Ynew,1)
END DO
Tau = T + ros_Alpha(istage)*Direction*H
CALL FunTemplate( Tau, Ynew, Fcn )
ISTATUS(Nfun) = ISTATUS(Nfun) + 1
END IF ! if istage == 1 elseif ros_NewF(istage)
!slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1)
K(ioffset+1:ioffset+N) = Fcn(1:N)
DO j = 1, istage-1
HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1)
END DO
IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
HG = Direction*H*ros_Gamma(istage)
CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1)
END IF
CALL ros_Solve(Ghimj, Pivot, K(ioffset+1))
END DO Stage
!~~~> Compute the new solution
!slim: CALL WCOPY(N,Y,1,Ynew,1)
Ynew(1:N) = Y(1:N)
DO j=1,ros_S
CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1)
END DO
!~~~> Compute the error estimation
!slim: CALL WSCAL(N,ZERO,Yerr,1)
Yerr(1:N) = ZERO
DO j=1,ros_S
CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1)
END DO
Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
Hnew = H*Fac
!~~~> Check the error magnitude and adjust step size
ISTATUS(Nstp) = ISTATUS(Nstp) + 1
IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
ISTATUS(Nacc) = ISTATUS(Nacc) + 1
IF (ICNTRL(16) == 1) THEN
! new value is non-negative:
Y = MAX(Ynew,ZERO)
ELSE
!slim: CALL WCOPY(N,Ynew,1,Y,1)
Y(1:N) = Ynew(1:N)
ENDIF
T = T + Direction*H
Hnew = MAX(Hmin,MIN(Hnew,Hmax))
IF (RejectLastH) THEN ! No step size increase after a rejected step
Hnew = MIN(Hnew,H)
END IF
RSTATUS(Nhexit) = H
RSTATUS(Nhnew) = Hnew
RSTATUS(Ntexit) = T
RejectLastH = .FALSE.
RejectMoreH = .FALSE.
H = Hnew
EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
ELSE !~~~> Reject step
IF (RejectMoreH) THEN
Hnew = H*FacRej
END IF
RejectMoreH = RejectLastH
RejectLastH = .TRUE.
H = Hnew
IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
END IF ! Err <= 1
END DO UntilAccepted
END DO TimeLoop
!~~~> Succesful exit
IERR = 1 !~~~> The integration was successful
END SUBROUTINE ros_Integrator
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, &
AbsTol, RelTol, &
!~~~> Integration parameters
Autonomous, VectorTol, Max_no_steps, &
Roundoff, Hmin, Hmax, Hstart, &
FacMin, FacMax, FacRej, FacSafe, &
!~~~> Autoreduce threshold
threshold, AR_target_spc, AR_thr_ratio, &
!~~~> Error indicator
IERR )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Template for the implementation of a generic Rosenbrock method
! defined by ros_S (no of stages)
! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Alternative micro-optimized implementation, hplin, 4/10/22
! which does not resize arrays to rNVAR, instead always keeping to full N size
! and skipping using DO_SLV
!
! All compiled assembly code was verified just short of linking to a proper BLAS.
USE gckpp_Global, ONLY : cNONZERO, rNVAR
use gckpp_Monitor, ONLY : SPC_NAMES
USE gckpp_JacobianSP
IMPLICIT NONE
!~~~> Input: the initial condition at Tstart; Output: the solution at T
REAL(kind=dp), INTENT(INOUT) :: Y(N)
!~~~> Input: integration interval
REAL(kind=dp), INTENT(IN) :: Tstart,Tend
!~~~> Output: time at which the solution is returned (T=Tend if success)
REAL(kind=dp), INTENT(OUT) :: T
!~~~> Input: tolerances
REAL(kind=dp), INTENT(IN) :: AbsTol(N), RelTol(N)
!~~~> Input: integration parameters
LOGICAL, INTENT(IN) :: Autonomous, VectorTol
REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
INTEGER, INTENT(IN) :: Max_no_steps
REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
!~~~> Autoreduction threshold
REAL(kind=dp), INTENT(IN) :: threshold
!~~~> Output: Error indicator
INTEGER, INTENT(OUT) :: IERR
! ~~~~ Local variables
REAL(kind=dp) :: Ynew(N), Fcn0(N), Fcn(N), Prod(N), Loss(N), LossY(N)
REAL(kind=dp) :: K(NVAR*ros_S), dFdT(N)
#ifdef FULL_ALGEBRA
REAL(kind=dp) :: Jac0(N,N), Ghimj(N,N)
#else
REAL(kind=dp) :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
REAL(kind=dp) :: cGhimj(LU_NONZERO) ! not known at this point what cNONZERO will be
#endif
REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau
REAL(kind=dp) :: Err, Yerr(N), Yerrsub(NVAR)
INTEGER :: Pivot(N), Direction, ioffset, j, istage
LOGICAL :: RejectLastH, RejectMoreH, Singular, Reduced
!~~~> Local parameters
REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
INTEGER :: SPC
REAL(kind=dp) :: alpha_factor ! hplin 4/10/22
! Inline local parameters for AR.
INTEGER :: II, III, idx, nrmv, s
REAL(kind=dp), INTENT(IN) :: AR_thr_ratio
REAL(kind=dp) :: AR_thr
INTEGER, INTENT(IN) :: AR_target_spc
!~~~> Initial preparations
DO_SLV = .true.
DO_FUN = .true.
DO_JVS = .true.
Reduced = .false.
T = Tstart
RSTATUS(Nhexit) = ZERO
H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
IF (Tend >= Tstart) THEN
Direction = +1
ELSE
Direction = -1
END IF
H = Direction*H
RejectLastH=.FALSE.
RejectMoreH=.FALSE.
! reset K - hplin. this is a multiplier that gets applied to ros_A and ros_M
K = 0.0_dp
Ghimj = 0.0_dp
!~~~> Time loop begins below
TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
CALL ros_ErrorMsg(-6,T,H,IERR)
RETURN
ENDIF
IF ( ((T+0.1_dp*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
CALL ros_ErrorMsg(-7,T,H,IERR)
RETURN
ENDIF
!~~~> Limit H if necessary to avoid going beyond Tend
H = MIN(H,ABS(Tend-T))
! ... 0.40% ... 5.2 %
!~~~> Compute the function at current time
! this is necessary anyway for the rest of the time loop. do not optimize out.
IF (T .eq. Tstart) THEN
CALL FunSplitF(T,Y,Fcn0,Prod,Loss,LossY) ! always calculates PL.
ELSE ! ELSE or ne seems reasonably close in performance.
CALL FunSplitN(T,Y,Fcn0)
ENDIF
! The above Prod, Loss were updated at every TimeLoop, so they no longer
! reflect initial condition. The copy operation is extremely expensive,
! so we instead run it once at the beginning of the timeloop and use a IF
! to switch it to a wholy different codepath...
ISTATUS(Nfun) = ISTATUS(Nfun) + 1
! ... 1.02% ... 13.00%
!~~~> Parse species for reduced computation
if (.not. reduced) then
! Inline the entire reduction operation here.
iSPC_MAP = 0
NRMV = 0
S = 1
AR_thr = threshold
! Target species?
if(AR_target_spc .gt. 0) then
AR_thr = AR_thr_ratio * max(LossY(AR_target_spc), Prod(AR_target_spc)) ! Lin et al., 2022 in prep.
RSTATUS(NARthr) = AR_thr
endif
! Checks should be kept out of tight inner loops.
IF(keepActive) THEN
DO i=1,NVAR
! Short-circuiting using SKIP is very important here.
if (.not. keepSpcActive(i) .and. &
abs(LossY(i)).lt.AR_thr .and. abs(Prod(i)).lt.AR_thr) then ! per Shen et al., 2020
NRMV=NRMV+1
! RMV(NRMV) = i ! not needed unless in append version.
DO_SLV(i) = .false.
! DO_FUN(i) = .false.
cycle
endif
SPC_MAP(S) = i ! Add to full spc map.
iSPC_MAP(i) = S
S=S+1
ENDDO
ENDIF
IF (.not. keepActive) THEN
DO i=1,NVAR
! Short-circuiting using SKIP is very important here.
if (abs(LossY(i)).lt.AR_thr .and. abs(Prod(i)).lt.AR_thr) then ! per Shen et al., 2020
NRMV=NRMV+1
! RMV(NRMV) = i ! not needed unless in append version.
DO_SLV(i) = .false.
! DO_FUN(i) = .false.
cycle
endif
SPC_MAP(S) = i ! Add to full spc map.
iSPC_MAP(i) = S
S=S+1
ENDDO
ENDIF
rNVAR = NVAR-NRMV ! Number of active species in the reduced mechanism
II = 1
III = 1
idx = 0
! This loop can be unrolled into two branches.
ScanFirstNonZero: DO i = 1,LU_NONZERO
IF ((DO_SLV(LU_IROW(i))).and.(DO_SLV(LU_ICOL(i)))) THEN
idx = 1
cLU_IROW(1) = iSPC_MAP(LU_IROW(i))
cLU_ICOL(1) = iSPC_MAP(LU_ICOL(i))
JVS_MAP(1) = i
EXIT ScanFirstNonZero
ENDIF
DO_JVS(i) = .false.
ENDDO ScanFirstNonZero
DO i = i+1,LU_NONZERO ! There is no escape for looping through LU_NONZERO here.
IF ((DO_SLV(LU_IROW(i))).and.(DO_SLV(LU_ICOL(i)))) THEN
idx=idx+1 ! counter for the number of non-zero elements in the reduced Jacobian
cLU_IROW(idx) = iSPC_MAP(LU_IROW(i))
cLU_ICOL(idx) = iSPC_MAP(LU_ICOL(i))
JVS_MAP(idx) = i
IF (cLU_IROW(idx).ne.cLU_IROW(idx-1)) THEN
II=II+1
cLU_CROW(II) = idx
ENDIF
IF (cLU_IROW(idx).eq.cLU_ICOL(idx)) THEN
III=III+1
cLU_DIAG(III) = idx
ENDIF
CYCLE
ENDIF
DO_JVS(i) = .false.
ENDDO
cNONZERO = idx
cLU_CROW(1) = 1 ! 1st index = 1