← Back to branch summary

~fluidity-core/fluidity/exorcised

« back to all changes in this revision

Viewing changes to legacy_reservoir_prototype/src/multi_eos.F90

  • Committer: saunde01
  • Date: 2011-03-23 13:15:31 UTC
  • Revision ID: svn-v4:5bf5533e-7014-46e3-b1bb-cce4b9d03719:trunk:3310
As discussed in the Dev meeting and in emails, this commit renames the multiphase directory to legacy_reservoir_prototype, and adds an obvious warning to the top of the README file to warn unsuspecting users of the experimental nature of the code. Iterations on the exact wording very welcome

Show diffs side-by-side

added added

removed removed

Lines of Context:
legacy_reservoir_prototype/src/multi_eos.F90
 
1
!    Copyright (C) 2006 Imperial College London and others.
 
2
!    
 
3
!    Please see the AUTHORS file in the main source directory for a full list
 
4
!    of copyright holders.
 
5
!
 
6
!    Prof. C Pain
 
7
!    Applied Modelling and Computation Group
 
8
!    Department of Earth Science and Engineering
 
9
!    Imperial College London
 
10
!
 
11
!    amcgsoftware@imperial.ac.uk
 
12
!    
 
13
!    This library is free software; you can redistribute it and/or
 
14
!    modify it under the terms of the GNU Lesser General Public
 
15
!    License as published by the Free Software Foundation,
 
16
!    version 2.1 of the License.
 
17
!
 
18
!    This library is distributed in the hope that it will be useful,
 
19
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
!    Lesser General Public License for more details.
 
22
!
 
23
!    You should have received a copy of the GNU Lesser General Public
 
24
!    License along with this library; if not, write to the Free Software
 
25
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
!    USA
 
27
 
 
28
module multiphase_EOS
 
29
 
 
30
contains
 
31
 
 
32
  SUBROUTINE CAL_BULK_DENSITY( NPHASE, NCOMP, CV_NONODS, CV_PHA_NONODS, NCOEF, DEN, DERIV, &
 
33
       T, P, MASSFRAC, EOS_COEFS, EOS_OPTION )
 
34
    IMPLICIT NONE
 
35
    INTEGER, intent( in ) :: NPHASE, NCOMP, CV_NONODS, CV_PHA_NONODS, NCOEF
 
36
    REAL, DIMENSION( CV_PHA_NONODS ), intent( inout ) :: DEN, DERIV
 
37
    REAL, DIMENSION( CV_PHA_NONODS ), intent( in ) :: T
 
38
    REAL, DIMENSION( CV_NONODS ), intent( in ) :: P
 
39
    REAL, DIMENSION( CV_PHA_NONODS * NCOMP ), intent( in ) :: MASSFRAC
 
40
    REAL, DIMENSION( NPHASE, NCOEF ), intent( in ) :: EOS_COEFS
 
41
    INTEGER, DIMENSION( NPHASE ), intent( in ) ::  EOS_OPTION
 
42
    ! Local
 
43
    INTEGER :: ICOMP, CV_NOD, IPHASE, NOD, NCOMP2
 
44
 
 
45
    write(357,*) 'In CAL_BULK_DENSITY'
 
46
 
 
47
    IF( NCOMP == 0 ) THEN
 
48
       NCOMP2 = 1
 
49
    ELSE
 
50
       NCOMP2 = NCOMP
 
51
    END IF
 
52
 
 
53
    Loop_Comp: DO ICOMP = 1, NCOMP2
 
54
 
 
55
       CALL CAL_DENSITY( NPHASE, ICOMP, CV_NONODS, CV_PHA_NONODS, NCOEF, DEN, DERIV, &
 
56
            T, P, EOS_COEFS, EOS_OPTION )
 
57
 
 
58
    END DO Loop_Comp
 
59
 
 
60
    write(357,*) 'Leaving CAL_BULK_DENSITY'
 
61
 
 
62
    RETURN
 
63
  END SUBROUTINE CAL_BULK_DENSITY
 
64
 
 
65
  SUBROUTINE CAL_DENSITY( NPHASE, ICOMP, CV_NONODS, CV_PHA_NONODS, NCOEF, DEN, DERIV, &
 
66
       T, P, EOS_COEFS, EOS_OPTION ) 
 
67
    ! This sub calculates the DEN ,DERIV for multiphase flows i.e. the EoS 
 
68
    ! and the derivative. 
 
69
    IMPLICIT NONE
 
70
    INTEGER, intent( in ) :: NPHASE, ICOMP, CV_NONODS, CV_PHA_NONODS, NCOEF
 
71
    REAL, DIMENSION( CV_PHA_NONODS ), intent( inout ) :: DEN, DERIV
 
72
    REAL, DIMENSION( CV_PHA_NONODS ), intent( in ) :: T
 
73
    REAL, DIMENSION( CV_NONODS ), intent( in ) :: P
 
74
    REAL, DIMENSION( NPHASE, NCOEF ), intent( in ) :: EOS_COEFS
 
75
    INTEGER, DIMENSION( NPHASE ), intent( in ) ::  EOS_OPTION
 
76
    ! Local
 
77
    REAL, PARAMETER :: TOLER = 1.0E-10
 
78
    REAL :: DEN_P1, DEN_M1, PERT_P
 
79
    INTEGER :: IPHASE, CV_NOD, NODE
 
80
 
 
81
    write(357,*) 'In CAL_DENSITY'
 
82
 
 
83
    Loop_Phase: DO IPHASE = 1, NPHASE
 
84
 
 
85
       Loop_CV: DO CV_NOD = 1, CV_NONODS
 
86
 
 
87
          NODE = CV_NOD + ( IPHASE - 1 ) * CV_NONODS
 
88
 
 
89
          CALL CAL_DENSITY_EOS( NCOEF, DEN( NODE ), T( NODE ), P( CV_NOD ), &
 
90
               EOS_COEFS( IPHASE, : ), EOS_OPTION( IPHASE ) )  
 
