-
Notifications
You must be signed in to change notification settings - Fork 0
/
mode_soatinit.f
155 lines (113 loc) · 4.63 KB
/
mode_soatinit.f
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
module mode_soatinit
use modd_glo
use modd_glodef
implicit none
!Purpose: Contains subroutines to initialize all parameters which are
!only Temperature-dependent
contains
!**********************************
subroutine VP_GET(
& TEMPK !I [K] Temperature
& ,VP !I [torr] saturation vapor pressures
& )
!Purpose: Find the saturation vapor pressures corrected for temperature
implicit none
!Input
REAL, DIMENSION(:), INTENT(IN) :: TEMPK !I [K] Temperature
!Output
REAL, DIMENSION(:,:), INTENT(OUT) :: VP !O [torr] saturation vapor pressures
!Local
REAL, DIMENSION(SIZE(VP,2)) :: DELCP
REAL, DIMENSION(SIZE(VP,2)) :: DELSB
REAL, DIMENSION(SIZE(VP,1),SIZE(VP,2)) :: TERM1
REAL, DIMENSION(SIZE(VP,1),SIZE(VP,2)) :: TERM2
INTEGER :: I ![idx] counter for component
!The following is cut and paste from Griffin's code
!I have to check the units here!
DO I = 1,NBSP
!Calculate DELSB
DELSB(I)= 86.0 + 0.4*TAUVP(I)+1421.*HBN(I)
!Calculate DELCP
DELCP(I) = -90.-2.1*TAUVP(I)
!Calculate Intermediate terms 1
TERM1(:,I) = DELSB(I)*(TEMPK(:)-TBOIL(I))
& /19.1/TEMPK(:)
!Calcualte intermediate term 2
TERM2(:,I) = ((TBOIL(I)-TEMPK(:))/TEMPK(:)
& -LOG(TBOIL(I)/TEMPK(:)))
& *DELCP(I)/19.1
!Get the temperature-dependent
VP(:,I) = 760.*(10.**(TERM1(:,I)+TERM2(:,I)))
ENDDO
!Tuning: Correct vp so that they match chamber data
VP(:,1) = VP(:,1)/1.5
VP(:,2) = VP(:,2)/1.4
VP(:,3) = VP(:,3)/2.4
VP(:,5) = VP(:,5)/66./3.0 !by james
VP(:,8) = VP(:,8)/33.
VP(:,9) = VP(:,9)/3.3
VP(:,10) = VP(:,10)/15.
VP(:,11) = VP(:,11)/5000. !by james
end subroutine VP_GET
!*********************************
subroutine AQCONST_GET(
& TEMPK !I [K] temperature
& ,K !O [units] Henry's law and dissociation constants
& )
!Purpose: Calculate
implicit none
!Input
REAL, INTENT(IN), DIMENSION(:) :: TEMPK !I [K] temperature
!Output
REAL, INTENT(OUT), DIMENSION(:,:) :: K !O [units] temperature dependent value of K_298
!Local
INTEGER :: I ![idx] counter for components
INTEGER :: COMP_IDX ![idx] pointer to right component
INTEGER :: J ![idx] counter for sub components (ions)
COMP_IDX=1
DO I=1,NBSP
!Adjust Henry's law constant
K(:,COMP_IDX) = K_298(COMP_IDX)
& *EXP(15.*((1./TEMPK(:))-(1./298.15))/0.001987)
!Prepare for next component
COMP_IDX=COMP_IDX+1
!Calculate dissociation constants of sub-components (ions)
DO J=2,NK(I)
K(:,COMP_IDX) = K_298(COMP_IDX)
& *EXP(0.5*((1./TEMPK(:))-(1./298.15))/0.001987)
!Prepare for next component
COMP_IDX=COMP_IDX + 1
ENDDO
ENDDO !loop on I
end subroutine AQCONST_GET
!*************************************
subroutine SI_GET(
& TEMPK !I [K] temperature
& ,A !I [units?] term in UNIFAC parameterization
& ,SI !O [units?] term in UNIFAC parameterization
& ,NFUNC !I [nbr] number of functional group in mixture
& )
!Purpose: Calculate the term in the unifac method which
!is dependent only on temperature. This term includes
!exponential of a large matrix, so it presumably takes
!up a lot of computer time.
implicit none
!Input
REAL, DIMENSION(:), INTENT(IN) :: TEMPK !I [K] temperature
REAL, DIMENSION(:,:), INTENT(IN) :: A !I [?] Term in UNIFAC parameterization
INTEGER, INTENT(IN) :: NFUNC !I [nbr] number of functional groups in mixture
!Output
REAL, DIMENSION(:,:,:), INTENT(OUT) :: SI !O [?] term in UNIFAC parameterization
!Local
REAL, DIMENSION(SIZE(TEMPK)) :: TEMPKINV ![1/K] inverse of temperature
INTEGER :: J1, J2 ![idx] counter for functional groups
!Get values of SI (temperature dependent), don't need to be in the iteration!
!eqn 9 in MaP05
TEMPKINV(:) = 1.d0/TEMPK(:)
DO J1=1,NFUNC
DO J2=1,NFUNC
SI(:,J1,J2) = EXP(-A(J1,J2)*TEMPKINV(:))
ENDDO
ENDDO
end subroutine SI_GET
end module mode_soatinit