-
Notifications
You must be signed in to change notification settings - Fork 0
/
input.f90
2862 lines (2534 loc) · 107 KB
/
input.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
!*********************************************************************************************************
module input
!This module contains input related subroutines. The species and reactions included in the perturbation are
!controlled by the spec_pert and rxn_pert vectors.
!*********************************************************************************************************
implicit none
contains
!*********************************************************************************************************
subroutine inco_tube
! This subroutine gets the program options, etc. contained in tube.inp. It is
! read once at the beginning of the program and never used again as the values
! read here are constant.
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles; use globalSA
implicit none
integer :: i, ierr
logical :: liso
double precision :: press
double precision :: logdt, logtout, omega_tmp
character(len=80) :: rxn_str_tmp
integer :: mem_stat, mem_stat_total=0
integer :: nlines
integer :: RxnID, NT, NVAL
double precision :: RVAL(1)
logical :: kerr
character(len=80) :: MARI_str, conv_str
integer :: SpecID
integer :: maxl
!-------------------------------------------------------------------------------------------------------
!======================================================================================================
! Section 1 -- Digest tube.inp
!======================================================================================================
!.........Open input file:
open(unit=itubeinp, file='INP.d/tube.inp', status='old', action='read')
!-------------------------------------------------------------------------------------------------------
!Reactor type, number of runs, and whether each run has unique parameters
call skip_comment(itubeinp)
read(itubeinp,*) irxtr, nruns, MultiInput
!Read operating conditions:
if (MultiInput) then
call skip_comment(itubeinp)
read(itubeinp,*) lstp
else
call skip_comment(itubeinp)
read(itubeinp,*) lstp,tin,press,velo,abyv,trise
endif
!More code options: gas/surface chemistry, isothermality, and reset reactor conditions
!to be consistent with the chosen options.
call skip_comment(itubeinp)
read(itubeinp,*)liso,itpd
if (liso) then
iiso=0 !Zero here will make the non-isothermal term drop out of the energy balance
else
iiso=1
end if
if (itpd==0) then
trpa_called=.true.
else
if (itpd==1) irxtr=0 !Must use a /UHV/ reactor for a UHV TPD
ltra=.true. !TPD experiments are inherently transient -- transient behavior is desired
end if
if (irxtr<=1) velo=0.0D0 !Batch and UHV reactors should have no flow in/out
!Heat transfer
call skip_comment(itubeinp)
read(itubeinp,*)text,aextbyv,htc,ramp
!Specify which species should be (1) printed during output and (2) is the primary reactant
!Use name/phase/ construction to simplify entry.
call skip_comment(itubeinp)
read(itubeinp,*) MARI_str, conv_str
!Convert species names into Species ID's
call sksnum(MARI_str,0,LOUT,knams,kmax,PhaseNames,nphases,SpecTot,SpecID,NT,NVAL,RVAL(1),KERR)
MARI_index=SpecID
call sksnum(conv_str,0,LOUT,knams,kmax,PhaseNames,nphases,SpecTot,SpecID,NT,NVAL,RVAL(1),KERR)
conv_index=SpecID
if (MARI_index==0 .or. conv_index==0) then
write(*,*) 'Invalid species specified for the MARI or the reactant. Stopping.'
STOP
end if
!Options for controlling the number of reactors, the time steps, and how much
!output is desired
call skip_comment(itubeinp)
read(itubeinp,*) rlen, nnodes, ttout, rtime, ntdec, ltra
if (irxtr/=3) nnodes=1 !If not a PFR, only 1 node
if (irxtr>1) then
iflow=1 !CSTR and PFR are flow reactors
else
iflow=0 !Batch/UHV are not flow reactors
end if
!Solver parameters
call skip_comment(itubeinp)
read(itubeinp,*) ltol,abstol,reltol,NonNeg, idid_restart_count_max
call skip_comment(itubeinp)
read(itubeinp,*) iSolver, mu, ml
!Read the option flags for semi-empirical techniques
call skip_comment(itubeinp)
read(itubeinp,*) lcov,lStatpQ,lBEP,iScale,lEA,lomega,omega_default, Tref_beta
!Reaction path and sensitivity analysis
call skip_comment(itubeinp)
read(itubeinp,*)mrpa,verbose_rpa,trpa,lsenbrute,lDOE
close(itubeinp)
!======================================================================================================
! Section 2 -- Digest omega.inp if needed (or just set to default)
!======================================================================================================
if (lomega) then !Omega correction on, use default plus values from file
omega=omega_default
open(unit=iomegainp, file='INP.d/omega.inp', status='unknown',action='read') !proximity factor
call skip_comment(iomegainp)
read(iomegainp,*)nlines
call skip_comment(iomegainp)
do i=1,nlines
read(iomegainp,*)omega_tmp,rxn_str_tmp
call skcomp(rxn_str_tmp,rxn_str,nrxn,RxnID,NT)
if (RxnID==0) then !Reaction not found
write(*,*) 'Invalid reaction specified in omega.inp on line ', i, '. Stopping.'
STOP
end if
omega(RxnID)=omega_tmp
end do
close(iomegainp)
else !omega correction off
omega=0.0D0
end if
!======================================================================================================
! Section 3 -- Initialize some miscellaneous stuff
!======================================================================================================
!Set up the time vector
logdt=1/float(ntdec)
logtout=log10(ttout)
ndec=log10(rtime)-logtout !Total number of decades of magnitude spanned
ncalls=ceiling(ndec*ntdec) !Want to round UP to overshoot the stop time
if (.not. allocated(tvec)) allocate(tvec(ncalls+1), stat=ierr)
if (.not. allocated(tvec)) then
write(*,*) 'Time vector not allocated. Stopping.'
STOP
end if
tvec=0.
tvec(ncalls+1)=rtime
do i=2,ncalls
tvec(i)=10**logtout
logtout=logtout+logdt
end do
!Keep track of which reactions have user specied activation energies so
!that if BEPs are used which might overwrite them an error is printed
if (.not. allocated(lEA_user)) allocate(lEA_user(ngrxn+nsrxn))
lEA_user=.false.
!Allocate memory for solver arrays & set problem sizes
liw=40+neqns
select case (iSolver)
case (1) !Krylov solver, DASPK
maxl=max(5,neqns)
lwp=(2*ml+mu+1)*neqns+2*((neqns/(ml+mu+1))+1) !Preconditioner real storage
lrw=50+10*(neqns)+(neqns)**2+maxl+(maxl+3)*neqns+(maxl+3)*maxl+lwp
case (0) !Direct solver, DASPK
lrw=50+9*(neqns)+(neqns)**2
case default !Error condition
write(*,*) 'Invalid solver selection specified. Using direct method, DASPK.'
lrw=50+9*(neqns)+(neqns)**2
end select
if (.not. allocated(atol)) allocate(atol(neqns), stat=ierr)
if (.not. allocated(atol)) mem_stat_total=mem_stat_total+1
if (.not. allocated(rtol)) allocate(rtol(neqns), stat=ierr)
if (.not. allocated(rtol)) mem_stat_total=mem_stat_total+1
if (.not. allocated(iwork)) allocate(iwork(liw), stat=ierr)
if (.not. allocated(iwork)) mem_stat_total=mem_stat_total+1
if (.not. allocated(rwork)) allocate(rwork(lrw), stat=ierr)
if (.not. allocated(rwork)) mem_stat_total=mem_stat_total+1
if (mem_stat_total>0) then
write(*,*) 'Errors initializing storage arrays for DASPK. Stopping.'
STOP
end if
!Initialize solver options
info=0 !Use defaults unless otherwise specified
info(2)=1 !Specifies vector tolerances
if (NonNeg) info(10)=2 !Specifies that Non-negativity constraint is ON during integration
if (iSolver==1) then
info(12)=1 !This turns the Krylov iterative solver /on/
info(15)=1 !This specifies that there is a subroutine to calculate the preconditioner
iwork(27)=lwp
iwork(28)=neqns
end if
!Initialize the BEP routines
if (lBEP) then
call inco_BEP
call zero_BEP_EA
else
if (.not. allocated(BEP_rxn)) allocate(BEP_rxn(nrxn,1))
BEP_rxn=0 !This ensures that we can easily tell which rxns use BEPs and which don't
end if
!Initialize scaling relations
if (iScale>0) then
call inco_Scale
else
nBE_coords=1
end if
!Read the sensitivity analysis file
if (lsenbrute) call inco_tube_SA
!Initialize the coverage parameters
if (lcov) then !Coverage effects (constant values)
call inco_cov
else
cov_matrix=0.0
end if
cov_factor=1.0 !This always starts at 1 and is adjusted as needed
if (iScale<3) cov_factor_adjust=.false. !Only adjust interaction parameters for maps
!Apply the StatpQ correction
if (lStatpQ) then !Approximate corrections for T to constant BE
call inco_StatpQ
else
StatpQ=0.0
end if
!If DOE or global SA perturbation is enabled, allocate memory and save the original correlations
if (lDOE .or. lsenbrute) then
if (perturb_param(globalSA_BEP)) then
if (.not. allocated(BEP_def_save)) allocate(BEP_def_save(nBEP,6), stat=mem_stat)
if (.not. allocated(BEP_def_save)) then
write(*,*) 'Failure to allocate BEP_def_save. Stopping.'
STOP
end if
BEP_def_save=BEP_def
end if
if (perturb_param(globalSA_LSR)) then
if (.not. allocated(scale_slope_save)) allocate(scale_slope_save(nscale,2), stat=mem_stat)
if (.not. allocated(scale_slope_save)) then
write(*,*) 'Failure to allocate scale_slope_save. Stopping.'
STOP
end if
scale_slope_save=scale_slope
end if
end if
!======================================================================================================
! Section 4 -- Write chosen options to screen and sort out possibly conflicting options
!======================================================================================================
!-------------------------------------------------------------------------------------------------------
!Write operating conditions:
select case (irxtr)
case (0)
write(*,*) 'UHV OR MOLECULAR BEAM REACTOR'
write(*,*) ' VELOCITY RESET TO 0 AUTOMATICALLY'
write(*,*) ' TOF CALCULATIONS IN tube_conv.out ARE DUMMY VALUES'
case (1)
write(*,*)'BATCH REACTOR'
write(*,*) ' VELOCITY RESET TO 0 AUTOMATICALLY'
write(*,*) ' TOF CALCULATIONS IN tube_conv.out ARE DUMMY VALUES'
write(*,*) ' TRANSIENT RESULTS AUTOMATICALLY SAVED'
ltra=.true.
case (2)
write(*,*)'CONTINUOUS STIRRED TANK REACTOR'
case (3)
write(*,*)'PLUG FLOW REACTOR'
case default
write(*,*)'INVALID REACTOR SELECTION. STOPPING'
STOP
end select
if (.not. gas_chem) then
write(*,*)'NO GAS CHEMISTRY'
else
write(*,*)'GAS CHEMISTRY PRESENT'
endif
if (.not. surf_chem) then
write(*,*)'NO SURFACE CHEMISTRY'
else
write(*,*)'SURFACE CHEMISTRY PRESENT'
endif
if (itpd>0) then
iiso=1 !TPD specification overrides isothermality flag
liso=.false.
write(*,*) 'TPD SPECIFICATION OVERRIDES ISOTHERMAL SPECIFICATION.'
end if
if (liso) then
write(*,*)'ISOTHERMAL RUN'
else
write(*,*)'NON-ISOTHERMAL RUN'
if (itpd>0) then
write(*,*) 'TPD -- TRANSIENT RESULTS AUTOMATICALLY SAVED'
ltra=.true.
if (itpd==1) then
write(*,*) 'UHV RUN, PRESSURE SET TO 1E-9 TORR. COMPOSITION SET TO 100% INERT.'
write(*,*) 'UHV RUN, UHV REACTOR MODEL USED.'
else
write(*,*) 'HIGH PRESSURE RUN'
end if
else if (htc/=0.0) then
write(*,*)'EXTERNAL HEAT TRANSFER'
else
write(*,*)'ADIABATIC RUN'
endif
endif
if (lomega) write(*,*) 'OMEGA CORRECTION ON'
if (lStatpQ) write(*,*) 'STATPQ T CORRECTION APPLIED TO BINDING ENERGIES'
if (lBEP) write(*,*) 'BEP RELATIONS USED TO ESTIMATE BARRIERS'
if (iScale>0) write(*,*) 'SCALING RELATIONS USED TO ESTIMATE HEATS OF FORMATION'
if (iScale==3) write(*,*) 'SCALING RELATIONS USED TO PERFORM CATALYST SCREENING'
if (lDOE) write(*,*) 'PARAMETER PERTURBATIONS IN EFFECT AFTER INITIAL RUN'
if (lEA) write(*,*) 'EXTERNAL USER-SPECIFIED ACTIVATION ENERGIES'
if (lstp) write(*,*)'INLET VELOCITY CHANGES WITH TEMPERATURE AND PRESSURE'
if (.not. ltra) write(*,*)'!!NOTE!! TRANSIENT RESULTS NOT SAVED'
select case (mrpa)
case (0)
write(*,*)'REACTION PATH ANALYSIS OFF'
verbose_rpa=.false. !No RPA implies no verbose RPA
case (1)
write(*,*)'REACTION PATH ANALYSIS ON'
trpa_called=.true. !No TPD RPA should be performed
case (2)
write(*,*) 'TEMPERATURE SPECIFIC TPD REACTION PATH ANALYSIS REQUESTED'
case default
mrpa=0
verbose_rpa=.false. !No RPA implies no verbose RPA
write(*,*) 'INVALID REACTION PATH ANALYSIS OPTION SELECTED'
write(*,*) 'REACTION PATH ANALYSIS NOT WRITTEN'
end select
if (mrpa>0 .and. verbose_rpa .and. irxtr==3) then
write(*,*) 'REACTION PATH ANALYSIS PERFORMED AT EVERY REACTOR NODE'
else if (mrpa>0) then
write(*,*) 'REACTION PATH ANALYSIS PERFORMED AT REACTOR EXIT'
verbose_rpa=.false. !Batch/CST reactors only have RPA at the exit
end if
if (lDOE .and. lsenbrute) then
write(*,*) 'DESIGN OF EXPERIMENTS AND SENSITIVITY ANALYSIS RUNS INCOMPATIBLE. STOPPING.'
STOP
end if
if (.not. lsenbrute) isenbrute=0
select case (isenbrute)
case (0)
write(*,*)'BRUTE FORCE SENSITIVITY ANALYSIS OFF'
case (1)
write(*,*) 'APPROXIMATE FIM-BASED SENSITIVITY ANALYSIS ON PRE-EXPONENTIALS'
case (2)
write(*,*) 'LOCAL BRUTE FORCE SENSITIVITY ANALYSIS SELECTED'
if (rel_pert) then
write(*,*) 'RELATIVE PERTURBATIONS TO PARAMETER VALUES USED'
else
write(*,*) 'ABSOLUTE PERTURBATIONS TO PARAMETER VALUES USED'
end if
if (count(spec_pert)>0) then
if (fix_EA) then
write(*,*) 'ACTIVATION ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
else
write(*,*) 'TRANSITION STATE ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
end if
end if
case (3)
write(*,*) 'VARIANCE-BASED GLOBAL SENSITIVITY ANALYSIS SELECTED'
if (count(spec_pert)>0) then
if (fix_EA) then
write(*,*) 'ACTIVATION ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
else
write(*,*) 'TRANSITION STATE ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
end if
end if
case (4)
write(*,*) 'DERIVATIVE-BASED GLOBAL SENSITIVITY ANALYSIS SELECTED'
if (count(spec_pert)>0) then
if (fix_EA) then
write(*,*) 'ACTIVATION ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
else
write(*,*) 'TRANSITION STATE ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
end if
end if
case (5)
write(*,*) 'UNCERTAINTY QUANTIFICATION ONLY'
if (count(spec_pert)>0) then
if (fix_EA) then
write(*,*) 'ACTIVATION ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
else
write(*,*) 'TRANSITION STATE ENERGIES FIXED WHEN SPECIES PROPERTIES PERTURBED'
end if
end if
end select
if (lsenbrute .and. isenbrute>2) then
if (.not. model_solve) write(*,*) 'NO MODEL SOLUTION -- RANDOM SAMPLING OF DESIGN POINTS ONLY'
end if
if (lcov) then
write(*,*) 'COVERAGE EFFECTS ON'
if (cov_factor_adjust) then
write(*,*) ' COVERAGE EFFECT PARAMETERS ADJUSTED FOR ATOMIC BINDING ENERGIES'
else
write(*,*) ' COVERAGE EFFECT PARAMETERS NOT ADJUSTED FOR ATOMIC BINDING ENERGIES'
end if
end if
if (nruns > 1) then
write(*,*)'MULTIPLE SETS OF REACTOR INPUT CONDITIONS'
if (.not. MultiInput) then
if (lDOE) write(*,*) 'DOE PERTURBATIONS BASED ON A SINGLE PARAMETER SET SPECIFIED.'
if (trise/=0) write(*,*) 'TEMPERATURE RAMP SPECIFIED.'
if (.not. lDOE .and. trise==0) then
write(*,*) 'TEMPERATURE RISE OF 0 SPECIFIED FOR NON-DOE RUN. STOPPING.'
STOP
end if
end if
else
write(*,*)'SINGLE SET OF REACTOR INPUT CONDITIONS'
endif
if (Tref_beta == 0) then
write(*,*)'T_ref FOR BETA TERM OF KINETIC RATE CONSTANTS = 300 K'
else
write(*,*)'T_ref FOR BETA TERM OF KINETIC RATE CONSTANTS = 1 K'
endif
write(*,*)'************ MultiSite Variable Check ************'
write(*,*)'number of site types:', nsphases
write(*,*)'site type IDs: ', site_type_lastID(2:nsphases+1)
!-------------------------------------------------------------------------------------------------------
if (itpd==1) press=1.32E-12 !1E-9 torr (typical UHV pressure)
p=press*patm ! atm --> dyne/cm2
!-------------------------------------------------------------------------------------------------------
!*********************************************************************************************************
end subroutine inco_tube
!*********************************************************************************************************
!*********************************************************************************************************
subroutine get_tube_conditions
!This subroutine is a wrapper which calls all the other subroutines to get
!run-specific parameters and store them in global arrays.
!*********************************************************************************************************
implicit none
call inco_Tflow !Operating conditions (T, P, etc.)
call inco_tube_mole !Feed composition
call inco_EA !User-specified activation energies
call inco_ddaspk !Solver parameters, tolerances, etc.
!*********************************************************************************************************
end subroutine get_tube_conditions
!*********************************************************************************************************
!*********************************************************************************************************
subroutine update_tube_conditions
!This subroutine extracts the run-specific parameters from the global arrays.
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles
implicit none
integer :: i, j, nrows
double precision :: addo, t
double precision :: EA_tmp
!Feed condition
tin=TPFSV(i_tube,1)
p=TPFSV(i_tube,2)
velo=TPFSV(i_tube,3)
abyv=TPFSV(i_tube,4)
t=tin
!Check to see if we need to adjust for STP
if (lstp) velo=velo*(t/273.15)*(patm/p) ! Account for T,P effect on velocity
!Calculate the density and the mass flux
call ckrhoy(p,t,yin,ickwrk,rckwrk,rhoin) ! Inlet density
fluxmass=rhoin*velo ! Mass flux
!Feed composition
!The rest of this is useful for converting the act_all values into the actual compositions...
actin(1:kgmax+ksmax)=0.0D0 !Initialize to zero and set only non-zero entries
!Extract the input compositiion
nrows=size(act_all,1)
if (itube_restart/=0) nrows=nrows-1
do i=1,nrows
actin(iMole_spec(i))=act_all(i,i_tube)
end do
!Correct for UHV TPD
if (itpd==1) then
forall (i=1:kgmax-1) actin(i)=0.0D0
actin(kgmax)=1.0D0
end if
!Normalize compositions for all phases to have activities summing to 1
!Activities are mole fractions for gas & coverages for each type of surface site
do_999: do i=1,nsphases+1 ! Loop over all surface site types
addo=sum(actin(site_type_firstID(i):site_type_lastID(i))) ! sums the initial values of the (i-1)th site
if (abs(1-addo)>1E-4) then
write(*,*) 'Warning: possible error in composition in phase ', &
PhaseNames(i),' due to renormalization. Continuing.'
end if
actin(site_type_firstID(i):site_type_lastID(i))=actin(site_type_firstID(i):site_type_lastID(i))/addo ! normalize coverages for (i+1)th site
end do do_999
actin(kmax-kbmax+1:kmax)=1.0D0 !Bulk species have activity of 1
act=actin
x(1:kgmax)=act(1:kgmax)
xin=x
call ckxty(x,ickwrk,rckwrk,y)
call ckxty(xin,ickwrk,rckwrk,yin)
call ckrhoy(p,t,yin,ickwrk,rckwrk,rhoin) ! Inlet density
win(1:kgmax)=yin
forall (i=kgmax+1:kmax-kbmax) win(i)=actin(i)
win(indext)=t !Temperature
win(indexrho)=rhoin !Density
fluxmass=rhoin*velo ! Mass flux
if (itube_restart/=0) fluxmass_inlet=act_all(nrows,i_tube)
!Activation energies
if (lEA) then
if (gas_chem) then
nrows=size(iEAg_rxn,1)
do i=1,nrows
j=iEAg_rxn(i)
EA_tmp=EAg_all(i,i_tube)*tin !Use Tin, not T since the EA values are divided by RTin
call ckreex(-j,rckwrk,EA_tmp)
end do
end if
if (surf_chem) then
nrows=size(iEAs_rxn,1)
do i=1,nrows
j=iEAs_rxn(i)
EA_tmp=EAs_all(i,i_tube)*tin
call skreex(-j,iskwrk,rskwrk,EA_tmp)
end do
end if
end if
!Solver tolerances
if (ltol) then
nrows=size(iTol_spec,1)
do i=1,nrows
atol(iTol_spec(i))=atol_all(i,i_tube)
rtol(iTol_spec(i))=rtol_all(i,i_tube)
end do
end if
!*********************************************************************************************************
end subroutine update_tube_conditions
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_Tflow
!This subroutine is responsible for reading conditions from T_flow.
!*********************************************************************************************************
use global_all; use file_handles; use global_tube1; use update_tube2
implicit none
integer :: i, ierr
double precision :: press
character(len=80) :: crap
if (.not. allocated(TPFSV)) allocate(TPFSV(nruns,4), stat=ierr)
if (.not. allocated(TPFSV)) then
write(*,*) 'Feed condition array not allocated. Stopping.'
STOP
end if
if (.not. MultiInput) then !Assign the tube.inp values to the TPFSV array
do i=1,nruns
TPFSV(i,1)=tin+(i-1)*trise
end do
TPFSV(:,2)=p !Already in dyne/cm^2 (converted in inco_tube)
TPFSV(:,3)=velo
TPFSV(:,4)=abyv
else !Read from T_flow.inp
open(unit=iTflowinp, file='INP.d/T_flow.inp', status='unknown', action='read')
call skip_comment(iTflowinp)
do i=1,nruns
read(iTflowinp,*) TPFSV(i,:) !tin,press,velo,abyv -- Current operating conditions
end do
TPFSV(:,2)=TPFSV(:,2)*patm ! atm --> dyne/cm2
close(iTflowinp)
end if
!*********************************************************************************************************
end subroutine inco_Tflow
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_tube_mole
!This subroutine is responsible for reading conditions from T_flow.
!*********************************************************************************************************
use global_all; use file_handles; use global_tube1; use update_tube2
implicit none
integer :: i, ii, nlines, ierr
character(len=80) :: crap
double precision :: act_local
!Read inlet mole fractions and coverages (ignore bulk species)
open(unit=itubemoleinp, file='INP.d/tube_mole.inp', status='unknown', action='read')
call skip_comment(itubemoleinp)
read(itubemoleinp,*)itube_restart
call skip_comment(itubemoleinp)
read(itubemoleinp,*)nlines
call skip_comment(itubemoleinp)
!Check for the restart and adjust the number of lines read
if (itube_restart>0) nlines=nlines-1
!Initialize the vector mapping species ID numbers to input file lines
if (.not. allocated(iMole_spec)) allocate(iMole_spec(nlines), stat=ierr)
if (.not. allocated(iMole_spec)) then
write(*,*) 'Species ID vector not allocated. Stopping.'
STOP
end if
call spec_vec_init(nlines,itubemoleinp,iMole_spec)
!Now allocate the memory for the stored feed composition array
if (itube_restart>0) nlines=nlines+1 !Account for mass flux for restart
if (.not. allocated(act_all)) allocate(act_all(nlines,nruns), stat=ierr)
if (.not. allocated(act_all)) then
write(*,*) 'Feed composition array not allocated. Stopping.'
STOP
end if
!Note: if itube_restart>0, then the last row in act_all is actually the mass flux!
do ii=1,nlines
if (MultiInput) then
read(itubemoleinp,*)crap, act_all(ii,:)
else
read(itubemoleinp,*)crap, act_local
act_all(ii,:)=act_local
end if
end do
close(itubemoleinp)
!*********************************************************************************************************
end subroutine inco_tube_mole
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_EA
!This subroutine gets user-specified activation energies from EA?.inp
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles
implicit none
integer :: irxn, jrxn, rxn_count, i
character(len=80) :: crap
if (.not. lEA) return
!Read in user specified activation energies for gas reactions and reset
!the corresponding value in the work arrays.
if (gas_chem) then
open(unit=iEAginp, file='INP.d/EAg.inp', status='unknown', action='read')
call skip_comment(iEAginp)
read(iEAginp,*) rxn_count
call skip_comment(iEAginp)
!Initialize the vector mapping reaction ID numbers to input file lines
if (.not. allocated(iEAg_rxn)) allocate(iEAg_rxn(rxn_count))
call rxn_vec_init(rxn_count,iEAginp,iEAg_rxn)
!Initialize the storage arrays for the EAg values
if (.not. allocated(EAg_all)) allocate(EAg_all(rxn_count,nruns))
do irxn=1,rxn_count
if (MultiInput) then
read(iEAginp,*) crap, EAg_all(irxn,:)
else
read(iEAginp,*) crap, EAg_all(irxn,1)
EAg_all(irxn,2:nruns)=EAg_all(irxn,1)
end if
jrxn=iEAg_rxn(irxn)
lEA_user(jrxn)=.true.
end do
close(iEAginp)
end if
!Read in user specified activation energies for surface reactions and reset
!the corresponding value in the work arrays.
if (surf_chem) then
open(unit=iEAsinp, file='INP.d/EAs.inp', status='unknown', action='read')
call skip_comment(iEAsinp)
read(iEAsinp,*) rxn_count
call skip_comment(iEAsinp)
!Initialize the vector mapping reaction ID numbers to input file lines
if (.not. allocated(iEAs_rxn)) allocate(iEAs_rxn(rxn_count))
call rxn_vec_init(rxn_count,iEAsinp,iEAs_rxn)
!Initialize the storage arrays for the EAs values
if (.not. allocated(EAs_all)) allocate(EAs_all(rxn_count,nruns))
!Accounts for the fact that the rxn_str vector has both gas and surface
!reaction ID numbers but the SK work arrays have only surface reactions.
iEAs_rxn=iEAs_rxn-ngrxn
do irxn=1,rxn_count
if (MultiInput) then
read(iEAsinp,*) crap, EAs_all(irxn,:)
else
read(iEAsinp,*) crap, EAs_all(irxn,1)
EAs_all(irxn,2:nruns)=EAs_all(irxn,1)
end if
jrxn=iEAs_rxn(irxn)
lEA_user(ngrxn+jrxn)=.true.
end do
close(iEAsinp)
end if
!*********************************************************************************************************
end subroutine inco_EA
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_ddaspk
!This subroutine gets initial parameters for ddaspk as well as allocates memory for arrays necessary
!to solving the problem.
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles
implicit none
double precision :: tol_tmp(2*nruns)
integer :: i, j, ierr, nlines, mem_stat_total=0
character(len=80) :: crap
logical :: lenergy
!-------------------------------------------------------------------------------------------------------
!Set solver tolerances
atol=abstol; rtol=reltol
if (.not. ltol) return !Only read the rest of the input for granular tolerances
open(unit=itolinp, file='INP.d/tol.inp', status='old', action='read')
call skip_comment(itolinp)
read(itolinp,*) nlines
if (nlines<=0) then !There's nothing else to do here
close(itolinp)
return
end if
call skip_comment(itolinp)
read(itolinp,*) lenergy
call skip_comment(itolinp)
!Initialize storage and species ID arrays if needed
if (.not. allocated(iTol_spec)) allocate(iTol_spec(nlines), stat=ierr)
if (.not. allocated(iTol_spec)) then
write(*,*) 'Species ID vector not allocated. Stopping.'
STOP
end if
if (.not. allocated(atol_all)) allocate(atol_all(nlines,nruns), stat=ierr)
if (.not. allocated(atol_all)) then
write(*,*) 'Absolute tolerance storage vector not allocated. Stopping.'
STOP
end if
if (.not. allocated(rtol_all)) allocate(rtol_all(nlines,nruns), stat=ierr)
if (.not. allocated(rtol_all)) then
write(*,*) 'Absolute tolerance storage vector not allocated. Stopping.'
STOP
end if
if (lenergy) then !Last line is energy balance
call spec_vec_init(nlines-1,itolinp,iTol_spec)
iTol_spec(nlines)=neqns
else
call spec_vec_init(nlines,itolinp,iTol_spec)
end if
!Read in tolerances
do j=1,nlines !Assume that the /last/ line is the energy balance
if (MultiInput) then
read(itolinp,*)crap,tol_tmp(:)
atol_all(j,:)=tol_tmp(1:2*nruns-1:2) !Odd elements
rtol_all(j,:)=tol_tmp(2:2*nruns:2) !Even elements
else
read(itolinp,*)crap,tol_tmp(1:2)
atol_all(j,:)=tol_tmp(1)
rtol_all(j,:)=tol_tmp(2)
end if
end do
close(itolinp)
!*********************************************************************************************************
end subroutine inco_ddaspk
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_BEP
!This subroutine reads BEP.inp and stores the necessary variables to global before being used in skBEP
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles; use globalSA
implicit none
integer :: del, i, mem_stat, mem_stat_total=0, ierr, nlines, j
integer :: RxnID, NT, BEP_rxn_tmp(nrxn,3)
character(len=80) :: crap, rxn_str_tmp
!-------------------------------------------------------------------------------------------------------
write(*,*)'************************************************'
!.........Open input files:
open(unit=iBEPinp, file='INP.d/BEP.inp', status='old', action='read')
!-------------------------------------------------------------------------------------------------------
call skip_comment(iBEPinp)
read(iBEPinp,*) nBEP !read number of correlations
!Allocate the necessary storage space
if (.not. allocated(BEP_def)) allocate(BEP_def(nBEP,6), stat=mem_stat)
if (.not. allocated(BEP_def)) mem_stat_total=mem_stat_total+1
if (.not. allocated(BEP_rxn)) allocate(BEP_rxn(nrxn,2), stat=mem_stat)
if (.not. allocated(BEP_rxn)) mem_stat_total=mem_stat_total+1
if (mem_stat_total>0) then
write(*,*) 'Memory allocation failure in inco_BEP, stopping'
STOP
end if
BEP_def=0.0D0
BEP_rxn=0
call skip_comment(iBEPinp)
do i=1,nBEP
!read in correlation info (type, m, b, direction)
read(iBEPinp,*) del, (BEP_def(i,j), j=1,6)
if ((BEP_def(i,1)<-1) .OR. BEP_def(i,1)>1) then
write(*,*) 'Invalid correlation type in BEP ',i,'. Stopping'
STOP
end if
if (int(abs(BEP_def(i,4)))/=1) then
write(*,*) 'Invalid correlation direction in BEP',i,'. Stopping'
STOP
end if
end do
!Normalize by R so that energies are stored internally in K
BEP_def(:,3)=BEP_def(:,3)/Rgas_kcal
BEP_def(:,5:6)=BEP_def(:,5:6)/Rgas_kcal
call skip_comment(iBEPinp)
read(iBEPinp,*) nlines
call skip_comment(iBEPinp)
call BEP_rxn_read(nlines,nBEP,iBEPinp,BEP_rxn_tmp(1:nlines,:))
!Assign the data
do i=1,nlines
BEP_rxn(BEP_rxn_tmp(i,3),:)=BEP_rxn_tmp(i,1:2)
end do
close(iBEPinp)
!*********************************************************************************************************
end subroutine inco_BEP
!*********************************************************************************************************
!*********************************************************************************************************
subroutine zero_BEP_EA
!This subroutine reinitializes the activation energies for reactions
!using BEP correlations to 0.
!*********************************************************************************************************
use global_all;
implicit none
integer :: i
double precision :: zero
zero=0.0
if (nBEP==0) return
!Go through the reactions resetting original activation energies of reactions
!using BEP correlations to zero
do i=1,ngrxn
if (BEP_rxn(i,1)/=0) then !This reaction does not use a BEP
call ckreex(-i,rckwrk,zero)
reg(i)=0.0
end if
end do
do i=1,nsrxn
if (BEP_rxn(ngrxn+i,1)/=0) then !This reaction does not use a BEP
call skreex(-i,iskwrk,rskwrk,zero)
re(i)=0.0
end if
end do
!*********************************************************************************************************
end subroutine zero_BEP_EA
!*********************************************************************************************************
!*********************************************************************************************************
subroutine inco_Scale
!This subroutine reads Scale.inp and stores the necessary variables to global before being used in skScale
!*********************************************************************************************************
use global_all; use global_tube1; use update_tube2; use file_handles; use globalSA
implicit none
double precision, allocatable :: BE_range(:,:)
logical, allocatable :: single_point(:)
integer, allocatable :: spec(:)
integer :: i, grid_type, nspec, skip
character(len=80) :: spec_name
double precision :: BE_delta
!This maps the species IDs to entries in the scale_ref array
integer :: specID(kmax)
!-------------------------------------------------------------------------------------------------------
open(unit=iScaleinp, file='INP.d/Scale.inp',status='old', action='read')
!Read the number of BE descriptors, their ranges, and references
call skip_comment(iScaleinp)
read(iScaleinp,*) natoms_scale
!Allocate memory
if (.not. allocated(scale_ref)) allocate(scale_ref(ksmax))
if (.not. allocated(scale_targ)) allocate(scale_targ(natoms_scale,2))
if (.not. allocated(BE_range)) allocate(BE_range(2,natoms_scale))
if (.not. allocated(single_point)) allocate(single_point(natoms_scale))
if (.not. allocated(spec)) allocate(spec(natoms_scale))
call spec_vec_init(natoms_scale,iScaleinp,spec) !This gets the atom names
scale_targ(:,1)=dble(spec) !Save the species ID numbers
scale_ref=1.0
!This gets the data
specID=0
do i=1,natoms_scale
specID(spec(i))=i
read(iScaleinp,*) spec_name, BE_range(:,i), scale_ref(spec(i)-kgmax), &
single_point(i)
!The following ensures that a single point calculation for an atomic
!binding energy will result in a range restricted to the lower value.
if (single_point(i) .or. iScale<3) BE_range(2,i)=BE_range(1,i)
end do
deallocate(spec) !Deallocate so we can reallocate later
!Read the grid options and setup the grid
call skip_comment(iScaleinp)
if (all(single_point) .or. iScale<3) then !Single point only
nBE_coords=1
if (.not. allocated(BE_coords)) allocate(BE_coords(natoms_scale,nBE_coords))
BE_coords(:,1)=BE_range(1,:)
read(iScaleinp,*) grid_type !Skip grid type
if (grid_type/=4) then
read(iScaleinp,*) !Skip grid parameters
else
read(iScaleinp,*) nBE_coords
do i=1,nBE_coords
read(iScaleinp,*) !Skip specified coordinates
end do
nBE_coords=1
end if
else !Binding energy map
read(iScaleinp,*) grid_type
if (grid_type==0) then !Default
if (natoms_scale>3) then !More than 3 dimensions -- use Sobol'
grid_type=3
else if (natoms_scale==1) then
grid_type=1 !Use rectangular grid for 1-D case
else !Use hex-grid for 2-D/3-D cases
grid_type=2
end if
else if (grid_type==2) then
if (natoms_scale>3) then
write(*,*) 'Fatal error: Hex grid not implemented for more than three dimensions. Stopping.'
STOP
end if
if (natoms_scale==1) grid_type=1 !Switch to rectangular grid
end if
select case (grid_type)
case (1, 2) !Rectangular, Hex
read(iScaleinp,*) BE_delta !BE increment
if (grid_type==1) then
call make_rect_grid(BE_delta,BE_range)
else
call make_hex_grid(BE_delta,BE_range)
end if
case (3) !Sobol'
read(iScaleinp,*) nBE_coords, skip !Number of binding energy coordinates