91
 
 
92
          IF( EOS_OPTION( IPHASE ) == 1) THEN
 
93
             PERT_P = 0.001 * ( ABS(P( CV_NOD )) + ABS(EOS_COEFS( IPHASE, 1 )) )
 
94
          ELSE
 
95
             PERT_P = 0.001 * ABS(P( CV_NOD ))
 
96
          ENDIF
 
97
          PERT_P =MAX(TOLER, PERT_P)
 
98
 
 
99
          CALL CAL_DENSITY_EOS( NCOEF, DEN_P1, T( NODE ), P( CV_NOD ) + PERT_P, &
 
100
               EOS_COEFS( IPHASE, : ), EOS_OPTION( IPHASE ) ) 
 
101
          CALL CAL_DENSITY_EOS( NCOEF, DEN_M1, T( NODE ), P( CV_NOD ) - PERT_P, &
 
102
               EOS_COEFS( IPHASE, : ), EOS_OPTION( IPHASE ) ) 
 
103
 
 
104
          ! Calculating d(DEN) / dP
 
105
          DERIV( NODE ) = ( DEN_P1 - DEN_M1 ) / ( 2. * PERT_P )
 
106
 
 
107
       END DO Loop_CV
 
108
 
 
109
    END DO Loop_Phase
 
110
 
 
111
    write(357,*) 'Leaving CAL_DENSITY'
 
112
 
 
113
    RETURN
 
114
 
 
115
  END SUBROUTINE CAL_DENSITY
 
116
 
 
117
 
 
118
 
 
119
  SUBROUTINE CAL_DENSITY_EOS( NCOEF, DEN, T, P , &
 
120
       EOS_COEFS, EOS_OPTION )
 
121
    IMPLICIT NONE
 
122
    INTEGER, intent( in ) :: NCOEF
 
123
    REAL, intent( inout ) :: DEN
 
124
    REAL, intent( in ) :: T, P
 
125
    REAL, DIMENSION( NCOEF ), intent( in ) :: EOS_COEFS
 
126
    INTEGER, intent( in ) :: EOS_OPTION
 
127
 
 
128
    IF( EOS_OPTION == 1 ) THEN
 
129
       ! stiffened gas EoS...
 
130
       DEN = ( P + EOS_COEFS( 1 )) * EOS_COEFS( 2 ) / T
 
131
    ELSEIF( EOS_OPTION == 2 ) THEN
 
132
       ! Linear EoS in pressure and temperature...
 
133
 
 
134
       DEN =  EOS_COEFS( 1 )                 + EOS_COEFS( 2 ) * P             + EOS_COEFS( 3 ) * T           &
 
135
            + EOS_COEFS( 4 ) * P * T         + EOS_COEFS( 5 ) * ( P **2 )     + EOS_COEFS( 6 ) * ( T **2 )   &
 
136
            + EOS_COEFS( 7 ) * ( P **2 ) * T + EOS_COEFS( 8 ) * P * ( T **2 ) + EOS_COEFS( 9 ) * (( P * T ) **2 )
 
137
 
 
138
    ENDIF
 
139
 
 
140
 
 
141
    RETURN
 
142
  END SUBROUTINE CAL_DENSITY_EOS
 
143
 
 
144
 
 
145
 
 
146
 
 
147
  SUBROUTINE CAL_CPDEN( NPHASE, CV_NONODS, CV_PHA_NONODS, CPDEN, DEN, NCP_COEFS, CP_COEFS, CP_OPTION, STOTEL, CV_SNLOC, SUF_CPD_BCU, SUF_D_BCU ) 
 
148
 
 
149
    ! This sub calculates the CPDEN ie. CP*DEN
 
150
    IMPLICIT NONE
 
151
    INTEGER, intent( in ) :: NPHASE, CV_NONODS, CV_PHA_NONODS, NCP_COEFS, STOTEL, CV_SNLOC
 
152
    REAL, DIMENSION( CV_PHA_NONODS ), intent( in ) :: DEN
 
153
    REAL, DIMENSION( CV_PHA_NONODS ), intent( inout ) :: CPDEN
 
154
    REAL, DIMENSION( STOTEL * CV_SNLOC * NPHASE ), intent( in ) :: SUF_D_BCU
 
155
    REAL, DIMENSION( STOTEL * CV_SNLOC * NPHASE ), intent( inout ) :: SUF_CPD_BCU
 
156
    REAL, DIMENSION( NPHASE, NCP_COEFS ), intent( in ) :: CP_COEFS
 
157
    INTEGER, DIMENSION( NPHASE ), intent( in ) ::  CP_OPTION
 
158
    ! Local
 
159
    INTEGER :: IPHASE, CV_NOD, CV_NOD_IPHA, IS_IPHA, II
 
160
 
 
161
    Loop_Phase1: DO IPHASE = 1, NPHASE
 
162
 
 
163
       Loop_CV: DO CV_NOD = 1, CV_NONODS
 
164
 
 
165
          CV_NOD_IPHA = CV_NOD + ( IPHASE - 1 ) * CV_NONODS
 
166
 
 
167
          IF( CP_OPTION( IPHASE ) == 1 ) THEN
 
168
             CPDEN( CV_NOD_IPHA ) = CP_COEFS( IPHASE, 1 ) * DEN( CV_NOD_IPHA )
 
169
          ELSE
 
170
             STOP 3931
 
171
          ENDIF
 
172
 
 
173
       END DO Loop_CV
 
174
 
 
175
    END DO Loop_Phase1
 
176
 
 
177
    ! For the b.c's
 
178
    Loop_Surfaces: DO II = 1, STOTEL * CV_SNLOC
 
179
 
 
180
       Loop_Phase2: DO IPHASE = 1, NPHASE
 
181
          IS_IPHA = II + ( IPHASE - 1 ) * NPHASE
 
182
 
 
183
          IF( CP_OPTION( IPHASE ) == 1 ) THEN
 
184
             SUF_CPD_BCU( IS_IPHA ) = CP_COEFS( IPHASE, 1 ) * SUF_D_BCU( IS_IPHA )
 
185
          ELSE
 
186
             STOP 3932
 
187
          ENDIF
 
188
 
 
189
       END DO Loop_Phase2
 
190
    END DO Loop_Surfaces
 
191
 
 
192
    RETURN
 
193
 
 
194
  END SUBROUTINE CAL_CPDEN
 
195
 
 
196
 
 
197
 
 
198
  SUBROUTINE CAL_U_ABSORB( MAT_NONODS, CV_NONODS, NPHASE, NDIM, SATURA, TOTELE, CV_NLOC, MAT_NLOC, &
 
199
       CV_NDGLN, MAT_NDGLN, &
 
200
       NUABS_COEFS, UABS_COEFS, UABS_OPTION, U_ABSORB, VOLFRA_PORE, PERM, &
 
201
       X, X_NLOC, X_NONODS, X_NDGLN, &
 
202
       OPT_VEL_UPWIND_COEFS, NOPT_VEL_UPWIND_COEFS )
 
203
    ! Calculate absorption for momentum eqns
 
204
    use matrix_operations
 
205
    !    use cv_advection
 
206
    implicit none
 
207
    INTEGER, intent( in ) :: MAT_NONODS, CV_NONODS, NPHASE, NDIM, TOTELE, CV_NLOC,MAT_NLOC, &
 
208
         NUABS_COEFS, X_NLOC, X_NONODS, NOPT_VEL_UPWIND_COEFS 
 
209
    REAL, DIMENSION( CV_NONODS * NPHASE ), intent( in ) :: SATURA
 
210
    INTEGER, DIMENSION( TOTELE * CV_NLOC ), intent( in ) :: CV_NDGLN
 
211
    INTEGER, DIMENSION( TOTELE * X_NLOC ), intent( in ) :: X_NDGLN
 
212
    INTEGER, DIMENSION( TOTELE * MAT_NLOC ), intent( in ) :: MAT_NDGLN
 
213
    REAL, DIMENSION( NPHASE, NUABS_COEFS ), intent( in ) :: UABS_COEFS
 
214
    REAL, DIMENSION( X_NONODS ), intent( in ) :: X
 
215
    INTEGER, DIMENSION( NPHASE ), intent( in ) ::  UABS_OPTION
 
216
    REAL, DIMENSION( MAT_NONODS, NDIM * NPHASE, NDIM * NPHASE ), intent( inout ) :: U_ABSORB
 
217
    REAL, DIMENSION( TOTELE ), intent( in ) :: VOLFRA_PORE
 
218
    REAL, DIMENSION( TOTELE, NDIM, NDIM ), intent( in ) :: PERM
 
219
    REAL, DIMENSION( NOPT_VEL_UPWIND_COEFS ), intent( inout ) :: OPT_VEL_UPWIND_COEFS
 
220
    ! local variable...
 
221
    REAL, DIMENSION( :, :, :), allocatable :: U_ABSORB2
 
222
    REAL, DIMENSION( : ), allocatable :: SATURA2
 
223
    REAL :: PERT
 
224
    INTEGER :: ELE, IMAT, ICV, IPHASE, CV_ILOC, IDIM, JDIM, IJ
 
225
 
 
226
    write(357,*) 'In CAL_U_ABSORB'
 
227
 
 
228
    ALLOCATE( U_ABSORB2( MAT_NONODS, NDIM * NPHASE, NDIM * NPHASE ))
 
229
    ALLOCATE( SATURA2( CV_NONODS * NPHASE ))
 
230
    U_ABSORB2 = 0.
 
231
    SATURA2 = 0.
 
232
 
 
233
    write(357,*)'b4 in cal_u_absorb2, SATURA0',size(SATURA),SATURA
 
234
    CALL CAL_U_ABSORB2( MAT_NONODS, CV_NONODS, NPHASE, NDIM, SATURA, TOTELE, CV_NLOC, MAT_NLOC, &
 
235
         CV_NDGLN, MAT_NDGLN, &
 
236
         NUABS_COEFS, UABS_COEFS, UABS_OPTION, U_ABSORB, PERM, &
 
237
         X, X_NLOC, X_NONODS, X_NDGLN )
 
238
 
 
239
 
 
240
    PERT=0.0001
 
241
    SATURA2( 1 : CV_NONODS ) = SATURA( 1 : CV_NONODS ) + PERT
 
242
    IF (NPHASE > 1) SATURA2( 1 + CV_NONODS : 2 * CV_NONODS ) = SATURA( 1 + CV_NONODS : 2 * CV_NONODS ) - PERT
 
243
 
 
244
    !write(357,*)'b4 in cal_u_absorb2, MAT_NONODS, CV_NONODS, NPHASE, NDIM:',MAT_NONODS, CV_NONODS, NPHASE, NDIM
 
245
    !write(357,*)'b4 in cal_u_absorb2, TOTELE, CV_NLOC, MAT_NLOC, X_NLOC, X_NONODS:',TOTELE, CV_NLOC, MAT_NLOC, X_NLOC, X_NONODS
 
246
    !write(357,*)'b4 in cal_u_absorb2, SATURA',size(SATURA),SATURA
 
247
    !write(357,*)'b4 in cal_u_absorb2, SATURA2',size(SATURA2),SATURA2
 
248
    !write(357,*)'b4 in cal_u_absorb2, CV_NDGLN',size(CV_NDGLN),CV_NDGLN
 
249
    !write(357,*)'b4 in cal_u_absorb2, MAT_NDGLN',size(MAT_NDGLN),MAT_NDGLN
 
250
    !write(357,*)'b4 in cal_u_absorb2, UABS_COEFS',size(UABS_COEFS),UABS_COEFS
 
251
    !write(357,*)'b4 in cal_u_absorb2, UABS_OPTION',size(UABS_OPTION),UABS_OPTION
 
252
    !write(357,*)'b4 in cal_u_absorb2, VOLFRA_PORE',size(VOLFRA_PORE),VOLFRA_PORE
 
253
    !write(357,*)'b4 in cal_u_absorb2, PERM',size(PERM),PERM
 
254
    !write(357,*)'b4 in cal_u_absorb2, X',size(X),X
 
255
    !write(357,*)'b4 in cal_u_absorb2, X_NDGLN',size(X_NDGLN),X_NDGLN
 
256
    !write(357,*)'b4 in cal_u_absorb2, U_ABSORB2:', size(U_ABSORB2),U_ABSORB2
 
257
    !write(357,*)'b4 in cal_u_absorb2, OPT_VEL_UPWIND_COEFS:', size(OPT_VEL_UPWIND_COEFS),OPT_VEL_UPWIND_COEFS
 
258
 
 
259
    CALL CAL_U_ABSORB2( MAT_NONODS, CV_NONODS, NPHASE, NDIM, SATURA2, TOTELE, CV_NLOC, MAT_NLOC, &
 
260
         CV_NDGLN, MAT_NDGLN, &
 
261
         NUABS_COEFS, UABS_COEFS, UABS_OPTION, U_ABSORB2, PERM, &
 
262
         X, X_NLOC, X_NONODS, X_NDGLN )
 
263
 
 
264
    write(357,*)'after in cal_u_absorb, U_ABSORB2:', size(U_ABSORB2),U_ABSORB2
 
265
 
 
266
    DO ELE=1,TOTELE
 
267
       DO CV_ILOC=1,CV_NLOC
 
268
          IMAT=MAT_NDGLN((ELE-1)*MAT_NLOC+CV_ILOC)
 
269
          ICV=CV_NDGLN((ELE-1)*CV_NLOC+CV_ILOC)
 
270
          DO IPHASE=1,NPHASE
 
271
             DO IDIM=1,NDIM
 
272
                DO JDIM=1,NDIM
 
273
                   IJ=(IPHASE-1)*MAT_NONODS*NDIM*NDIM + (IMAT-1)*NDIM*NDIM + (IDIM-1)*NDIM +JDIM
 
274
                   OPT_VEL_UPWIND_COEFS(IJ) &
 
275
                        = U_ABSORB(IMAT, IDIM+(IPHASE-1)*NDIM, JDIM+(IPHASE-1)*NDIM) 
 
276
                   ! This is the gradient...
 
277
                   OPT_VEL_UPWIND_COEFS(IJ+NPHASE*MAT_NONODS*NDIM*NDIM) &
 
278
                        = (U_ABSORB2(IMAT, IDIM+(IPHASE-1)*NDIM, JDIM+(IPHASE-1)*NDIM) &
 
279
                        -U_ABSORB(IMAT, IDIM+(IPHASE-1)*NDIM, JDIM+(IPHASE-1)*NDIM))  &
 
280
                        / (SATURA2(ICV+(IPHASE-1)*CV_NONODS)-SATURA(ICV+(IPHASE-1)*CV_NONODS)) 
 
281
                END DO
 
282
             END DO
 
283
          END DO
 
284
       END DO
 
285
    END DO
 
286
 
 
287
    write(357,*)'in cal_u_absorb, U_ABSORB:', size(U_ABSORB),U_ABSORB
 
288
    write(357,*)'in cal_u_absorb, OPT_VEL_UPWIND_COEFS:', size(OPT_VEL_UPWIND_COEFS),OPT_VEL_UPWIND_COEFS
 
289
    write(357,*) 'Leaving CAL_U_ABSORB'
 
290
    RETURN
 
291
 
 
292
  END SUBROUTINE CAL_U_ABSORB
 
293
 
 
294
 
 
295
 
 
296
 
 
297
  SUBROUTINE CAL_U_ABSORB2( MAT_NONODS, CV_NONODS, NPHASE, NDIM, SATURA, TOTELE, CV_NLOC, MAT_NLOC, &
 
298
       CV_NDGLN, MAT_NDGLN, &
 
299
       NUABS_COEFS, UABS_COEFS, UABS_OPTION, U_ABSORB, PERM, &
 
300
       X, X_NLOC, X_NONODS, X_NDGLN ) 
 
301
    ! Calculate absorption for momentum eqns
 
302
    use matrix_operations
 
303
    !    use cv_advection
 
304
    implicit none
 
305
    INTEGER, intent( in ) :: MAT_NONODS, CV_NONODS, NPHASE, NDIM, TOTELE, CV_NLOC,MAT_NLOC, &
 
306
         NUABS_COEFS, X_NLOC, X_NONODS
 
307
    REAL, DIMENSION( CV_NONODS * NPHASE ), intent( in ) :: SATURA
 
308
    INTEGER, DIMENSION( TOTELE * CV_NLOC ), intent( in ) :: CV_NDGLN
 
309
    INTEGER, DIMENSION( TOTELE * X_NLOC ), intent( in ) :: X_NDGLN
 
310
    INTEGER, DIMENSION( TOTELE * MAT_NLOC ), intent( in ) :: MAT_NDGLN
 
311
    REAL, DIMENSION( NPHASE, NUABS_COEFS ), intent( in ) :: UABS_COEFS
 
312
    REAL, DIMENSION( X_NONODS ), intent( in ) :: X
 
313
    INTEGER, DIMENSION( NPHASE ), intent( in ) ::  UABS_OPTION
 
314
    REAL, DIMENSION( MAT_NONODS, NDIM * NPHASE, NDIM * NPHASE ), intent( inout ) :: U_ABSORB
 
315
    REAL, DIMENSION( TOTELE, NDIM, NDIM ), intent( in ) :: PERM
 
316
    ! Local variable
 
317
    REAL, PARAMETER :: TOLER = 1.E-10, REF_MOBILITY = 10.0
 
318
    INTEGER :: ELE, CV_ILOC, CV_NOD, CV_PHA_NOD, MAT_NOD, JPHA_JDIM, &
 
319
         IPHA_IDIM, IDIM, JDIM, IPHASE, II, jphase, X_NODI, MAT_NOD1, MAT_NOD2, MAT_NOD3, &
 
320
         CV_NOD1, CV_NOD2, ELE2, CV_PHA_NOD1, CV_PHA_NOD2, &
 
321
         MAT_NOD_MID, CV_ILOC2
 
322
    REAL :: SATURATION, ABS_SUM, RSUM, &
 
323
         U_ABS1, U_ABS2, U_ABS3
 
324
    REAL, DIMENSION( :, :, :), allocatable :: INV_PERM
 
325
    LOGICAL :: MEAN_ELE
 
326
 
 
327
    write(357,*) 'In CAL_U_ABSORB2'
 
328
 
 
329
    ALLOCATE( INV_PERM( TOTELE, NDIM, NDIM ))
 
330
 
 
331
    CALL PHA_BLOCK_INV( INV_PERM, PERM, TOTELE, NDIM )
 
332
 
 
333
    U_ABSORB = 0.0
 
334
 
 
335
    Loop_ELE: DO ELE = 1, TOTELE
 
336
 
 
337
       Loop_CVNLOC: DO CV_ILOC = 1, CV_NLOC
 
338
 
 
339
          MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
340
          CV_NOD = CV_NDGLN(( ELE - 1) * CV_NLOC + CV_ILOC )
 
341
 
 
342
          Loop_NPHASE: DO IPHASE = 1, NPHASE
 
343
             CV_PHA_NOD = CV_NOD + ( IPHASE - 1 ) * CV_NONODS
 
344
 
 
345
             Loop_DimensionsI: DO IDIM = 1, NDIM
 
346
 
 
347
                Loop_DimensionsJ: DO JDIM = 1, NDIM
 
348
 
 
349
                   SATURATION = SATURA( CV_PHA_NOD ) 
 
350
                   IPHA_IDIM = ( IPHASE - 1 ) * NDIM + IDIM 
 
351
                   JPHA_JDIM = ( IPHASE - 1 ) * NDIM + JDIM 
 
352
                   !JPHA_JDIM = IPHA_IDIM
 
353
 
 
354
                   Case_UABSOption: SELECT CASE( UABS_OPTION( IPHASE ) )
 
355
                   CASE( 0 ) 
 
356
                      ! no absorption option
 
357
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = 0.0
 
358
 
 
359
                   CASE( 1 )
 
360
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = INV_PERM( ELE, IDIM, JDIM )
 
361
 
 
362
                   CASE( 2 ) ! A standard polynomial representation of relative permeability                   
 
363
                      ABS_SUM = 0.
 
364
                      DO II = 1, NUABS_COEFS
 
365
                         ABS_SUM = ABS_SUM + UABS_COEFS( IPHASE, II) * SATURATION** (II - 1 )
 
366
                      END DO
 
367
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = INV_PERM( ELE, IDIM, JDIM ) / &
 
368
                           MAX( TOLER, ABS_SUM )
 
369
 
 
370
                   CASE( 3 ) ! From Tara and Martin's notes for relative permeability
 
371
 
 
372
                      CALL ABS3( U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ), &
 
373
                           INV_PERM( ELE, IDIM, JDIM ), max(0.0,SATURA(CV_NOD)), IPHASE)
 
374
 
 
375
                   CASE( 4 ) 
 
376
                      ABS_SUM = 0.0
 
377
                      DO II = 1, NUABS_COEFS
 
378
                         ABS_SUM = ABS_SUM + UABS_COEFS( IPHASE, II) * SATURATION** (II - 1 )
 
379
                      END DO
 
380
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = ABS_SUM
 
381
 
 
382
                   CASE DEFAULT
 
383
                      write(357,*) 'Unknown Option -- UABS_OPTION( IPHASE )'
 
384
                      ABS_SUM = 0.0
 
385
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = ABS_SUM
 
386
                      !                      stop 1023
 
387
 
 
388
                   END SELECT Case_UABSOption
 
389
 
 
390
 
 
391
                   !                   write(357,*)'ELE, IDIM, JDIM, iphase, U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ):', &
 
392
                   !                        ELE, IDIM, JDIM, iphase, U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )
 
393
                   !U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) = INV_PERM( ELE, IDIM, JDIM ) / &
 
394
                   !     MAX( TOLER, UABS_COEFS( IPHASE, 1 ) + UABS_COEFS( IPHASE, 2) * SATURATION + &
 
395
                   !     UABS_COEFS( IPHASE, 3 ) * SATURATION **2 )
 
396
                   !write(357,*) 'MAT_NOD, IPHA, JPHA, sat., U_ABSORB:', MAT_NOD, IPHA_IDIM, JPHA_JDIM, saturation, U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )
 
397
 
 
398
                END DO Loop_DimensionsJ
 
399
 
 
400
             END DO Loop_DimensionsI
 
401
 
 
402
          END DO Loop_NPHASE
 
403
 
 
404
       END DO Loop_CVNLOC
 
405
 
 
406
    END DO Loop_ELE
 
407
    !    stop 23
 
408
 
 
409
    Conditional_NotUsed1: IF(.false.) THEN 
 
410
       !    IF(.true.) THEN
 
411
       ! Use upwind absorptions. 
 
412
       !       IF((CV_NONODS==CV_NLOC*TOTELE).OR.(CV_NLOC==1)) THEN
 
413
       IF(.true.) THEN
 
414
          DO ELE = 2,TOTELE
 
415
             DO CV_ILOC=1,1
 
416
                MAT_NOD=MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
417
                CV_NOD=CV_NDGLN(( ELE - 1 ) * CV_NLOC + CV_ILOC)
 
418
 
 
419
                DO IPHA_IDIM=1,NDIM*NPHASE
 
420
                   DO JPHA_JDIM=1,NDIM*NPHASE
 
421
                      if(IPHA_IDIM==JPHA_JDIM) then
 
422
                         CALL ABS3(rsum, &
 
423
                              1.0, 0.5*(SATURA(CV_NOD)+SATURA(CV_NOD-1)), IPHA_IDIM)
 
424
                         U_ABSORB( MAT_NOD,   IPHA_IDIM, JPHA_JDIM )=RSUM
 
425
                         U_ABSORB( MAT_NOD-1, IPHA_IDIM, JPHA_JDIM )=RSUM
 
426
                      endif
 
427
                   END DO
 
428
                END DO
 
429
             END DO
 
430
          END DO
 
431
       ELSE
 
432
          DO ELE = TOTELE,1,-1
 
433
             DO CV_ILOC=CV_NLOC,1,-1
 
434
                MAT_NOD=MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
435
                ELE2=ELE
 
436
                CV_ILOC2=CV_ILOC-1
 
437
                IF(CV_ILOC==1) THEN
 
438
                   CV_ILOC2=CV_NLOC-1
 
439
                   ELE2=ELE-1
 
440
                   IF(ELE2==0) THEN
 
441
                      CV_ILOC2=CV_ILOC
 
442
                      ELE2=ELE
 
443
                   ENDIF
 
444
                ENDIF
 
445
                MAT_NOD2=MAT_NDGLN(( ELE2 - 1 ) * MAT_NLOC + CV_ILOC2)
 
446
                DO IPHA_IDIM=1,NDIM*NPHASE
 
447
                   DO JPHA_JDIM=1,NDIM*NPHASE
 
448
                      U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM )
 
449
                   END DO
 
450
                END DO
 
451
             END DO
 
452
          END DO
 
453
       ENDIF
 
454
    ENDIF Conditional_NotUsed1
 
455
 
 
456
    Conditional_NotUsed2: IF(.FALSE.) THEN
 
457
       DO ELE = 1, TOTELE
 
458
          DO IPHA_IDIM=1,NDIM*NPHASE
 
459
             DO JPHA_JDIM=1,NDIM*NPHASE
 
460
                RSUM=0.0
 
461
                !          RSUM=1.e+10
 
462
                DO CV_ILOC = 1, CV_NLOC
 
463
                   MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
464
                   X_NODI = X_NDGLN(( ELE - 1 ) * X_NLOC + CV_ILOC)
 
465
                   U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=X(X_NODI)
 
466
                END DO
 
467
             END DO
 
468
          END DO
 
469
       END DO
 
470
    ENDIF Conditional_NotUsed2
 
471
 
 
472
    MEAN_ELE=.false.
 
473
    IF( MEAN_ELE ) THEN
 
474
       ! Average absoption throughout element -using the harmonic average
 
475
       DO ELE = 1, TOTELE
 
476
          DO IPHA_IDIM=1,NDIM*NPHASE
 
477
             DO JPHA_JDIM=1,NDIM*NPHASE
 
478
                RSUM=0.0
 
479
                !          RSUM=1.e+10
 
480
                DO CV_ILOC = 1, CV_NLOC
 
481
                   MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
482
                   !         RSUM=RSUM+ (1./max(1.e-10,U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )))/REAL(CV_NLOC+1)
 
483
                   !         IF(CV_ILOC==2) RSUM=RSUM+ (1./max(1.e-10,U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )))/REAL(CV_NLOC+1)
 
484
                   !         RSUM=RSUM+ U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )/REAL(CV_NLOC)
 
485
                   RSUM=RSUM+ U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )/REAL(CV_NLOC+1)
 
486
                   IF(CV_ILOC==2) RSUM=RSUM+ U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )/REAL(CV_NLOC+1)
 
487
                   !         RSUM=max(RSUM, U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ))
 
488
                   !         RSUM=min(RSUM, U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ))
 
489
                END DO
 
490
                MAT_NOD1 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 1)
 
491
                MAT_NOD2 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 2)
 
492
                MAT_NOD3 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 3)
 
493
 
 
494
                U_ABS1= &
 
495
                     0.66666* U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ) &
 
496
                     +0.33333* U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM )
 
497
 
 
498
                U_ABS2= &
 
499
                     0.5*0.5* U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ) &
 
500
                     +0.5* U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ) &
 
501
                     +0.5*0.5* U_ABSORB( MAT_NOD3, IPHA_IDIM, JPHA_JDIM ) 
 
502
 
 
503
                U_ABS3= &
 
504
                     0.33333* U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ) &
 
505
                     +0.66666* U_ABSORB( MAT_NOD3, IPHA_IDIM, JPHA_JDIM )
 
506
 
 
507
                U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM )=U_ABS1
 
508
                U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM )=U_ABS2
 
509
                U_ABSORB( MAT_NOD3, IPHA_IDIM, JPHA_JDIM )=U_ABS3
 
510
 
 
511
                DO CV_ILOC = 1, -CV_NLOC
 
512
                   MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
513
                   MAT_NOD1 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 3)
 
514
                   !         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=  &
 
515
                   !       0.25*U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM )  &
 
516
                   !       +0.25*U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )
 
517
                   U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=RSUM
 
518
                   !         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=0.65*1./RSUM &
 
519
                   !                 + 0.35*U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )
 
520
                   !         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=0.5*1./RSUM &
 
521
                   !                 + 0.5*U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )
 
522
                   !         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )=1./RSUM 
 
523
                END DO
 
524
             END DO
 
525
          END DO
 
526
       END DO
 
527
 
 
528
    ENDIF
 
529
 
 
530
    !    write(357,*)'CV_NONODs,TOTELE*CV_NLOC,CV_NLOC:',CV_NONODs,TOTELE*CV_NLOC,CV_NLOC
 
531
    !     stop 383
 
532
 
 
533
    !    IF(CV_NONODS/=TOTELE*CV_NLOC) THEN
 
534
    Conditional_NotUsed3: IF(.false.) THEN
 
535
       IF(CV_NLOC==3) THEN
 
536
          !    IF(.false.) THEN
 
537
          DO ELE = 1, -TOTELE
 
538
             DO IPHASE = 1, NPHASE
 
539
                DO IDIM = 1, NDIM
 
540
                   DO JDIM = 1, NDIM
 
541
                      IPHA_IDIM = ( IPHASE - 1 ) * NDIM + IDIM 
 
542
                      JPHA_JDIM = ( IPHASE - 1 ) * NDIM + JDIM 
 
543
                      DO CV_ILOC=1,CV_NLOC
 
544
                         MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
545
                         CV_NOD = CV_NDGLN(( ELE - 1 ) * CV_NLOC + CV_ILOC)
 
546
                         SATURATION=SATURA(CV_NOD + ( IPHASE - 1 ) * CV_NONODS )
 
547
                         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )  &
 
548
                              = U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) &
 
549
                              / max(0.03,SATURATION)
 
550
                      END DO
 
551
                   END DO
 
552
                END DO
 
553
             END DO
 
554
          END DO
 
555
          DO ELE = 1, TOTELE
 
556
             DO IPHASE = 1, NPHASE
 
557
                DO IDIM = 1, NDIM
 
558
                   DO JDIM = 1, NDIM
 
559
                      IPHA_IDIM = ( IPHASE - 1 ) * NDIM + IDIM 
 
560
                      JPHA_JDIM = ( IPHASE - 1 ) * NDIM + JDIM 
 
561
                      MAT_NOD1 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + MAT_NLOC)
 
562
                      MAT_NOD2 = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 1)
 
563
                      MAT_NOD_MID = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + 2)
 
564
 
 
565
                      CV_NOD1 = CV_NDGLN(( ELE - 1 ) * CV_NLOC + CV_NLOC)
 
566
                      CV_NOD2 = CV_NDGLN(( ELE - 1 ) * CV_NLOC + 1)
 
567
 
 
568
                      CV_PHA_NOD1 = CV_NOD1 + ( IPHASE - 1 ) * CV_NONODS
 
569
                      CV_PHA_NOD2 = CV_NOD2 + ( IPHASE - 1 ) * CV_NONODS
 
570
                      if(.true.) then
 
571
                         !           stop 56
 
572
                         ! max...
 
573
                         IF(U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )  &
 
574
                              >U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ) ) THEN
 
575
                            IF(U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )  &
 
576
                                 >U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ) ) THEN
 
577
                               U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )= &
 
578
                                    MAX(U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ), &
 
579
                                    U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ))
 
580
                            ENDIF
 
581
                         ENDIF
 
582
                         ! min...
 
583
                         IF(U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )  &
 
584
                              <U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ) ) THEN
 
585
                            IF(U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )  &
 
586
                                 <U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ) ) THEN
 
587
                               U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )= &
 
588
                                    MIN(U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ), &
 
589
                                    U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM ))
 
590
                            ENDIF
 
591
                         ENDIF
 
592
                      ENDIF
 
593
                      if(.false.) then
 
594
                         U_ABSORB( MAT_NOD_MID, IPHA_IDIM, JPHA_JDIM )= &
 
595
                              0.5*U_ABSORB( MAT_NOD1, IPHA_IDIM, JPHA_JDIM ) &
 
596
                              +0.5*U_ABSORB( MAT_NOD2, IPHA_IDIM, JPHA_JDIM )
 
597
                      ENDIF
 
598
                   END DO
 
599
                END DO
 
600
             END DO
 
601
          END DO
 
602
          DO ELE = 1, -TOTELE
 
603
             DO IPHASE = 1, NPHASE
 
604
                DO IDIM = 1, NDIM
 
605
                   DO JDIM = 1, NDIM
 
606
                      IPHA_IDIM = ( IPHASE - 1 ) * NDIM + IDIM 
 
607
                      JPHA_JDIM = ( IPHASE - 1 ) * NDIM + JDIM 
 
608
                      DO CV_ILOC=1,CV_NLOC
 
609
                         MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + CV_ILOC)
 
610
                         CV_NOD = CV_NDGLN(( ELE - 1 ) * CV_NLOC + CV_ILOC)
 
611
                         SATURATION=SATURA(CV_NOD + ( IPHASE - 1 ) * CV_NONODS )
 
612
                         U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM )  &
 
613
                              = U_ABSORB( MAT_NOD, IPHA_IDIM, JPHA_JDIM ) &
 
614
                              * max(0.03,SATURATION)
 
615
                      END DO
 
616
                   END DO
 
617
                END DO
 
618
             END DO
 
619
          END DO
 
620
       ENDIF
 
621
    ENDIF Conditional_NotUsed3
 
622
 
 
623
    write(357,*)'NUABS_COEFS,UABS_COEFS:',NUABS_COEFS,UABS_COEFS
 
624
    write(357,*)'U_ABSORB:'
 
625
    do mat_nod =1, -mat_nonods
 
626
       do iphase=1,nphase
 
627
          write(357,*)mat_nod,iphase,(U_ABSORB( MAT_NOD, iphase, jphase ),jphase=1,2) 
 
628
       end do
 
629
    end do
 
630
    !    write(357,*)'U_ABSORB:',U_ABSORB
 
631
    !    stop 4941
 
632
    DEALLOCATE( INV_PERM )
 
633
 
 
634
    write(357,*) 'Leaving CAL_U_ABSORB2'
 
635
 
 
636
    RETURN
 
637
 
 
638
  END SUBROUTINE CAL_U_ABSORB2
 
639
 
 
640
 
 
641
 
 
642
 
 
643
  SUBROUTINE ABS3( ABSP, INV_PERM, SAT, IPHASE )
 
644
    IMPLICIT NONE
 
645
    REAL, intent( inout ) :: ABSP
 
646
    REAL, intent( in ) :: SAT, INV_PERM
 
647
    INTEGER, intent( in ) :: IPHASE
 
648
    ! Local variables...
 
649
    REAL :: VISC1, VISC2, S_GC, S_OR, REF_MOBILITY, &
 
650
         KR1, KR2, KR, VISC, SATURATION, ABS_SUM, SAT2
 
651
 
 
652
    VISC1 = 1.0
 
653
    S_GC = 0.1
 
654
    S_OR = 0.3
 
655
    REF_MOBILITY = 10.0
 
656
 
 
657
    SATURATION = SAT
 
658
    IF( IPHASE == 2 ) SATURATION = 1. - SAT
 
659
 
 
660
    IF( SAT < S_GC ) THEN
 
661
       KR1 = 0.0
 
662
    ELSE IF( SAT > 1. -S_OR ) THEN
 
663
       KR1 = 1.0
 
664
    ELSE
 
665
       KR1 = ( ( SAT - S_GC) / ( 1. - S_GC - S_OR ))**2
 
666
    ENDIF
 
667
 
 
668
    SAT2 = 1.0 - SAT
 
669
    IF( SAT2 < S_OR ) THEN
 
670
       KR2 = 0.0
 
671
    ELSEIF( SAT2 > 1. - S_GC ) THEN
 
672
       KR2 = 1.0
 
673
    ELSE
 
674
       KR2 = ( ( SAT2 - S_OR ) / ( 1. - S_GC - S_OR ))**2
 
675
    ENDIF
 
676
    VISC2 = REF_MOBILITY
 
677
 
 
678
    IF( IPHASE == 1 ) THEN
 
679
       KR = KR1
 
680
       VISC = VISC1
 
681
    ELSE
 
682
       KR = KR2
 
683
       VISC = VISC2
 
684
    ENDIF
 
685
 
 
686
    ABS_SUM = KR / MAX( 1.e-6, VISC * max( 0.01, SATURATION ))
 
687
 
 
688
    ABSP = INV_PERM / MAX( 1.e-6, ABS_SUM )
 
689
 
 
690
    if( iphase == 1 ) then
 
691
       ABSP =  min( 1.e+4, ABSP )
 
692
       if( saturation < 0.1 ) then
 
693
          ABSP = ( 1. + max( 100. * ( 0.1 - saturation ), 0.0 )) * ABSP
 
694
       endif
 
695
    else
 
696
       ABSP = min( 1.e+5, ABSP )
 
697
       if( saturation < 0.3 ) then
 
698
          ABSP = ( 1. + max( 100. * ( 0.3 - saturation ), 0.0 )) * ABSP
 
699
       endif
 
700
    endif
 
701
 
 
702
    RETURN
 
703
  END SUBROUTINE ABS3
 
704
 
 
705
 
 
706
 
 
707
  SUBROUTINE CALC_CAPIL_PRES( CAPIL_PRES_OPT, CV_NONODS, NPHASE, NCAPIL_PRES_COEF, &
 
708
       CAPIL_PRES_COEF, IPLIKE_GRAD_SOU, PLIKE_GRAD_SOU_COEF, PLIKE_GRAD_SOU_GRAD, SATURA )
 
709
 
 
710
    ! CAPIL_PRES_OPT is the capillary pressure option for deciding what form it might take. 
 
711
    ! CAPIL_PRES_COEF( NCAPIL_PRES_COEF, NPHASE, NPHASE ) are the coefficients 
 
712
    ! Capillary pressure coefs have the dims CAPIL_PRES_COEF( NCAPIL_PRES_COEF, NPHASE,NPHASE )
 
713
    ! used to calulcate the capillary pressure. 
 
714
 
 
715
    IMPLICIT NONE
 
716
    INTEGER, intent( in ) :: IPLIKE_GRAD_SOU, CAPIL_PRES_OPT, CV_NONODS, NPHASE, NCAPIL_PRES_COEF
 
717
    REAL, DIMENSION( NCAPIL_PRES_COEF, NPHASE, NPHASE ), intent( in ) :: CAPIL_PRES_COEF
 
718
    REAL, DIMENSION( IPLIKE_GRAD_SOU*CV_NONODS * NPHASE ), intent( inout ) :: &
 
719
         PLIKE_GRAD_SOU_COEF, PLIKE_GRAD_SOU_GRAD
 
720
    REAL, DIMENSION( CV_NONODS * NPHASE ), intent( in ) :: SATURA
 
721
    ! Local Variables
 
722
    INTEGER :: IPHASE, JPHASE
 
723
 
 
724
 
 
725
    IF( IPLIKE_GRAD_SOU /= 1 ) THEN
 
726
       WRITE(357,*),'PROBLEM WITH THE CAPILLARY OPTION SET UP' 
 
727
       STOP 3831
 
728
    ENDIF
 
729
 
 
730
 
 
731
    Case_CapillaryPressure: SELECT CASE( CAPIL_PRES_OPT )
 
732
 
 
733
    CASE DEFAULT
 
734
 
 
735
       WRITE(357,*),'PROBLEM WITH THE CAPILLARY OPTION SET UP - OPTION NOT AVAILABLE' 
 
736
       STOP 3832
 
737
 
 
738
    CASE( 1 )
 
739
       ! Capillary pressure coefs 
 
740
       ! In this form PLIKE_GRAD_SOU_GRAD is the actual capillary pressure. 
 
741
       PLIKE_GRAD_SOU_COEF = 1.0
 
742
       PLIKE_GRAD_SOU_GRAD = 0.0
 
743
 
 
744
       DO IPHASE = 1, NPHASE
 
745
 
 
746
          DO JPHASE = 1, NPHASE
 
747
 
 
748
             PLIKE_GRAD_SOU_GRAD( 1 + ( IPHASE - 1 ) * CV_NONODS : IPHASE * CV_NONODS ) = &
 
749
                  PLIKE_GRAD_SOU_GRAD( 1 + ( IPHASE - 1 ) * CV_NONODS : IPHASE * CV_NONODS ) + &
 
750
                  CAPIL_PRES_COEF( 1, IPHASE, JPHASE ) * &
 
751
                  MAX( SATURA( 1 + ( JPHASE - 1 ) * CV_NONODS : JPHASE * CV_NONODS ), 0.0 ) &
 
752
                  ** CAPIL_PRES_COEF( 2, IPHASE, JPHASE )
 
753
 
 
754
          END DO
 
755
 
 
756
       END DO
 
757
 
 
758
    end SELECT Case_CapillaryPressure
 
759
 
 
760
    RETURN
 
761
  END SUBROUTINE CALC_CAPIL_PRES
 
762
 
 
763
 
 
764
end module multiphase_EOS