1
! Copyright (C) 2006 Imperial College London and others.
3
! Please see the AUTHORS file in the main source directory for a full list
4
! of copyright holders.
7
! Applied Modelling and Computation Group
8
! Department of Earth Science and Engineering
9
! Imperial College London
11
! amcgsoftware@imperial.ac.uk
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.
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.
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
34
subroutine open_output_file(output_channel,output_name,output_name_length,file_format)
36
integer, intent(in) :: output_channel,output_name_length
37
character(len=12), intent(in) :: file_format
38
character(len=output_name_length), intent(in) :: output_name
43
write(357,*) 'In open_output_file'
45
open(output_channel,file=trim(output_name),status='replace',form=trim(file_format),action='write',iostat=ierror)
47
if (ierror .ne. 0) then
48
write(*,*) 'Problem opening output file ',trim(output_name)
50
end if ! if (ierror .ne. 0) then
52
write(357,*) 'Leaving open_output_file'
54
end subroutine open_output_file
56
!----------------------------------------------------------------
58
subroutine write_integer_to_string(integer_variable,string_variable,len_string_variable)
60
! write an integer to a string
62
integer, intent(in) :: integer_variable,len_string_variable
63
character(len=len_string_variable), intent(inout) :: string_variable
65
write(357,*) 'In write_integer_to_string'
67
if (integer_variable .lt. 10) then
69
if (len_string_variable .lt. 1) then
71
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
72
len_string_variable,integer_variable
75
end if ! if (len_string_variable .lt. 1) then
77
write(unit=string_variable,fmt='(I1)') integer_variable
79
else if ((integer_variable .ge. 10) .and. (integer_variable .lt. 100)) then
81
if (len_string_variable .lt. 2) then
83
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
84
len_string_variable,integer_variable
87
end if ! if (len_string_variable .lt. 2) then
89
write(unit=string_variable,fmt='(I2)') integer_variable
91
else if ((integer_variable .ge. 100) .and. (integer_variable .lt. 1000)) then
93
if (len_string_variable .lt. 3) then
95
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
96
len_string_variable,integer_variable
99
end if ! if (len_string_variable .lt. 3) then
101
write(unit=string_variable,fmt='(I3)') integer_variable
103
else if ((integer_variable .ge. 1000) .and. (integer_variable .lt. 10000)) then
105
if (len_string_variable .lt. 4) then
107
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
108
len_string_variable,integer_variable
111
end if ! if (len_string_variable .lt. 4) then
113
write(unit=string_variable,fmt='(I4)') integer_variable
115
else if ((integer_variable .ge. 10000) .and. (integer_variable .lt. 100000)) then
117
if (len_string_variable .lt. 5) then
119
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
120
len_string_variable,integer_variable
123
end if ! if (len_string_variable .lt. 5) then
125
write(unit=string_variable,fmt='(I5)') integer_variable
127
else if ((integer_variable .ge. 100000) .and. (integer_variable .lt. 1000000)) then
129
if (len_string_variable .lt. 6) then
131
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
132
len_string_variable,integer_variable
135
end if ! if (len_string_variable .lt. 6) then
137
write(unit=string_variable,fmt='(I6)') integer_variable
139
else if ((integer_variable .ge. 1000000) .and. (integer_variable .lt. 10000000)) then
141
if (len_string_variable .lt. 7) then
143
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
144
len_string_variable,integer_variable
147
end if ! if (len_string_variable .lt. 7) then
149
write(unit=string_variable,fmt='(I7)') integer_variable
151
else if ((integer_variable .ge. 10000000) .and. (integer_variable .lt. 100000000)) then
153
if (len_string_variable .lt. 8) then
155
write(*,*) 'Write integer to string error - string to short,len_string_variable,integer_variable', &
156
len_string_variable,integer_variable
159
end if ! if (len_string_variable .lt. 8) then
161
write(unit=string_variable,fmt='(I8)') integer_variable
165
write(*,*) 'Write integer to string error - integer to big,integer_variable',integer_variable
168
end if ! if (integer_variable .lt. 10) then
170
write(357,*) 'Leaving write_integer_to_string'
172
end subroutine write_integer_to_string
174
! -----------------------------------------------------------------------------------------
176
subroutine output_fem_sol_of_cv( unit, totele, cv_nonods, x_nonods, nphase, &
177
cv_nloc, x_nloc, cv_ndgln, x_ndgln, x, femt )
181
integer, intent( in ) :: unit, totele, cv_nonods, x_nonods, nphase, cv_nloc, x_nloc
182
integer, dimension( totele * cv_nloc ), intent( in ) :: cv_ndgln
183
integer, dimension( totele * x_nloc ), intent( in ) :: x_ndgln
184
real, dimension( x_nonods ), intent( in ) :: x
185
real, dimension( cv_nonods * nphase ), intent( in ) :: femt
188
integer :: ele, cv_iloc
190
write(357,*) 'In output_fem_sol_of_cv'
193
do cv_iloc = 1, cv_nloc
194
write( unit ,* ) x( x_ndgln(( ele - 1 ) * x_nloc + cv_iloc )), &
195
femt( cv_ndgln(( ele - 1 ) * cv_nloc + cv_iloc ))
199
write(357,*) 'Leaving output_fem_sol_of_cv'
201
end subroutine output_fem_sol_of_cv
203
! -----------------------------------------------------------------------------------------
205
subroutine generate_name_dump( itime, unit, field, field_no )
207
integer, intent( in ) :: itime
208
integer, intent( inout ) :: unit
209
character( len = 50 ), intent( in ) :: field
210
integer, intent( in ) :: field_no
213
character( len = 50 ) :: file_name_in, file_name_out, dump
214
integer :: iaux, k, k1, k2, k3
216
write(357,*) 'In generate_name_dump'
219
file_name_in = 'test_'
220
k1 = index( file_name_in, ' ' ) - 1
221
k2 = index( field, ' ' ) - 1
223
k3 = index( dump, ' ' ) - 1
224
file_name_in = file_name_in( 1 : k1 ) // field( 1 : k2 ) // dump( 1 : k3 )
225
!file_name_in = trim( trim( file_name_in )//trim( field )//trim( dump )
227
open( iaux, file = 'tempfile', status = 'unknown' )
228
k = index( file_name_in, ' ' ) - 1
229
write( iaux, 222 ) file_name_in( 1 : k ), itime
230
!trim( file_name_in ), '.d.', itime
232
rewind( unit = iaux )
233
read( iaux, * ) file_name_out
234
k = index( file_name_out, ' ' ) - 1
235
!unit = itime + field_no
236
unit = 100000 + field_no
237
open( unit, file = trim( file_name_out ), status = 'unknown' )
243
!222 format( a, a3, i0 )
245
write(357,*) 'Leaving generate_name_dump'
247
end subroutine generate_name_dump
249
! -----------------------------------------------------------------------------------------
251
subroutine printing_field_array( unit, totele, &
252
cv_nonods, x_nonods, x_nloc, x_ndgln, cv_nloc, cv_ndgln, &
253
pos_x, field_length, field, iphase )
255
integer, intent( in ) :: unit, totele, x_nloc, cv_nonods, x_nonods, cv_nloc, iphase
256
integer, dimension( totele * x_nloc ), intent( in ) :: x_ndgln
257
integer, dimension( totele * cv_nloc ), intent( in ) :: cv_ndgln
258
real, dimension( x_nonods ), intent( in ) :: pos_x
259
integer, intent( in ) :: field_length
260
real, dimension( field_length ), intent( in ) :: field
263
integer :: cv_iloc, ele, xi_nod, xi_nod_plus, xi_nod_minus, field_nod
266
write(357,*) 'In printing_field_array'
268
Loop_Elements: do ele = 1, totele
271
xi_nod = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc )
272
xi_nod_plus = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc + 1 )
273
field_nod = cv_ndgln(( ele - 1 ) * cv_nloc + cv_iloc ) + ( iphase - 1 ) * cv_nonods
274
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_plus )
275
x_coord = 0.5 * x_coord
276
write( unit, * ) pos_x( xi_nod ), field( field_nod )
277
write( unit, * ) x_coord, field( field_nod )
279
do cv_iloc = 2, cv_nloc - 1
281
xi_nod = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc )
282
xi_nod_plus = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc + 1 )
283
xi_nod_minus = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc - 1 )
284
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_minus )
285
x_coord = 0.5 * x_coord
286
field_nod = cv_ndgln(( ele - 1 ) * cv_nloc + cv_iloc ) + ( iphase - 1 ) * cv_nonods
287
write( unit, * ) x_coord, field( field_nod )
288
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_plus )
289
x_coord = 0.5 * x_coord
290
write( unit, * ) x_coord, field( field_nod )
295
xi_nod = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc )
296
xi_nod_minus = x_ndgln(( ele - 1 ) * x_nloc + cv_iloc - 1 )
297
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_minus )
298
x_coord = 0.5 * x_coord
299
field_nod = cv_ndgln(( ele - 1 ) * cv_nloc + cv_iloc ) + ( iphase - 1 ) * cv_nonods
300
write( unit, * ) x_coord, field( field_nod )
301
write( unit, * ) pos_x( xi_nod ), field( field_nod )
305
write(357,*) 'Leaving printing_field_array'
307
end subroutine printing_field_array
309
! -----------------------------------------------------------------------------------------
311
subroutine printing_veloc_field( unit, totele, &
312
xu_nonods, xu_nloc, xu_ndgln, u_nloc, u_ndgln, &
313
pos_x, u_nonods, field_length, field, iphase )
315
integer, intent( in ) :: unit, totele, xu_nonods, xu_nloc
316
integer, dimension( totele * xu_nloc ), intent( in ) :: xu_ndgln
317
integer, intent( in ) :: u_nloc
318
integer, dimension( totele * u_nloc ), intent( in ) :: u_ndgln
319
real, dimension( xu_nonods ), intent( in ) :: pos_x
320
integer, intent( in ) :: u_nonods
321
integer, intent( in ) :: field_length
322
real, dimension( field_length ), intent( in ) :: field
323
integer, intent( in ) :: iphase
326
integer :: x_iloc, ele, xi_nod, xi_nod_plus, &
327
xi_nod_minus, field_nod
330
write(357,*) 'In printing_veloc_field'
332
Loop_Elements: do ele = 1, totele
335
xi_nod = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc )
336
xi_nod_plus = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc + 1 )
337
field_nod = u_ndgln(( ele - 1 ) * u_nloc + x_iloc ) + ( iphase - 1 ) * u_nonods
338
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_plus )
339
x_coord = 0.5 * x_coord
340
write( unit, * ) pos_x( xi_nod ), field( field_nod )
341
write( unit, * ) x_coord, field( field_nod )
343
do x_iloc = 2, xu_nloc - 1
345
xi_nod = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc )
346
xi_nod_plus = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc + 1 )
347
xi_nod_minus = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc - 1 )
348
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_minus )
349
x_coord = 0.5 * x_coord
350
field_nod = u_ndgln(( ele - 1 ) * u_nloc + x_iloc ) + ( iphase - 1 ) * u_nonods ! Is this correct? check later!
351
write( unit, * ) x_coord, field( field_nod )
352
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_plus )
353
x_coord = 0.5 * x_coord
354
write( unit, * ) x_coord, field( field_nod )
359
xi_nod = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc )
360
xi_nod_minus = xu_ndgln(( ele - 1 ) * xu_nloc + x_iloc - 1 )
361
x_coord = pos_x( xi_nod ) + pos_x( xi_nod_minus )
362
x_coord = 0.5 * x_coord
363
field_nod = u_ndgln(( ele - 1 ) * u_nloc + x_iloc ) + ( iphase - 1 ) * u_nonods
364
write( unit, * ) x_coord, field( field_nod )
365
write( unit, * ) pos_x( xi_nod ), field( field_nod )
369
write(357,*) 'Leaving printing_veloc_field'
371
end subroutine printing_veloc_field
373
! -----------------------------------------------------------------------------------------
375
subroutine printing_fw_field( unit, totele, &
376
cv_nonods, cv_nloc, cv_ndgln, &
377
field_length, field )
379
integer, intent( in ) :: unit, totele, cv_nonods, cv_nloc
380
integer, dimension( totele * cv_nloc ), intent( in ) :: cv_ndgln
381
integer, intent( in ) :: field_length
382
real, dimension( field_length ), intent( in ) :: field
385
integer :: cv_iloc, ele, cv_nod
386
real :: visc1, visc2, s_gc, s_or, fw, kr1, kr2, sat
388
write(357,*) 'In printing_fw_field'
390
Loop_ELE: DO ELE = 1, TOTELE
392
Loop_CVNLOC: DO CV_ILOC = 1, CV_NLOC
394
CV_NOD = CV_NDGLN(( ELE - 1) * CV_NLOC + CV_ILOC )
398
S_GC = 0.2 ! here S_gc --> S_wi
401
SAT = max( 0.0, FIELD( CV_NOD + CV_NONODS )) ! Second phase
402
KR2 = ( 1. - ( 1. - SAT - S_GC ) / MAX( 1.E-5, 1. - S_GC - S_OR )) **2
403
KR2 = KR2 / REAL( CV_NLOC )
405
SAT = MAX( 0.0, FIELD( CV_NOD )) ! First phase
406
KR1 = 0.7 * (( SAT - S_GC ) / MAX( 1.E-5, 1. - S_GC - S_OR )) **2
407
KR1 = KR1 / REAL( CV_NLOC )
409
FW = FW + 1. / ( 1. + VISC1 / KR1 * KR2 / VISC2 )
410
FW = FW / REAL( CV_NLOC )
414
write( unit, * ) SAT, KR1
418
write(357,*) 'Leaving printing_fw_field'
420
end subroutine printing_fw_field
422
! -----------------------------------------------------------------------------------------
424
subroutine check_sparsity( &
425
u_pha_nonods, cv_pha_nonods, &
426
u_nonods, cv_nonods, totele, &
427
mx_ncolacv, ncolacv, finacv, colacv, midacv, & ! CV multi-phase eqns (e.g. vol frac, temp)
428
nlenmcy, mx_ncolmcy, ncolmcy, finmcy, colmcy, midmcy, & ! Force balance plus cty multi-phase eqns
429
mxnele, ncolele, midele, finele, colele, & ! Element connectivity
430
mx_ncoldgm_pha, ncoldgm_pha, coldgm_pha, findgm_pha, middgm_pha, & ! Force balance sparsity
431
mx_nct, ncolct, findct, colct, & ! CT sparsity - global cty eqn
432
mx_nc, ncolc, findc, colc, & ! C sparsity operating on pressure in force balance
433
mx_ncolcmc, ncolcmc, findcmc, colcmc, midcmc, & ! pressure matrix for projection method
434
mx_ncolm, ncolm, findm, colm, midm )
438
integer, intent( in ) :: u_pha_nonods, cv_pha_nonods, u_nonods, cv_nonods, totele
439
integer, intent ( in ) :: mx_ncolacv, ncolacv
440
integer, dimension( cv_pha_nonods + 1 ), intent (in ) :: finacv
441
integer, dimension( mx_ncolacv ), intent (in ) :: colacv
442
integer, dimension( cv_pha_nonods ), intent (in ) :: midacv
443
integer, intent ( in ) :: nlenmcy, mx_ncolmcy, ncolmcy
444
integer, dimension( nlenmcy + 1 ), intent (in ) :: finmcy
445
integer, dimension( mx_ncolmcy ), intent (in ) :: colmcy
446
integer, dimension( nlenmcy ), intent (in ) :: midmcy
447
integer, intent ( in ) :: mxnele, ncolele
448
integer, dimension( totele ), intent (in ) :: midele
449
integer, dimension( totele + 1 ), intent (in ) :: finele
450
integer, dimension( mxnele ), intent (in ) :: colele
451
integer, intent ( in ) :: mx_ncoldgm_pha, ncoldgm_pha
452
integer, dimension( mx_ncoldgm_pha ), intent (in ) :: coldgm_pha
453
integer, dimension( u_pha_nonods + 1 ), intent (in ) :: findgm_pha
454
integer, dimension( u_pha_nonods ), intent (in ) :: middgm_pha
455
integer, intent ( in ) :: mx_nct, ncolct
456
integer, dimension( cv_nonods + 1 ), intent (in ) :: findct
457
integer, dimension( mx_nct ), intent (in ) :: colct
458
integer, intent ( in ) :: mx_nc, ncolc
459
integer, dimension( u_nonods + 1 ), intent (in ) :: findc
460
integer, dimension( mx_nc ), intent (in ) :: colc
461
integer, intent ( in ) :: mx_ncolcmc, ncolcmc
462
integer, dimension( cv_nonods + 1 ), intent (in ) :: findcmc
463
integer, dimension( mx_ncolcmc ), intent (in ) :: colcmc
464
integer, dimension( cv_nonods ), intent (in ) :: midcmc
465
integer, intent ( in ) :: mx_ncolm, ncolm
466
integer, dimension( cv_nonods + 1 ), intent (in ) :: findm
467
integer, dimension( ncolm ), intent (in ) :: colm
468
integer, dimension( cv_nonods ), intent (in ) :: midm
471
integer, dimension( : ), allocatable :: dummy
473
write(357,*) 'In check_sparsity'
475
open( 15, file = 'CheckSparsityMatrix.dat', status = 'unknown' )
476
write( 15, * )'########## FINMCY, MIDMCY, COLMCY ##################'
477
write(15, * )'NCOLMCY:', NCOLMCY
478
call checksparsity( .true., 15, NCOLMCY, NLENMCY, MX_NCOLMCY, FINMCY, MIDMCY, COLMCY )
480
write( 15, * )'########## FINACV, COLACV, MIDACV ##################'
481
write(15, * )'NCOLACV:', NCOLACV
482
call checksparsity( .true., 15, NCOLACV, CV_PHA_NONODS, MX_NCOLACV, FINACV, MIDACV, COLACV )
484
write( 15, * )'########## FINELE, MIDELE, COLELE ##################'
485
write(15, * )'NCOLELE:',NCOLELE
486
call checksparsity( .true., 15, NCOLELE, TOTELE, MXNELE, FINELE, MIDELE, COLELE )
488
allocate( dummy( CV_NONODS ))
489
write( 15, * )'########## FINDCT, COLCT ##################'
490
write(15, * )'NCOLCT:', NCOLCT
491
call checksparsity( .false., 15, NCOLCT, CV_NONODS, MX_NCT, FINDCT, dummy, COLCT )
494
allocate( dummy( U_NONODS ))
495
write( 15, * )'########## FINDC, COLC ##################'
496
write(15, * )'NCOLC:', NCOLC
497
call checksparsity( .false., 15, NCOLC, U_NONODS, MX_NC, FINDC, dummy, COLC )
500
write( 15, * )'########## FINDGM_PHA, MIDDGM_PHA, COLDGM_PHA ##################'
501
write(15, * )'NCOLDGM_PHA:',NCOLDGM_PHA
502
call checksparsity( .true., 15, NCOLDGM_PHA, U_PHA_NONODS, MX_NCOLDGM_PHA, FINDGM_PHA, MIDDGM_PHA, COLDGM_PHA )
504
write( 15, * )'########## FINDCMC, MIDCMC, COLCMC ##################'
505
write(15, * )'NCOLCMC:',NCOLCMC
506
call checksparsity( .true., 15, NCOLCMC, CV_NONODS, MX_NCOLCMC, FINDCMC, MIDCMC, COLCMC )
508
write( 15, * )'########## FINDM, MIDM, COLM ##################'
509
write(15, * )'NCOLM:',NCOLM
510
call checksparsity( .true., 15, NCOLM, CV_NONODS, MX_NCOLM, FINDM, MIDM, COLM )
514
write(357,*) 'Leaving check_sparsity'
518
end subroutine check_sparsity
520
subroutine mirror_data( unit_debug, problem, nphase, ncomp, totele, ndim, nlev, &
521
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, &
522
cv_snloc, p_snloc, stotel, &
523
ncoef, nuabs_coefs, &
524
u_ele_type, p_ele_type, mat_ele_type, cv_ele_type, &
525
cv_sele_type, u_sele_type, &
526
ntime, nits, ndpset, &
527
dt, patmos, p_ini, t_ini, &
528
t_beta, v_beta, t_theta, v_theta, u_theta, &
529
t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt, &
530
u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt, &
531
domain_length, u_snloc, mat_nloc, cv_nonods, u_nonods, &
532
p_nonods, mat_nonods, ncp_coefs, x_nonods, xu_nonods, &
534
nopt_vel_upwind_coefs, &
535
u_ndgln, xu_ndgln, cv_ndgln, x_ndgln, p_ndgln, &
536
mat_ndgln, u_sndgln, cv_sndgln, x_sndgln, p_sndgln, &
537
wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc, &
538
suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc, suf_p_bc, &
540
suf_u_bc_rob1, suf_u_bc_rob2, &
541
opt_vel_upwind_coefs, &
543
uabs_option, u_abs_stab, u_absorb, &
546
den, satura, comp, p, cv_p, volfra_pore, perm )
549
integer, intent( in ) :: unit_debug,problem, nphase, ncomp, totele, ndim, nlev, &
550
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, &
551
cv_snloc, p_snloc, stotel, &
552
ncoef, nuabs_coefs, &
553
u_ele_type, p_ele_type, mat_ele_type, cv_ele_type, &
554
cv_sele_type, u_sele_type, &
555
ntime, nits, ndpset, &
556
t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt, &
557
u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt, &
558
nopt_vel_upwind_coefs, &
559
u_snloc, mat_nloc, cv_nonods, u_nonods, &
560
p_nonods, mat_nonods, ncp_coefs, x_nonods, xu_nonods, &
562
real, intent( in ) :: dt, patmos, p_ini, t_ini, t_beta, v_beta, &
563
t_theta, v_theta, u_theta, domain_length
564
integer, dimension( stotel * nphase ), intent( in ) :: wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc
565
real, dimension( stotel * cv_snloc * nphase ), intent( in ) :: suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc
566
real, dimension( stotel * p_snloc * nphase ), intent( in ) :: suf_p_bc
567
real, dimension( stotel * u_snloc * nphase ), intent( in ) :: suf_u_bc
568
real, dimension( stotel * u_snloc * nphase ), intent( in ) :: suf_u_bc_rob1, suf_u_bc_rob2
569
integer, dimension( totele * u_nloc ), intent( in ) :: u_ndgln
570
integer, dimension( totele * xu_nloc ), intent( in ) :: xu_ndgln
571
integer, dimension( totele * cv_nloc ), intent( in ) :: cv_ndgln
572
integer, dimension( totele * cv_nloc ), intent( in ) :: x_ndgln
573
integer, dimension( totele * p_nloc ), intent( in ) :: p_ndgln
574
integer, dimension( totele * mat_nloc ), intent( in ) :: mat_ndgln
575
integer, dimension( stotel * u_snloc ), intent( in ) :: u_sndgln
576
integer, dimension( stotel * cv_snloc ), intent( in ) :: cv_sndgln, x_sndgln
577
integer, dimension( stotel * p_snloc ), intent( in ) :: p_sndgln
578
real, dimension( nopt_vel_upwind_coefs ), intent( in ) :: opt_vel_upwind_coefs
579
real, dimension( x_nonods ), intent( in ) :: x
580
real, dimension( xu_nonods ), intent( in ) :: xu
581
real, dimension( u_nonods * nphase ), intent( in ) :: nu, ug
582
integer, dimension( nphase ), intent( in ) :: uabs_option
583
real, dimension( mat_nonods, ndim * nphase, ndim * nphase ), intent( in ) :: u_abs_stab, u_absorb
584
real, dimension( u_nonods * nphase ), intent( in ) :: u_source
585
real, dimension( u_nonods * nphase ), intent( in ) :: u
586
real, dimension( cv_nonods * nphase ), intent( in ) :: den, satura
587
real, dimension( cv_nonods * nphase * ncomp ), intent( in ) :: comp
588
real, dimension( cv_nonods ), intent( in ) :: p, cv_p
589
real, dimension( totele ), intent( in ) :: volfra_pore
590
real, dimension( totele , ndim, ndim ), intent( in ) :: perm
594
write(357,*) 'In mirror_data'
596
write( unit_debug, 201 ) 'problem, nphase, ncomp, totele, ndim, nlev: ', &
597
problem, nphase, ncomp, totele, ndim, nlev
599
write( unit_debug, 201 ) 'u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc: ' , &
600
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc
602
write( unit_debug, 201 ) 'ncoef, nuabs_coefs, u_ele_type, p_ele_type, mat_ele_type: ', &
603
ncoef, nuabs_coefs, u_ele_type, p_ele_type, mat_ele_type
605
write( unit_debug, 201 ) 'cv_ele_type, cv_sele_type, u_sele_type, ntime, nits: ', &
606
cv_ele_type, cv_sele_type, u_sele_type, ntime, nits
608
write( unit_debug, 201 ) 'ndpset, t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt: ', &
609
ndpset, t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt
611
write( unit_debug, 201 ) 'u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt, u_snloc, mat_nloc: ', &
612
u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt, u_snloc, mat_nloc
614
write( unit_debug, 201 ) 'cv_nonods, u_nonods, p_nonods, mat_nonods, ncp_coefs: ', &
615
cv_nonods, u_nonods, p_nonods, mat_nonods, ncp_coefs
617
write( unit_debug, 202 ) 'x_nonods, xu_nonods, nlenmcy: ', &
618
x_nonods, xu_nonods, nlenmcy
620
write( unit_debug, 203 ) 'dt, patmos, p_ini, t_ini: ', dt, patmos, p_ini, t_ini
622
write( unit_debug, 204 ) 't_beta, v_beta, t_theta, v_theta, u_theta: ', &
623
t_beta, v_beta, t_theta, v_theta, u_theta
625
write( unit_debug, * ) 'domain_length: ', domain_length
627
write( unit_debug, * ) '############################################'
629
write( unit_debug, * ) 'wic_vol_bc( stotel * nphase ):', ( wic_vol_bc( i ), i = 1, stotel * nphase )
631
write( unit_debug, * ) 'wic_d_bc( stotel * nphase ):', ( wic_d_bc( i ), i = 1, stotel * nphase )
633
write( unit_debug, * ) 'wic_u_bc( stotel * nphase ):', ( wic_u_bc( i ), i = 1, stotel * nphase )
635
write( unit_debug, * ) 'wic_p_bc( stotel * nphase ):', ( wic_p_bc( i ), i = 1, stotel * nphase )
637
write( unit_debug, * ) 'wic_t_bc( stotel * nphase ):', ( wic_t_bc( i ), i = 1, stotel * nphase )
639
write( unit_debug, * ) '############################################'
641
write( unit_debug, * ) 'suf_vol_bc( stotel * cv_snloc * nphase ):', &
642
( suf_vol_bc( i ), i = 1, stotel * cv_snloc * nphase )
644
write( unit_debug, * ) 'suf_d_bc( stotel * cv_snloc * nphase ):', &
645
( suf_d_bc( i ), i = 1, stotel * cv_snloc * nphase )
647
write( unit_debug, * ) 'suf_cpd_bc( stotel * cv_snloc * nphase ):', &
648
( suf_cpd_bc( i ), i = 1, stotel * cv_snloc * nphase )
650
write( unit_debug, * ) 'suf_t_bc( stotel * cv_snloc * nphase ):', &
651
( suf_t_bc( i ), i = 1, stotel * cv_snloc * nphase )
653
write( unit_debug, * ) 'suf_p_bc( stotel * p_snloc * nphase ):', &
654
( suf_p_bc( i ), i = 1, stotel * p_snloc * nphase )
656
write( unit_debug, * ) 'suf_u_bc( stotel * u_snloc * nphase ):', &
657
( suf_u_bc( i ), i = 1, stotel * u_snloc * nphase )
659
write( unit_debug, * ) 'suf_u_bc_rob1( stotel * u_snloc * nphase ):', &
660
( suf_u_bc_rob1( i ), i = 1, stotel * u_snloc * nphase )
662
write( unit_debug, * ) 'suf_u_bc_rob2( stotel * u_snloc * nphase ):', &
663
( suf_u_bc_rob2( i ), i = 1, stotel * u_snloc * nphase )
665
write( unit_debug, * ) 'suf_t_bc_rob1( stotel * cv_snloc * nphase ):', &
666
( suf_u_bc_rob1( i ), i = 1, stotel * cv_snloc * nphase )
668
write( unit_debug, * ) 'suf_t_bc_rob2( stotel * cv_snloc * nphase ):', &
669
( suf_u_bc_rob2( i ), i = 1, stotel * cv_snloc * nphase )
671
write( unit_debug, * ) '############################################'
673
write( unit_debug, * ) 'opt_vel_upwind_coefs( nopt_vel_upwind_coefs ):', &
674
( opt_vel_upwind_coefs( i ), i = 1, nopt_vel_upwind_coefs )
676
write( unit_debug, * ) 'x( x_nonods ):', ( x( i ), i = 1, x_nonods )
678
write( unit_debug, * ) 'xu( xu_nonods ):', ( xu( i ), i = 1, xu_nonods )
680
write( unit_debug, * ) 'nu( u_nonods * nphase ):', &
681
( nu( i ), i = 1, u_nonods * nphase )
683
write( unit_debug, * ) 'ug( u_nonods * nphase ):', &
684
( ug( i ), i = 1, u_nonods * nphase )
686
write( unit_debug, * ) 'uabs_option( nphase ):', &
687
( uabs_option( i ), i = 1, nphase )
689
write( unit_debug, * ) '############################################'
691
write( unit_debug, * ) 'u_ndgln', size( u_ndgln ), ( u_ndgln( i ), i = 1, totele * u_nloc )
692
write( unit_debug, * ) 'xu_ndgln', size( xu_ndgln ), ( xu_ndgln( i ), i = 1, totele * xu_nloc )
693
write( unit_debug, * ) 'cv_ndgln', size( cv_ndgln ), ( cv_ndgln( i ), i = 1, totele * cv_nloc )
694
write( unit_debug, * ) 'x_ndgln', size( x_ndgln ), ( x_ndgln( i ), i = 1, totele * cv_nloc )
695
write( unit_debug, * ) 'p_ndgln', size( p_ndgln ), ( p_ndgln( i ), i = 1, totele * p_nloc )
696
write( unit_debug, * ) 'mat_ndgln', size( mat_ndgln ), ( mat_ndgln( i ), i = 1, totele * mat_nloc )
697
write( unit_debug, * ) 'u_sndgln', size( u_sndgln), ( u_sndgln( i ), i = 1, stotel * u_snloc )
698
write( unit_debug, * ) 'cv_sndgln', size( cv_sndgln ), ( cv_sndgln( i ), i = 1, stotel * cv_snloc )
699
write( unit_debug, * ) 'x_sndgln', size( x_sndgln ),( x_sndgln( i ), i = 1, stotel * cv_snloc )
700
write( unit_debug, * ) 'p_sndgln', size( p_sndgln ), ( p_sndgln( i ), i = 1, stotel * p_snloc )
702
write( unit_debug, * ) '############################################'
704
write( unit_debug, * ) 'u_abs_stab( mat_nonods, ndim * nphase, ndim * nphase ):'
706
do j = 1, ndim * nphase
707
write( unit_debug, * ) i , j, ( u_abs_stab( i, j, k ), k = 1, ndim * nphase )
711
write( unit_debug, * ) '############################################'
713
write( unit_debug, * ) 'u_absorb( mat_nonods, ndim * nphase, ndim * nphase ):'
715
do j = 1, ndim * nphase
716
write( unit_debug, * ) i , j, ( u_absorb( i, j, k ), k = 1, ndim * nphase )
720
write( unit_debug, * ) '############################################'
722
write( unit_debug, * ) 'u_source( u_nonods * nphase ):', &
723
( u_source( i ), i = 1, u_nonods * nphase )
725
write( unit_debug, * ) '############################################'
728
write( unit_debug, * ) 'u( u_nonods * nphase ):', &
729
( u( i ), i = 1, u_nonods * nphase )
731
write( unit_debug, * ) 'den( cv_nonods * nphase ):', &
732
( den( i ), i = 1, cv_nonods * nphase )
734
write( unit_debug, * ) 'satura( cv_nonods * nphase ):', &
735
( satura( i ), i = 1, cv_nonods * nphase )
737
write( unit_debug, * ) 'comp( cv_nonods * nphase * ncomp ):', &
738
( comp( i ), i = 1, cv_nonods * nphase * ncomp )
740
write( unit_debug, * ) 'p( cv_nonods ):', &
741
( p( i ), i = 1, cv_nonods )
743
write( unit_debug, * ) 'cv_p( cv_nonods ):', &
744
( cv_p( i ), i = 1, cv_nonods )
746
write( unit_debug, * ) 'volfra_pore( totele ):', &
747
( volfra_pore( i ), i = 1, totele )
749
write( unit_debug, * ) 'perm( totele, ndim, ndim ):'
752
write( unit_debug, * ) i , j, ( perm( i, j, k ), k = 1, ndim )
756
!write( unit_debug, * ) '( ):', &
759
write( unit_debug, * ) '############################################'
762
201 format( a, 1x, 5i5 )
763
202 format( a, 1x, 3i5 )
764
203 format( a, 1x, 4g10.4 )
765
204 format( a, 1x, 5g10.4 )
767
write(357,*) 'Leaving mirror_data'
770
end subroutine mirror_data
772
subroutine mirror_array_int( unit_debug, name, ndim, array )
774
integer, intent( in ) :: unit_debug
775
character( len = 50 ), intent( in ) :: name
776
integer, intent( in ) :: ndim
777
integer, dimension( ndim ), intent( in ) :: array
782
k = index( name, ' ' ) - 1
783
write( unit_debug, * ) name( 1 : k ), ( array( idim ), idim = 1, ndim )
786
end subroutine mirror_array_int
789
subroutine mirror_array_real( unit_debug, name, ndim, array )
791
integer, intent( in ) :: unit_debug
792
character( len = 50 ), intent( in ) :: name
793
integer, intent( in ) :: ndim
794
real, dimension( ndim ), intent( in ) :: array
799
k = index( name, ' ' ) - 1
800
write( unit_debug, * ) name( 1 : k ), ( array( idim ), idim = 1, ndim )
803
end subroutine mirror_array_real
806
subroutine mirror_matrix_int( unit_debug, name, ndim1, ndim2, array )
808
integer, intent( in ) :: unit_debug
809
character( len = 50 ), intent( in ) :: name
810
integer, intent( in ) :: ndim1, ndim2
811
integer, dimension( ndim1, ndim2 ), intent( in ) :: array
814
integer :: idim1, idim2, k
816
k = index( name, ' ' ) - 1
818
write( unit_debug, * ) name( 1 : k )
820
write( unit_debug, * ) ( array( idim1, idim2 ), idim2 = 1, ndim2 )
824
end subroutine mirror_matrix_int
826
subroutine mirror_matrix_real( unit_debug, name, ndim1, ndim2, array )
828
integer, intent( in ) :: unit_debug
829
character( len = 50 ), intent( in ) :: name
830
integer, intent( in ) :: ndim1, ndim2
831
real, dimension( ndim1, ndim2 ), intent( in ) :: array
834
integer :: idim1, idim2, k
836
k = index( name, ' ' ) - 1
838
write( unit_debug, * ) name( 1 : k )
840
write( unit_debug, * ) ( array( idim1, idim2 ), idim2 = 1, ndim2 )
844
end subroutine mirror_matrix_real
855
subroutine get_entry( unit, len_name, ior, ifile, fcn_name, value_real, value_bool )
857
integer, intent( in ) :: unit, len_name
858
integer, intent( inout ) :: ior ! ior=13 indicates a non-scalar, ie, array, matrix or tensor
859
character( len = len_name ), intent( inout ) :: ifile, fcn_name
860
real, intent( inout ) :: value_real
861
logical, intent( inout ) :: value_bool
863
character( len = len_name ) :: name
864
integer, parameter :: nexit = 4
865
character( len = nexit ), parameter :: exit_file = 'exit'
869
!write(357,*) 'In get_entry'
879
if( ( name( 1 : 1 ) == '#' ) .or. ( name( 1 : 1 ) == ' ' ) ) then
882
k = index( name, ' ' ) - 1
883
if( name( 1 : nexit ) == exit_file ) then
886
elseif(( name( 1 : 9 ) == 'lump_eqns' ) .or. &
887
( name( 1 : 21 ) == 'volfra_use_theta_flux' ) .or. ( name( 1 : 21 ) == 'volfra_get_theta_flux' ) .or. &
888
( name( 1 : 19 ) == 'comp_use_theta_flux' ) .or. ( name( 1 : 19 ) == 'comp_get_theta_flux' )) then
890
read( unit, * ) name, value_bool
891
k = index( name, ' ' ) - 1
892
ifile( 1 : k ) = name( 1 : k )
896
read( unit, * ) name, value_real
898
if( value_real < -1000. ) then
900
read( unit, * ) name, dummy_real, fcn_name
901
k = index( name, ' ' ) - 1
902
ifile( 1 : k ) = name( 1 : k )
903
fcn_name = trim( fcn_name )
905
elseif( abs( value_real + 1000. ) < 1.e-6 ) then ! so value is -1000
906
! Real or integer array
908
read( unit, * ) name, dummy_real, fcn_name
909
k = index( name, ' ' ) - 1
910
ifile( 1 : k ) = name( 1 : k )
911
fcn_name = trim( fcn_name )
913
elseif( abs( value_real + 999. ) < 1.e-6 ) then ! so value is -999
916
read( unit, * ) name, dummy_real, fcn_name
917
k = index( name, ' ' ) - 1
918
ifile( 1 : k ) = name( 1 : k )
919
fcn_name = trim( fcn_name )
921
elseif( abs( value_real + 998. ) < 1.e-6 ) then ! so value is -998
924
read( unit, * ) name, dummy_real, fcn_name
925
k = index( name, ' ' ) - 1
926
ifile( 1 : k ) = name( 1 : k )
927
fcn_name = trim( fcn_name )
929
elseif( abs( value_real + 997 ) < 1.e-6 ) then ! so value is -997
932
read( unit, * ) name, dummy_real, fcn_name
933
k = index( name, ' ' ) - 1
934
ifile( 1 : k ) = name( 1 : k )
935
fcn_name = trim( fcn_name )
938
k = index( name, ' ' ) - 1
939
ifile( 1 : k ) = name( 1 : k )
947
end subroutine get_entry
949
subroutine Scalar_Int_Val( unit, len_name, fcn_name, scalar_int )
951
integer, intent( in ) :: unit, len_name
952
character( len = len_name ), intent( in ) :: fcn_name
953
integer, intent ( inout ) :: scalar_int
957
k = index( fcn_name, ' ' ) - 1
958
call system( 'rm -f fvalues' )
959
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
960
open( unit + 1, file = 'fvalues', status = 'unknown' )
962
read( unit + 1, * ) scalar_int
967
end subroutine Scalar_Int_Val
969
subroutine Scalar_Real_Val( unit, len_name, fcn_name, scalar_real )
971
integer, intent( in ) :: unit, len_name
972
character( len = len_name ), intent( in ) :: fcn_name
973
real, intent ( inout ) :: scalar_real
977
k = index( fcn_name, ' ' ) - 1
978
call system( 'rm -f fvalues' )
979
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
980
open( unit + 1, file = 'fvalues', status = 'unknown' )
982
read( unit + 1, * ) scalar_real
987
end subroutine Scalar_Real_Val
990
subroutine Array_Int_Set( unit, ior, ndim, len_name, fcn_name, array )
992
integer, intent( in ) :: unit, ior, ndim, len_name
993
character( len = len_name ), intent( in ) :: fcn_name
994
integer, dimension( ndim ), intent( inout ) :: array
1000
call system( 'rm -f filedim' )
1001
open( unit + 1, file = 'filedim', status = 'unknown' )
1002
write( unit + 1, * ) ndim
1007
k = index( fcn_name, ' ' ) - 1
1008
call system( 'rm -f fvalues' )
1009
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
1010
open( unit + 1, file = 'fvalues', status = 'unknown' )
1013
read( unit + 1, * ) array( idim )
1019
end subroutine Array_Int_Set
1022
subroutine Array_Real_Set( unit, ior, ndim, len_name, fcn_name, array )
1024
integer, intent( in ) :: unit, ior, ndim, len_name
1025
character( len = len_name ), intent( in ) :: fcn_name
1026
real, dimension( ndim ), intent( inout ) :: array
1030
print*, 'ndim:', ndim
1031
if (ior == 15 ) then
1033
call system( 'rm -f filedim' )
1034
open( unit + 699, file = 'filedim', status = 'unknown' )
1035
write( unit + 699, * ) ndim
1040
k = index( fcn_name, ' ' ) - 1
1041
call system( 'rm -f fvalues' )
1042
call system( './' // fcn_name( 1 : k ) // ' > fvalues' )
1043
open( unit + 699, file = 'fvalues', status = 'unknown' )
1046
read( unit + 699, * ) array( idim )
1052
end subroutine Array_Real_Set
1055
subroutine Array_Matrix2_Set( unit, ior, ndim1, ndim2, len_name, fcn_name, matrix )
1057
integer, intent( in ) :: unit, ior, ndim1, ndim2, len_name
1058
character( len = len_name ), intent( in ) :: fcn_name
1059
real, dimension( ndim1, ndim2 ), intent( inout ) :: matrix
1061
integer :: idim1, idim2, k
1063
if (ior == 16 ) then
1065
open( unit + 563, file = 'filedim', status = 'unknown' )
1066
write( unit + 563, * ) ndim1, ndim2
1071
k = index( fcn_name, ' ' ) - 1
1072
call system( 'rm -f fvalues' )
1073
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
1074
open( unit + 563, file = 'fvalues', status = 'unknown' )
1078
read( unit + 563, * ) matrix( idim1, idim2 )
1085
end subroutine Array_Matrix2_Set
1089
subroutine Array_Matrix3_Set( unit, ior, ndim1, ndim2, ndim3, len_name, fcn_name, matrix )
1091
integer, intent( in ) :: unit, ior, ndim1, ndim2, ndim3, len_name
1092
character( len = len_name ), intent( in ) :: fcn_name
1093
real, dimension( ndim1, ndim2, ndim3 ), intent( inout ) :: matrix
1095
integer :: idim1, idim2, idim3, k
1097
if (ior == 17 ) then
1099
open( unit + 1, file = 'filedim', status = 'unknown' )
1100
write( unit + 1, * ) ndim1, ndim2, ndim3
1104
k = index( fcn_name, ' ' ) - 1
1105
call system( 'rm -f fvalues' )
1106
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
1107
open( unit + 1, file = 'fvalues', status = 'unknown' )
1112
read( unit + 1, * ) matrix( idim1, idim2, idim3 )
1120
end subroutine Array_Matrix3_Set
1123
subroutine Array_Matrix4_Set( unit, ior, ndim1, ndim2, ndim3, ndim4, len_name, fcn_name, matrix )
1125
integer, intent( in ) :: unit, ior, ndim1, ndim2, ndim3, ndim4, len_name
1126
character( len = len_name ), intent( in ) :: fcn_name
1127
real, dimension( ndim1, ndim2, ndim3, ndim4 ), intent( inout ) :: matrix
1129
integer :: idim1, idim2, idim3, idim4, k
1131
if (ior == 18 ) then
1133
open( unit + 1, file = 'filedim', status = 'unknown' )
1134
write( unit + 1, * ) ndim1, ndim2, ndim3, ndim4
1139
k = index( fcn_name, ' ' ) - 1
1140
call system( 'rm -f fvalues' )
1141
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
1142
open( unit + 1, file = 'fvalues', status = 'unknown' )
1148
read( unit + 1, * ) matrix( idim1, idim2, idim3, idim4 )
1157
end subroutine Array_Matrix4_Set
1160
subroutine read_scalar( unit, option_debug, problem, nphase, ncomp, totele, ndim, nlev, &
1161
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, &
1162
cv_snloc, u_snloc, p_snloc, x_snloc, stotel, &
1163
ncoef, nuabs_coefs, &
1164
u_ele_type, p_ele_type, mat_ele_type, cv_ele_type, &
1165
cv_sele_type, u_sele_type, ntime, nits, nits_internal, noit_dim, &
1166
nits_flux_lim_volfra, nits_flux_lim_comp, &
1168
dt, patmos, p_ini, t_ini, &
1169
t_beta, v_beta, t_theta, v_theta, u_theta, &
1170
t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt, &
1171
u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt, &
1173
volfra_use_theta_flux, volfra_get_theta_flux, comp_use_theta_flux, comp_get_theta_flux, &
1174
capil_pres_opt, ncapil_pres_coef, comp_diffusion_opt, ncomp_diff_coef, &
1178
integer, intent( in ) :: unit
1179
integer, intent( inout ) :: option_debug, problem, nphase, ncomp, totele, ndim, nlev, &
1180
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, &
1181
cv_snloc, u_snloc, p_snloc, x_snloc, stotel, &
1182
ncoef, nuabs_coefs, &
1183
u_ele_type, p_ele_type, mat_ele_type, cv_ele_type, &
1184
cv_sele_type, u_sele_type, ntime, nits, nits_internal, noit_dim, &
1185
nits_flux_lim_volfra, nits_flux_lim_comp, &
1187
real, intent( inout ) :: dt, patmos, p_ini, t_ini, &
1188
t_beta, v_beta, t_theta, v_theta, u_theta
1189
integer, intent( inout ) :: t_disopt, u_disopt, v_disopt, t_dg_vel_int_opt, &
1190
u_dg_vel_int_opt, v_dg_vel_int_opt, w_dg_vel_int_opt
1191
logical, intent( inout ) :: lump_eqns, &
1192
volfra_use_theta_flux, volfra_get_theta_flux, comp_use_theta_flux, comp_get_theta_flux
1193
integer, intent( inout ) :: capil_pres_opt, ncapil_pres_coef, &
1194
comp_diffusion_opt, ncomp_diff_coef
1195
real, intent( inout ) :: domain_length
1197
integer, parameter :: len_name = 50
1198
character( len = len_name ) :: ifile, fcn_name
1201
logical :: value_bool
1206
do while ( .not. ( ior == -1000 ))
1207
call get_entry( unit, len_name, ior, ifile, fcn_name, value_real, value_bool )
1209
Conditional_IOR: if ( ior > 0 ) then
1210
k = index( ifile, ' ') - 1
1212
Select Case ( ifile( 1 : k ))
1214
Case( 'option_debug' );
1215
option_debug = real( value_real )
1218
problem = real( value_real )
1221
nphase = int( value_real )
1224
ncomp = int( value_real )
1227
totele = int( value_real )
1230
ndim = int( value_real )
1233
nlev = int( value_real )
1236
u_nloc = int( value_real )
1239
xu_nloc = int( value_real )
1242
cv_nloc = int( value_real )
1245
x_nloc = int( value_real )
1248
p_nloc = int( value_real )
1251
cv_snloc = int( value_real )
1254
u_snloc = int( value_real )
1257
p_snloc = int( value_real )
1260
x_snloc = int( value_real )
1263
stotel = int( value_real )
1266
ncoef = int( value_real )
1268
Case( 'nuabs_coefs' );
1269
nuabs_coefs = int( value_real )
1271
Case( 'u_ele_type' );
1272
u_ele_type = int( value_real )
1274
Case( 'p_ele_type' );
1275
p_ele_type = int( value_real )
1277
Case( 'mat_ele_type' );
1278
mat_ele_type = int( value_real )
1280
Case( 'cv_ele_type' );
1281
cv_ele_type = int( value_real )
1283
Case( 'cv_sele_type' );
1284
cv_sele_type = int( value_real )
1286
Case( 'u_sele_type' );
1287
u_sele_type = int( value_real )
1290
ntime = int( value_real )
1293
nits = int( value_real )
1295
Case( 'nits_internal' );
1296
nits_internal = int( value_real )
1299
noit_dim = int( value_real )
1301
Case( 'nits_flux_lim_volfra' );
1302
nits_flux_lim_volfra = int( value_real )
1304
Case( 'nits_flux_lim_comp' );
1305
nits_flux_lim_comp = int( value_real )
1308
ndpset = int( value_real )
1311
v_disopt = int( value_real )
1314
t_disopt = int( value_real )
1317
u_disopt = int( value_real )
1319
Case( 't_dg_vel_int_opt' );
1320
t_dg_vel_int_opt = int( value_real )
1322
Case( 'u_dg_vel_int_opt' );
1323
u_dg_vel_int_opt = int( value_real )
1325
Case( 'v_dg_vel_int_opt' );
1326
v_dg_vel_int_opt = int( value_real )
1328
Case( 'w_dg_vel_int_opt' );
1329
w_dg_vel_int_opt = int( value_real )
1331
Case( 'lump_eqns' );
1332
lump_eqns = value_bool
1334
Case( 'volfra_use_theta_flux' );
1335
volfra_use_theta_flux = value_bool
1337
Case( 'volfra_get_theta_flux' );
1338
volfra_get_theta_flux = value_bool
1340
Case( 'comp_use_theta_flux' );
1341
comp_use_theta_flux = value_bool
1343
Case( 'comp_get_theta_flux' );
1344
comp_get_theta_flux = value_bool
1346
Case( 'capil_pres_opt' );
1347
capil_pres_opt = int( value_real )
1349
Case( 'ncapil_pres_coef' );
1350
ncapil_pres_coef = int( value_real )
1352
Case( 'comp_diffusion_opt' );
1353
comp_diffusion_opt = int( value_real )
1355
Case( 'ncomp_diff_coef' );
1356
ncomp_diff_coef = int( value_real )
1358
Case( 'domain_length' );
1359
domain_length = int( value_real )
1380
t_theta = value_real
1383
v_theta = value_real
1386
u_theta = value_real
1389
write( 357, * ) 'Option not found in read_scalar subrt.', ifile( 1 : k )
1394
end if Conditional_IOR
1400
end subroutine read_scalar
1403
subroutine read_all( unit, nphase, ncomp, totele, ndim, &
1404
cv_snloc, u_snloc, p_snloc, stotel, &
1405
ncoef, nuabs_coefs, ncp_coefs, &
1406
cv_nonods, u_nonods, &
1407
mat_nonods, x_nonods, xu_nonods, &
1408
nopt_vel_upwind_coefs, &
1409
wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc, wic_comp_bc, &
1410
suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc, suf_p_bc, &
1411
suf_u_bc, suf_v_bc, suf_w_bc, suf_one_bc, suf_comp_bc, &
1412
suf_u_bc_rob1, suf_u_bc_rob2, suf_v_bc_rob1, suf_v_bc_rob2, &
1413
suf_w_bc_rob1, suf_w_bc_rob2, suf_t_bc_rob1, suf_t_bc_rob2, &
1414
suf_vol_bc_rob1, suf_vol_bc_rob2, &
1415
suf_comp_bc_rob1, suf_comp_bc_rob2, &
1416
opt_vel_upwind_coefs, &
1417
volfra_error, volfra_relax, volfra_relax_diag, volfra_relax_row, volfra_relax_number_iterations, &
1418
scalar_error, scalar_relax, scalar_relax_diag, scalar_relax_row, scalar_relax_number_iterations, &
1419
global_error, global_relax, global_relax_diag, global_relax_row, global_relax_number_iterations, &
1420
velocity_error, velocity_relax, velocity_relax_diag, velocity_relax_row, velocity_relax_number_iterations, &
1421
pressure_error, pressure_relax, pressure_relax_diag, pressure_relax_row, pressure_relax_number_iterations, &
1422
mass_matrix_error, mass_matrix_relax, mass_matrix_relax_diag, mass_matrix_relax_row, mass_matrix_relax_number_iterations, &
1423
in_ele_upwind, dg_ele_upwind, &
1424
x, y, z, xu, yu, zu, nu, nv, nw, ug, vg, wg, &
1425
uabs_option, uabs_coefs, u_abs_stab, &
1426
u_absorb, t_absorb, v_absorb, comp_absorb, &
1427
u_source, t_source, v_source, comp_source, udiffusion, tdiffusion, &
1428
ncomp_diff_coef, comp_diffusion, comp_diff_coef, &
1429
ncapil_pres_coef, capil_pres_coef, &
1431
den, satura, comp, volfra, t, cv_one, p, cv_p, volfra_pore, perm, &
1432
K_Comp, alpha_beta, &
1433
eos_option, cp_option, eos_coefs, cp_coefs )
1436
integer, intent( in ) :: unit, nphase, ncomp, totele, ndim, &
1437
cv_snloc, u_snloc, p_snloc, stotel, &
1438
ncoef, nuabs_coefs, ncp_coefs, &
1439
cv_nonods, u_nonods, &
1440
mat_nonods, x_nonods, xu_nonods, &
1441
nopt_vel_upwind_coefs
1442
integer, dimension( stotel * nphase ), intent( inout ) :: wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc, wic_comp_bc
1443
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc
1444
real, dimension( stotel * p_snloc * nphase ), intent( inout ) :: suf_p_bc
1445
real, dimension( stotel * u_snloc * nphase ), intent( inout ) :: suf_u_bc, suf_v_bc, suf_w_bc
1446
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_one_bc
1447
real, dimension( stotel * cv_snloc * nphase * ncomp ), intent( inout ) :: suf_comp_bc
1448
real, dimension( stotel * u_snloc * nphase ), intent( inout ) :: suf_u_bc_rob1, suf_u_bc_rob2, suf_v_bc_rob1, suf_v_bc_rob2, &
1449
suf_w_bc_rob1, suf_w_bc_rob2
1450
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_t_bc_rob1, suf_t_bc_rob2
1451
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_vol_bc_rob1, suf_vol_bc_rob2
1452
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_comp_bc_rob1, suf_comp_bc_rob2
1453
real, dimension( nopt_vel_upwind_coefs ), intent( inout ) :: opt_vel_upwind_coefs
1454
real, intent( inout ) :: volfra_error, volfra_relax, volfra_relax_diag, volfra_relax_row, &
1455
scalar_error, scalar_relax, scalar_relax_diag, scalar_relax_row, &
1456
global_error, global_relax, global_relax_diag, global_relax_row, &
1457
velocity_error, velocity_relax, velocity_relax_diag, velocity_relax_row, &
1458
pressure_error, pressure_relax, pressure_relax_diag, pressure_relax_row, &
1459
mass_matrix_error, mass_matrix_relax, mass_matrix_relax_diag, mass_matrix_relax_row
1460
integer, intent( inout ) :: volfra_relax_number_iterations, scalar_relax_number_iterations, &
1461
global_relax_number_iterations, velocity_relax_number_iterations, &
1462
pressure_relax_number_iterations, mass_matrix_relax_number_iterations, &
1463
in_ele_upwind, dg_ele_upwind
1464
real, dimension( x_nonods ), intent( inout ) :: x, y, z
1465
real, dimension( xu_nonods ), intent( inout ) :: xu, yu, zu
1466
real, dimension( u_nonods * nphase ), intent( inout ) :: nu, nv, nw, ug, vg, wg
1467
integer, dimension( nphase ), intent( inout ) :: uabs_option
1468
real, dimension( nphase, nuabs_coefs ), intent( inout ) :: uabs_coefs
1469
real, dimension( mat_nonods, ndim * nphase, ndim * nphase ), intent( inout ) :: u_abs_stab, u_absorb
1470
real, dimension( cv_nonods, nphase, nphase ), intent( inout ) :: t_absorb, v_absorb
1471
real, dimension( cv_nonods * nphase, nphase, nphase ), intent( inout ) :: comp_absorb
1472
real, dimension( u_nonods * nphase ), intent( inout ) :: u_source
1473
real, dimension( cv_nonods * nphase ), intent( inout ) :: t_source, v_source, comp_source
1474
real, dimension( mat_nonods, ndim, ndim, nphase ), intent( inout ) :: udiffusion, tdiffusion
1476
integer, intent( in ) :: ncomp_diff_coef
1477
real, dimension( mat_nonods, ndim, ndim, nphase ), intent( inout ) :: comp_diffusion
1478
real, dimension( ncomp, ncomp_diff_coef, nphase ), intent( inout ) :: comp_diff_coef
1479
integer, intent( in ) :: ncapil_pres_coef
1480
real, dimension( ncapil_pres_coef, nphase, nphase ), intent( inout ) :: &
1483
real, dimension( u_nonods * nphase ), intent( inout ) :: u, v, w
1484
real, dimension( cv_nonods * nphase ), intent( inout ) :: den, satura
1485
real, dimension( cv_nonods * nphase * ncomp ), intent( inout ) :: comp
1486
real, dimension( cv_nonods * nphase ), intent( inout ) :: volfra, t, cv_one
1487
real, dimension( cv_nonods ), intent( inout ) :: p, cv_p
1488
real, dimension( totele ), intent( inout ) :: volfra_pore
1489
real, dimension( totele , ndim, ndim ), intent( inout ) :: perm
1490
real, dimension( ncomp, nphase, nphase ), intent( inout ) :: K_Comp
1491
real, intent( inout ) :: alpha_beta
1492
integer, dimension( nphase ), intent( inout ) :: eos_option, cp_option
1493
real, dimension( nphase, ncoef ), intent( inout ) :: eos_coefs
1494
real, dimension( nphase, ncp_coefs ), intent( inout ) :: cp_coefs
1498
integer, parameter :: len_name = 50
1499
character( len = len_name ) :: ifile, fcn_name
1502
logical :: value_bool
1504
write(357,*) 'In read_all'
1507
do while ( .not. ( ior == -1000 ))
1510
call get_entry( unit, len_name, ior, ifile, fcn_name, value_real, value_bool )
1512
Conditional_IOR: if ( ior > 0 ) then
1513
k = index( ifile, ' ') - 1
1514
Select Case ( ifile( 1 : k ))
1517
Case( 'alpha_beta' );
1518
if( ior == 13 ) then
1519
call Scalar_Real_Val( unit, len_name, fcn_name, alpha_beta )
1521
alpha_beta = value_real
1524
Case( 'volfra_error' );
1525
if( ior == 13 ) then
1526
call Scalar_Real_Val( unit, len_name, fcn_name, volfra_error )
1528
volfra_error = value_real
1531
Case( 'scalar_error' );
1532
if( ior == 13 ) then
1533
call Scalar_Real_Val( unit, len_name, fcn_name, scalar_error )
1535
scalar_error = value_real
1538
Case( 'global_error' );
1539
if( ior == 13 ) then
1540
call Scalar_Real_Val( unit, len_name, fcn_name, global_error )
1542
global_error = value_real
1545
Case( 'velocity_error' );
1546
if( ior == 13 ) then
1547
call Scalar_Real_Val( unit, len_name, fcn_name, velocity_error )
1549
velocity_error = value_real
1552
Case( 'pressure_error' );
1553
if( ior == 13 ) then
1554
call Scalar_Real_Val( unit, len_name, fcn_name, pressure_error )
1556
pressure_error = value_real
1559
Case( 'mass_matrix_error' );
1560
if( ior == 13 ) then
1561
call Scalar_Real_Val( unit, len_name, fcn_name, mass_matrix_error )
1563
mass_matrix_error = value_real
1566
Case( 'volfra_relax' );
1567
if( ior == 13 ) then
1568
call Scalar_Real_Val( unit, len_name, fcn_name, volfra_relax )
1570
volfra_relax = value_real
1573
Case( 'scalar_relax' );
1574
if( ior == 13 ) then
1575
call Scalar_Real_Val( unit, len_name, fcn_name, scalar_relax )
1577
scalar_relax = value_real
1580
Case( 'global_relax' );
1581
if( ior == 13 ) then
1582
call Scalar_Real_Val( unit, len_name, fcn_name, global_relax )
1584
global_relax = value_real
1587
Case( 'velocity_relax' );
1588
if( ior == 13 ) then
1589
call Scalar_Real_Val( unit, len_name, fcn_name, velocity_relax )
1591
velocity_relax = value_real
1594
Case( 'pressure_relax' );
1595
if( ior == 13 ) then
1596
call Scalar_Real_Val( unit, len_name, fcn_name, pressure_relax )
1598
pressure_relax = value_real
1601
Case( 'mass_matrix_relax' );
1602
if( ior == 13 ) then
1603
call Scalar_Real_Val( unit, len_name, fcn_name, mass_matrix_relax )
1605
mass_matrix_relax = value_real
1608
Case( 'volfra_relax_diag' );
1609
if( ior == 13 ) then
1610
call Scalar_Real_Val( unit, len_name, fcn_name, volfra_relax_diag )
1612
volfra_relax_diag = value_real
1615
Case( 'scalar_relax_diag' );
1616
if( ior == 13 ) then
1617
call Scalar_Real_Val( unit, len_name, fcn_name, scalar_relax_diag )
1619
scalar_relax_diag = value_real
1622
Case( 'global_relax_diag' );
1623
if( ior == 13 ) then
1624
call Scalar_Real_Val( unit, len_name, fcn_name, global_relax_diag )
1626
global_relax_diag = value_real
1629
Case( 'velocity_relax_diag' );
1630
if( ior == 13 ) then
1631
call Scalar_Real_Val( unit, len_name, fcn_name, velocity_relax_diag )
1633
velocity_relax_diag = value_real
1636
Case( 'pressure_relax_diag' );
1637
if( ior == 13 ) then
1638
call Scalar_Real_Val( unit, len_name, fcn_name, pressure_relax_diag )
1640
pressure_relax_diag = value_real
1643
Case( 'mass_matrix_relax_diag' );
1644
if( ior == 13 ) then
1645
call Scalar_Real_Val( unit, len_name, fcn_name, mass_matrix_relax_diag )
1647
mass_matrix_relax_diag = value_real
1650
Case( 'volfra_relax_row' );
1651
if( ior == 13 ) then
1652
call Scalar_Real_Val( unit, len_name, fcn_name, volfra_relax_row )
1654
volfra_relax_row = value_real
1657
Case( 'scalar_relax_row' );
1658
if( ior == 13 ) then
1659
call Scalar_Real_Val( unit, len_name, fcn_name, scalar_relax_row )
1661
scalar_relax_row = value_real
1664
Case( 'global_relax_row' );
1665
if( ior == 13 ) then
1666
call Scalar_Real_Val( unit, len_name, fcn_name, global_relax_row )
1668
global_relax_row = value_real
1671
Case( 'velocity_relax_row' );
1672
if( ior == 13 ) then
1673
call Scalar_Real_Val( unit, len_name, fcn_name, velocity_relax_row )
1675
velocity_relax_row = value_real
1678
Case( 'pressure_relax_row' );
1679
if( ior == 13 ) then
1680
call Scalar_Real_Val( unit, len_name, fcn_name, pressure_relax_row )
1682
pressure_relax_row = value_real
1685
Case( 'mass_matrix_relax_row' );
1686
if( ior == 13 ) then
1687
call Scalar_Real_Val( unit, len_name, fcn_name, mass_matrix_relax_row )
1689
mass_matrix_relax_row = value_real
1692
Case( 'volfra_relax_number_iterations' );
1693
if( ior == 13 ) then
1694
call Scalar_Int_Val( unit, len_name, fcn_name, volfra_relax_number_iterations )
1696
volfra_relax_number_iterations = int( value_real )
1699
Case( 'scalar_relax_number_iterations' );
1700
if( ior == 13 ) then
1701
call Scalar_Int_Val( unit, len_name, fcn_name, scalar_relax_number_iterations )
1703
scalar_relax_number_iterations = int( value_real )
1706
Case( 'global_relax_number_iterations' );
1707
if( ior == 13 ) then
1708
call Scalar_Int_Val( unit, len_name, fcn_name, global_relax_number_iterations )
1710
global_relax_number_iterations = int( value_real )
1713
Case( 'velocity_relax_number_iterations' );
1714
if( ior == 13 ) then
1715
call Scalar_Int_Val( unit, len_name, fcn_name, velocity_relax_number_iterations )
1717
velocity_relax_number_iterations = int( value_real )
1720
Case( 'pressure_relax_number_iterations' );
1721
if( ior == 13 ) then
1722
call Scalar_Int_Val( unit, len_name, fcn_name, pressure_relax_number_iterations )
1724
pressure_relax_number_iterations = int( value_real )
1727
Case( 'mass_matrix_relax_number_iterations' );
1728
if( ior == 13 ) then
1729
call Scalar_Int_Val( unit, len_name, fcn_name, mass_matrix_relax_number_iterations )
1731
mass_matrix_relax_number_iterations = int( value_real )
1734
Case( 'in_ele_upwind' );
1735
if( ior == 13 ) then
1736
call Scalar_Int_Val( unit, len_name, fcn_name, in_ele_upwind )
1738
in_ele_upwind = int( value_real )
1741
Case( 'dg_ele_upwind' );
1742
if( ior == 13 ) then
1743
call Scalar_Int_Val( unit, len_name, fcn_name, dg_ele_upwind )
1745
dg_ele_upwind = int( value_real )
1750
Case( 'wic_vol_bc' );
1752
if(( ior == 13 ) .or. ( ior == 15 )) then
1753
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_vol_bc )
1755
wic_vol_bc = int( value_real )
1759
if(( ior == 13 ) .or. ( ior == 15 )) then
1760
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_d_bc )
1762
wic_d_bc = int( value_real )
1766
if(( ior == 13 ) .or. ( ior == 15 )) then
1767
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_u_bc )
1769
wic_u_bc = int( value_real )
1773
if(( ior == 13 ) .or. ( ior == 15 )) then
1774
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_p_bc )
1776
wic_p_bc = int( value_real )
1780
if(( ior == 13 ) .or. ( ior == 15 )) then
1781
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_t_bc )
1783
wic_t_bc = int( value_real )
1786
Case( 'wic_comp_bc' );
1787
if(( ior == 13 ) .or. ( ior == 15 )) then
1788
call Array_Int_Set( unit, ior, stotel * nphase, len_name, fcn_name, wic_comp_bc )
1790
wic_comp_bc = int( value_real )
1793
Case( 'uabs_option' );
1794
if(( ior == 13 ) .or. ( ior == 15 )) then
1795
call Array_Int_Set( unit, ior, nphase, len_name, fcn_name, uabs_option )
1797
uabs_option = int( value_real )
1800
Case( 'eos_option' );
1801
if(( ior == 13 ) .or. ( ior == 15 )) then
1802
call Array_Int_Set( unit, ior, nphase, len_name, fcn_name, eos_option )
1804
eos_option = int( value_real )
1807
Case( 'cp_option' );
1808
if(( ior == 13 ) .or. ( ior == 15 )) then
1809
call Array_Int_Set( unit, ior, nphase, len_name, fcn_name, cp_option )
1811
cp_option = int( value_real )
1815
Case( 'suf_vol_bc' );
1816
if(( ior == 13 ) .or. ( ior == 15 )) then
1817
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_vol_bc )
1819
suf_vol_bc = value_real
1823
if(( ior == 13 ) .or. ( ior == 15 )) then
1824
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_d_bc )
1826
suf_d_bc = value_real
1829
Case( 'suf_cpd_bc' );
1830
if(( ior == 13 ) .or. ( ior == 15 )) then
1831
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_cpd_bc )
1833
suf_cpd_bc = value_real
1837
if(( ior == 13 ) .or. ( ior == 15 )) then
1838
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_t_bc )
1840
suf_t_bc = value_real
1844
if(( ior == 13 ) .or. ( ior == 15 )) then
1845
call Array_Real_Set( unit, ior, stotel * p_snloc * nphase, len_name, fcn_name, suf_p_bc )
1847
suf_p_bc = value_real
1851
if(( ior == 13 ) .or. ( ior == 15 )) then
1852
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_u_bc )
1854
suf_u_bc = value_real
1858
if(( ior == 13 ) .or. ( ior == 15 )) then
1859
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_v_bc )
1861
suf_v_bc = value_real
1865
if(( ior == 13 ) .or. ( ior == 15 )) then
1866
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_w_bc )
1868
suf_w_bc = value_real
1871
Case( 'suf_one_bc' );
1872
if(( ior == 13 ) .or. ( ior == 15 )) then
1873
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_one_bc )
1875
suf_one_bc = value_real
1878
Case( 'suf_comp_bc' );
1879
if(( ior == 13 ) .or. ( ior == 15 )) then
1880
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase * ncomp, len_name, fcn_name, suf_comp_bc )
1882
suf_comp_bc = value_real
1885
Case( 'suf_u_bc_rob1' );
1886
if(( ior == 13 ) .or. ( ior == 15 )) then
1887
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_u_bc_rob1 )
1889
suf_u_bc_rob1 = value_real
1892
Case( 'suf_u_bc_rob2' );
1893
if(( ior == 13 ) .or. ( ior == 15 )) then
1894
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_u_bc_rob2 )
1896
suf_u_bc_rob2 = value_real
1899
Case( 'suf_v_bc_rob1' );
1900
if(( ior == 13 ) .or. ( ior == 15 )) then
1901
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_v_bc_rob1 )
1903
suf_v_bc_rob1 = value_real
1906
Case( 'suf_v_bc_rob2' );
1907
if(( ior == 13 ) .or. ( ior == 15 )) then
1908
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_v_bc_rob2 )
1910
suf_v_bc_rob2 = value_real
1913
Case( 'suf_w_bc_rob1' );
1914
if(( ior == 13 ) .or. ( ior == 15 )) then
1915
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_w_bc_rob1 )
1917
suf_w_bc_rob1 = value_real
1920
Case( 'suf_w_bc_rob2' );
1921
if(( ior == 13 ) .or. ( ior == 15 )) then
1922
call Array_Real_Set( unit, ior, stotel * u_snloc * nphase, len_name, fcn_name, suf_w_bc_rob2 )
1924
suf_w_bc_rob2 = value_real
1927
Case( 'suf_t_bc_rob1' );
1928
if(( ior == 13 ) .or. ( ior == 15 )) then
1929
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_t_bc_rob1 )
1931
suf_t_bc_rob1 = value_real
1934
Case( 'suf_t_bc_rob2' );
1935
if(( ior == 13 ) .or. ( ior == 15 )) then
1936
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_t_bc_rob2 )
1938
suf_t_bc_rob2 = value_real
1941
Case( 'suf_vol_bc_rob1' );
1942
if(( ior == 13 ) .or. ( ior == 15 )) then
1943
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_vol_bc_rob1 )
1945
suf_vol_bc_rob1 = value_real
1948
Case( 'suf_vol_bc_rob2' );
1949
if(( ior == 13 ) .or. ( ior == 15 )) then
1950
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_vol_bc_rob2 )
1952
suf_vol_bc_rob2 = value_real
1955
Case( 'suf_comp_bc_rob1' );
1956
if(( ior == 13 ) .or. ( ior == 15 )) then
1957
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_comp_bc_rob1 )
1959
suf_comp_bc_rob1 = value_real
1962
Case( 'suf_comp_bc_rob2' );
1963
if(( ior == 13 ) .or. ( ior == 15 )) then
1964
call Array_Real_Set( unit, ior, stotel * cv_snloc * nphase, len_name, fcn_name, suf_comp_bc_rob2 )
1966
suf_comp_bc_rob2 = value_real
1969
Case( 'opt_vel_upwind_coefs' );
1970
if(( ior == 13 ) .or. ( ior == 15 )) then
1971
call Array_Real_Set( unit, ior, nopt_vel_upwind_coefs, len_name, fcn_name, opt_vel_upwind_coefs )
1973
opt_vel_upwind_coefs = value_real
1977
if(( ior == 13 ) .or. ( ior == 15 )) then
1978
call Array_Real_Set( unit, ior, x_nonods, len_name, fcn_name, x )
1984
if(( ior == 13 ) .or. ( ior == 15 )) then
1985
call Array_Real_Set( unit, ior, x_nonods, len_name, fcn_name, y )
1991
if(( ior == 13 ) .or. ( ior == 15 )) then
1992
call Array_Real_Set( unit, ior, x_nonods, len_name, fcn_name, z )
1998
if(( ior == 13 ) .or. ( ior == 15 )) then
1999
call Array_Real_Set( unit, ior, xu_nonods, len_name, fcn_name, xu )
2005
if(( ior == 13 ) .or. ( ior == 15 )) then
2006
call Array_Real_Set( unit, ior, xu_nonods, len_name, fcn_name, yu )
2012
if(( ior == 13 ) .or. ( ior == 15 )) then
2013
call Array_Real_Set( unit, ior, xu_nonods, len_name, fcn_name, zu )
2019
if(( ior == 13 ) .or. ( ior == 15 )) then
2020
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, nu )
2026
if(( ior == 13 ) .or. ( ior == 15 )) then
2027
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, nv )
2033
if(( ior == 13 ) .or. ( ior == 15 )) then
2034
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, nw )
2040
if(( ior == 13 ) .or. ( ior == 15 )) then
2041
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, ug )
2047
if(( ior == 13 ) .or. ( ior == 15 )) then
2048
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, vg )
2054
if(( ior == 13 ) .or. ( ior == 15 )) then
2055
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, wg )
2061
if(( ior == 13 ) .or. ( ior == 15 )) then
2062
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, u_source )
2064
u_source = value_real
2068
if(( ior == 13 ) .or. ( ior == 15 )) then
2069
call Array_Real_Set( unit, ior, cv_nonods * nphase, len_name, fcn_name, t_source )
2071
t_source = value_real
2075
if(( ior == 13 ) .or. ( ior == 15 )) then
2076
call Array_Real_Set( unit, ior, cv_nonods * nphase, len_name, fcn_name, v_source )
2078
v_source = value_real
2081
Case( 'comp_source' );
2082
if(( ior == 13 ) .or. ( ior == 15 )) then
2083
call Array_Real_Set( unit, ior, cv_nonods * nphase, len_name, fcn_name, comp_source )
2085
comp_source = value_real
2089
if(( ior == 13 ) .or. ( ior == 15 )) then
2090
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, u )
2096
if(( ior == 13 ) .or. ( ior == 15 )) then
2097
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, v )
2103
if(( ior == 13 ) .or. ( ior == 15 )) then
2104
call Array_Real_Set( unit, ior, u_nonods * nphase, len_name, fcn_name, w )
2110
if(( ior == 13 ) .or. ( ior == 15 )) then
2111
call Array_Real_Set( unit, ior, cv_nonods * nphase, len_name, fcn_name, den )
2117
if(( ior == 13 ) .or. ( ior == 15 )) then
2118
call Array_Real_Set( unit, ior, cv_nonods * nphase, len_name, fcn_name, satura )
2124
if(( ior == 13 ) .or. ( ior == 15 )) then
2125
call Array_Real_Set( unit, ior, cv_nonods * nphase * ncomp, len_name, fcn_name, comp )
2131
if(( ior == 13 ) .or. ( ior == 15 )) then
2132
call Array_Real_Set( unit, ior, cv_nonods * nphase , len_name, fcn_name, volfra )
2138
if(( ior == 13 ) .or. ( ior == 15 )) then
2139
call Array_Real_Set( unit, ior, cv_nonods * nphase , len_name, fcn_name, t )
2145
if(( ior == 13 ) .or. ( ior == 15 )) then
2146
call Array_Real_Set( unit, ior, cv_nonods * nphase , len_name, fcn_name, cv_one )
2152
if(( ior == 13 ) .or. ( ior == 15 )) then
2153
call Array_Real_Set( unit, ior, cv_nonods, len_name, fcn_name, p )
2159
if(( ior == 13 ) .or. ( ior == 15 )) then
2160
call Array_Real_Set( unit, ior, cv_nonods, len_name, fcn_name, cv_p )
2165
Case( 'volfra_pore' );
2166
if(( ior == 13 ) .or. ( ior == 15 )) then
2167
call Array_Real_Set( unit, ior, totele, len_name, fcn_name, volfra_pore )
2169
volfra_pore = value_real
2173
if(( ior == 13 ) .or. ( ior == 17 )) then
2174
call Array_Matrix3_Set( unit, ior, ncomp, nphase, nphase, len_name, fcn_name, K_Comp )
2179
! Matrices and Tensors ( x , y )
2180
Case( 'uabs_coefs' );
2181
if(( ior == 13 ) .or. ( ior == 16 )) then
2182
call Array_Matrix2_Set( unit, ior, nphase, nuabs_coefs, len_name, fcn_name, uabs_coefs )
2184
uabs_coefs = value_real
2187
Case( 'eos_coefs' );
2188
if(( ior == 13 ) .or. ( ior == 16 )) then
2189
call Array_Matrix2_Set( unit, ior, nphase, ncoef, len_name, fcn_name, eos_coefs )
2191
eos_coefs = value_real
2195
if(( ior == 13 ) .or. ( ior == 16 )) then
2196
call Array_Matrix2_Set( unit, ior, nphase, ncp_coefs, len_name, fcn_name, cp_coefs )
2198
cp_coefs = value_real
2201
! Matrices and Tensors ( x , y , z)
2203
Case( 'comp_diff_coef' );
2204
if(( ior == 13 ) .or. ( ior == 17 )) then
2205
call Array_Matrix3_Set( unit, ior, ncomp, ncomp_diff_coef, nphase, len_name, fcn_name, comp_diff_coef )
2207
comp_diff_coef = value_real
2210
Case( 'capil_pres_coef' );
2211
if(( ior == 13 ) .or. ( ior == 17 )) then
2212
call Array_Matrix3_Set( unit, ior, ncapil_pres_coef, nphase, nphase, len_name, fcn_name, capil_pres_coef )
2214
capil_pres_coef = value_real
2217
Case( 'u_abs_stab' );
2218
if(( ior == 13 ) .or. ( ior == 17 )) then
2219
call Array_Matrix3_Set( unit, ior, mat_nonods, ndim * nphase, ndim * nphase , len_name, fcn_name, u_abs_stab )
2221
u_abs_stab = value_real
2225
if(( ior == 13 ) .or. ( ior == 17 )) then
2226
call Array_Matrix3_Set( unit, ior, mat_nonods, ndim * nphase, ndim * nphase , len_name, fcn_name, u_absorb )
2228
u_absorb = value_real
2232
if(( ior == 13 ) .or. ( ior == 17 )) then
2233
call Array_Matrix3_Set( unit, ior, cv_nonods, nphase, nphase, len_name, fcn_name, t_absorb )
2235
t_absorb = value_real
2239
if(( ior == 13 ) .or. ( ior == 17 )) then
2240
call Array_Matrix3_Set( unit, ior, cv_nonods, nphase, nphase, len_name, fcn_name, v_absorb )
2242
v_absorb = value_real
2246
if(( ior == 13 ) .or. ( ior == 17 )) then
2247
call Array_Matrix3_Set( unit, ior, totele, ndim, ndim, len_name, fcn_name, perm )
2252
! Matrices and Tensors ( x , y , z , w)
2253
Case( 'udiffusion' );
2254
if(( ior == 13 ) .or. ( ior == 18 )) then
2255
call Array_Matrix4_Set( unit, ior, mat_nonods, ndim, ndim, nphase, len_name, fcn_name, udiffusion )
2257
udiffusion = value_real
2260
Case( 'tdiffusion' );
2261
if(( ior == 13 ) .or. ( ior == 18 )) then
2262
call Array_Matrix4_Set( unit, ior, mat_nonods, ndim, ndim, nphase, len_name, fcn_name, tdiffusion )
2264
tdiffusion = value_real
2267
Case( 'comp_diffusion' );
2268
if(( ior == 13 ) .or. ( ior == 18 )) then
2269
call Array_Matrix4_Set( unit, ior, mat_nonods, ndim, ndim, nphase, len_name, fcn_name, comp_diffusion )
2271
comp_diffusion = value_real
2275
write( 375, * ) 'Option not found - read_all subrt., ', ifile( 1 : k )
2280
end if Conditional_IOR
2284
write(357,*) 'Leaving read_all'
2287
end subroutine read_all
2290
subroutine readin_bc( unit_input, nphase, stotel, cv_snloc, p_snloc, u_snloc, &
2291
wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc, &
2292
suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc, suf_p_bc, &
2293
suf_u_bc, suf_v_bc, suf_w_bc, suf_one_bc, &
2294
suf_u_bc_rob1, suf_u_bc_rob2, suf_v_bc_rob1, suf_v_bc_rob2, &
2295
suf_w_bc_rob1, suf_w_bc_rob2, suf_t_bc_rob1, suf_t_bc_rob2 )
2297
integer, intent( in ) :: unit_input
2298
integer, intent( in ) :: nphase, stotel, cv_snloc, p_snloc, u_snloc
2299
integer, dimension( stotel * nphase ), intent( inout ) :: wic_vol_bc, wic_d_bc, wic_u_bc, wic_p_bc, wic_t_bc
2300
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_vol_bc, suf_d_bc, suf_cpd_bc, suf_t_bc
2301
real, dimension( stotel * p_snloc * nphase ), intent( inout ) :: suf_p_bc
2302
real, dimension( stotel * u_snloc * nphase ), intent( inout ) :: suf_u_bc, suf_v_bc, suf_w_bc
2303
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_one_bc
2304
real, dimension( stotel * u_snloc * nphase ), intent( inout ) :: suf_u_bc_rob1, suf_u_bc_rob2, suf_v_bc_rob1, suf_v_bc_rob2, &
2305
suf_w_bc_rob1, suf_w_bc_rob2
2306
real, dimension( stotel * cv_snloc * nphase ), intent( inout ) :: suf_t_bc_rob1, suf_t_bc_rob2
2309
character( len = 150) :: title
2311
write(357,*) 'In readin_bc'
2313
read( unit_input, * ) title( 1 : 150 )
2314
read( unit_input, * ) title( 1 : 150 )
2315
read( unit_input, * ) title( 1 : 150 )
2316
read( unit_input, * ) title( 1 : 150 )
2317
call input_int( unit_input, stotel * nphase, wic_vol_bc )
2318
call input_int( unit_input, stotel * nphase, wic_d_bc )
2319
call input_int( unit_input, stotel * nphase, wic_u_bc )
2320
call input_int( unit_input, stotel * nphase, wic_p_bc )
2321
call input_int( unit_input, stotel * nphase, wic_t_bc )
2323
call input_real( unit_input, stotel * cv_snloc * nphase, suf_vol_bc )
2324
call input_real( unit_input, stotel * cv_snloc * nphase, suf_d_bc )
2325
call input_real( unit_input, stotel * cv_snloc * nphase, suf_cpd_bc )
2326
call input_real( unit_input, stotel * cv_snloc * nphase, suf_t_bc )
2328
call input_real( unit_input, stotel * p_snloc * nphase, suf_p_bc )
2330
call input_real( unit_input, stotel * u_snloc * nphase, suf_u_bc )
2331
call input_real( unit_input, stotel * u_snloc * nphase, suf_v_bc )
2332
call input_real( unit_input, stotel * u_snloc * nphase, suf_w_bc )
2334
call input_real( unit_input, stotel * cv_snloc * nphase, suf_one_bc )
2336
call input_real( unit_input, stotel * u_snloc * nphase, suf_u_bc_rob1 )
2337
call input_real( unit_input, stotel * u_snloc * nphase, suf_u_bc_rob2 )
2338
call input_real( unit_input, stotel * u_snloc * nphase, suf_v_bc_rob1 )
2339
call input_real( unit_input, stotel * u_snloc * nphase, suf_v_bc_rob2 )
2340
call input_real( unit_input, stotel * u_snloc * nphase, suf_w_bc_rob1 )
2341
call input_real( unit_input, stotel * u_snloc * nphase, suf_w_bc_rob2 )
2342
call input_real( unit_input, stotel * cv_snloc * nphase, suf_t_bc_rob1 )
2343
call input_real( unit_input, stotel * cv_snloc * nphase, suf_t_bc_rob2 )
2345
write(357,*) 'Leaving readin_bc'
2348
end subroutine readin_bc
2350
subroutine reading_initial_position_velocities( unit_input, nphase, &
2351
x_nonods, xu_nonods, u_nonods, &
2352
x, y, z, xu, yu, zu, nu, nv, nw, ug, vg, wg )
2355
integer, intent( in ) :: unit_input, nphase, x_nonods, xu_nonods,u_nonods
2356
real, dimension( x_nonods ), intent( inout ) :: x, y, z
2357
real, dimension( xu_nonods ), intent( inout ) :: xu, yu, zu
2358
real, dimension( u_nonods * nphase ), intent( inout ) :: nu, nv, nw, ug, vg, wg
2360
write(357,*) 'In reading_initial_position_velocities'
2362
call input_real( unit_input, x_nonods, x )
2363
call input_real( unit_input, x_nonods, y )
2364
call input_real( unit_input, x_nonods, z )
2366
call input_real( unit_input, xu_nonods, xu )
2367
call input_real( unit_input, xu_nonods, yu )
2368
call input_real( unit_input, xu_nonods, zu )
2370
call input_real( unit_input, u_nonods * nphase, nu )
2371
call input_real( unit_input, u_nonods * nphase, nv )
2372
call input_real( unit_input, u_nonods * nphase, nw )
2374
call input_real( unit_input, u_nonods * nphase, ug )
2375
call input_real( unit_input, u_nonods * nphase, vg )
2376
call input_real( unit_input, u_nonods * nphase, wg )
2378
write(357,*) 'Leaving reading_initial_position_velocities'
2381
end subroutine reading_initial_position_velocities
2384
subroutine input_int( unit, n, array )
2386
integer, intent( in ) :: unit, n
2387
integer, dimension( n ), intent( inout ) :: array
2390
character( len = 100 ) :: description, fcn_name
2391
integer :: value, i, k, npos, length
2392
integer, dimension( : ), allocatable :: ind
2394
read( unit, 101 ) description
2395
read( unit, * ) value
2396
if( value < -1000 ) then ! Pre-defined function that can be used
2397
read( unit, * ) fcn_name
2398
k = index( fcn_name, ' ' ) - 1
2399
call system( 'rm -f fvalues' )
2400
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
2401
open( unit + 1, file = 'fvalues', status = 'unknown' )
2402
read( unit + 1, * ) length
2404
read( unit + 1 , * ) array( i )
2408
array( 1 : n ) = value
2409
! For specific values in the array, specify number of index (if negative, it is assumed
2410
! the array is filled with one single value as described by the previous 'read' statement.
2413
! 1 6 5 8 --> two index to be modified, 1 and 6, and array in these positions assumed
2415
read( unit, * ) npos
2417
print*, 'number of index is larger the size of the array'
2419
elseif( npos < 0 ) then
2422
allocate( ind( npos ))
2423
read( unit, * ) (ind( i ), i = 1, npos ), ( array ( ind( i )), i = 1, npos )
2430
end subroutine input_int
2432
subroutine input_real( unit, n, array )
2434
integer, intent( in ) :: unit, n
2435
real, dimension( n ), intent( inout ) :: array
2438
character( len = 100 ) :: description, fcn_name
2440
integer :: i, npos, k, length
2441
integer, dimension( : ), allocatable :: ind
2443
read( unit, 101 ) description
2444
read( unit, * ) value
2445
if( value < -1000. ) then ! pre-defined function that can be used
2446
read( unit, * ) fcn_name
2447
k = index( fcn_name, ' ' ) - 1
2448
call system( 'rm -f fvalues' )
2449
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
2450
open( unit + 1, file = 'fvalues', status = 'unknown' )
2451
read( unit + 1, * ) length
2453
read( unit + 1 , * ) array( i )
2457
array( 1 : n ) = value
2458
! For specific values in the array, specify number of index (if negative, it is assumed
2459
! the array is filled with one single value as described by the previous 'read' statement.
2462
! 1 6 5 8 --> two index to be modified, 1 and 6, and array in these positions assumed
2464
read( unit, * ) npos
2466
print*, 'number of index is larger the size of the array'
2468
elseif( npos < 0 ) then
2471
allocate( ind( npos ))
2472
read( unit, * ) (ind( i ), i = 1, npos ), ( array ( ind( i )), i = 1, npos )
2479
end subroutine input_real
2481
subroutine input_matrix( unit, ndim1, ndim2, array )
2482
! Here, user can either provide a single value for all components of
2483
! the tensor or use a external function to design the tensor with
2484
! input array( i, j, k )
2486
integer, intent( in ) :: unit, ndim1, ndim2
2487
real, dimension( ndim1, ndim2 ), intent( inout ) :: array
2489
character( len = 100 ) :: description, fcn_name
2490
integer :: i, j, k, length1, length2
2493
read( unit, 101 ) description
2494
read( unit, * ) value
2495
if( value < -1000. ) then ! pre-defined function that can be used
2496
read( unit, * ) fcn_name
2497
k = index( fcn_name, ' ' ) - 1
2498
call system( 'rm -f fvalues' )
2499
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
2500
open( unit + 1, file = 'fvalues', status = 'unknown' )
2501
read( unit + 1, * ) length1, length2
2502
if( ( length1 > ndim1 ) .or. ( length2 > ndim2 ))then
2503
print*, 'Matrix dimension exceeded'
2508
read( unit + 1 , * ) array( i, j )
2513
array( 1 : ndim1, 1 : ndim2 ) = value
2519
end subroutine input_matrix
2522
subroutine input_tensor( unit, ndim1, ndim2, ndim3, array )
2523
! Here, user can either provide a single value for all components of
2524
! the tensor or use a external function to design the tensor with
2525
! input array( i, j, k )
2527
integer, intent( in ) :: unit, ndim1, ndim2, ndim3
2528
real, dimension( ndim1, ndim2, ndim3 ), intent( inout ) :: array
2530
character( len = 100 ) :: description, fcn_name
2531
integer :: i, j, k, length1, length2, length3
2534
read( unit, 102 ) description
2535
read( unit, * ) value
2536
if( value < -1000. ) then ! pre-defined function that can be used
2537
read( unit, * ) fcn_name
2538
k = index( fcn_name, ' ' ) - 1
2539
call system( 'rm -f fvalues' )
2540
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
2541
open( unit + 1, file = 'fvalues', status = 'unknown' )
2542
read( unit + 1, * ) length1, length2, length3
2543
if( ( length1 > ndim1 ) .or. ( length2 > ndim2 ) .or. &
2544
( length3 > ndim3 ))then
2545
print*, 'Tensor dimension exceeded'
2551
read( unit + 1 , * ) array( i, j, k )
2557
array( 1 : ndim1, 1 : ndim2, 1 : ndim3 ) = value
2563
end subroutine input_tensor
2567
subroutine input_tensor2( unit, ndim1, ndim2, ndim3, ndim4, array )
2568
! Here, user can either provide a single value for all components of
2569
! the tensor or use a external function to design the tensor with
2570
! input array( i, j, k )
2572
integer, intent( in ) :: unit, ndim1, ndim2, ndim3, ndim4
2573
real, dimension( ndim1, ndim2, ndim3, ndim4 ), intent( inout ) :: array
2575
character( len = 100 ) :: description, fcn_name
2576
integer :: i, j, k, l, length1, length2, length3, length4
2579
read( unit, 101 ) description
2580
read( unit, * ) value
2581
if( value < -1000. ) then ! pre-defined function that can be used
2582
read( unit, * ) fcn_name
2583
k = index( fcn_name, ' ' ) - 1
2584
call system( 'rm -f fvalues' )
2585
call system( './' // fcn_name( 1 : k ) // ' > fvalues ' )
2586
open( unit + 1, file = 'fvalues', status = 'unknown' )
2587
read( unit + 1, * ) length1, length2, length3
2588
if( ( length1 > ndim1 ) .or. ( length2 > ndim2 ) .or. &
2589
( length3 > ndim3 ) .or. ( length4 > ndim4 ))then
2590
print*, 'Tensor dimension exceeded'
2597
read( unit + 1 , * ) array( i, j, k, l )
2604
array( 1 : ndim1, 1 : ndim2, 1 : ndim3, ndim4 ) = value
2610
end subroutine input_tensor2
2613
subroutine allocating_global_nodes( ndim, totele, domain_length, &
2614
u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, mat_nloc, &
2615
cv_snloc, u_snloc, p_snloc, stotel, &
2616
cv_nonods, u_nonods, x_nonods, xu_nonods, &
2617
u_ele_type, cv_ele_type, &
2619
u_ndgln, xu_ndgln, cv_ndgln, x_ndgln, p_ndgln, &
2620
mat_ndgln, u_sndgln, cv_sndgln, p_sndgln )
2624
integer, intent( in ) :: ndim, totele
2625
real, intent( in ) :: domain_length
2626
integer, intent( in ) :: u_nloc, xu_nloc, cv_nloc, x_nloc, p_nloc, mat_nloc, &
2627
cv_snloc, u_snloc, p_snloc, stotel, &
2628
cv_nonods, u_nonods, x_nonods, xu_nonods, &
2629
u_ele_type, cv_ele_type
2630
real, dimension( x_nonods ), intent( inout ) :: x
2631
real, dimension( xu_nonods ), intent( inout ) :: xu
2632
integer, dimension( totele * u_nloc ), intent( inout ) :: u_ndgln
2633
integer, dimension( totele * xu_nloc ), intent( inout ) :: xu_ndgln
2634
integer, dimension( totele * cv_nloc ), intent( inout ) :: cv_ndgln
2635
integer, dimension( totele * cv_nloc ), intent( inout ) :: x_ndgln
2636
integer, dimension( totele * p_nloc ), intent( inout ) :: p_ndgln
2637
integer, dimension( totele * mat_nloc ), intent( inout ) :: mat_ndgln
2638
integer, dimension( stotel * u_snloc ), intent( inout ) :: u_sndgln
2639
integer, dimension( stotel * cv_snloc ), intent( inout ) :: cv_sndgln
2640
integer, dimension( stotel * p_snloc ), intent( inout ) :: p_sndgln
2643
integer :: u_nloc2, u_nod, xu_nod, cv_nod, x_nod, mat_nod, p_nod
2645
integer :: ele, iloc, cv_iloc
2646
character( len = 100 ) :: name
2648
write(357,*) 'In allocating_global_nodes'
2650
Conditional_NDIM: if( ndim == 1 ) then ! This needs to be updated for 2-3D
2652
dx = domain_length / real( totele )
2654
cv_sndgln( 2 ) = cv_nonods
2656
p_sndgln( 2 ) = cv_nonods
2658
if( cv_ele_type == 2 ) then
2659
u_nloc2 = u_nloc / cv_nloc
2660
do cv_iloc = 1, cv_nloc
2661
u_sndgln( cv_iloc ) = 1 + ( cv_iloc -1 ) * u_nloc2
2662
u_sndgln( cv_iloc + cv_nloc ) = u_nonods - u_nloc + cv_iloc * u_nloc2
2666
u_sndgln( 2 ) = u_nonods
2668
name = '####u_sndgln####'
2669
call mirror_array_int( 357, name, stotel * u_snloc, u_sndgln )
2678
Loop_Elements: do ele = 1, totele
2680
Loop_U: do iloc = 1, u_nloc ! storing velocity nodes
2682
u_ndgln( ( ele - 1 ) * u_nloc + iloc ) = u_nod
2685
Loop_XU:do iloc = 1, xu_nloc
2687
xu_ndgln( ( ele - 1 ) * xu_nloc + iloc ) = xu_nod
2688
if ( xu_nloc == 1 ) then
2689
xu( xu_nod ) = ( real( ele - 1 ) + 0.5 ) * dx
2691
if( u_ele_type == 2 ) then
2692
xu( xu_nod ) = real( ele - 1 ) * dx + real ( mod( iloc, cv_nloc ) - 1 ) &
2693
* dx / real ( u_nloc - 1 )
2695
xu( xu_nod ) = real( ele - 1 ) * dx + real ( iloc - 1 ) &
2696
* dx / real ( u_nloc - 1 )
2700
if( xu_nloc /= 1 ) xu_nod = xu_nod - 1
2702
Loop_P: do iloc = 1, p_nloc
2704
p_ndgln( ( ele - 1 ) * p_nloc + iloc ) = p_nod
2706
!if( problem == 2 ) p_nod = p_nod - 1
2707
if( cv_nonods /= totele * cv_nloc ) p_nod = p_nod - 1
2709
Loop_CV: do iloc = 1, cv_nloc
2711
cv_ndgln( ( ele - 1 ) * cv_nloc + iloc ) = cv_nod
2713
!if( problem == 2 ) cv_nod = cv_nod - 1
2714
if( cv_nonods /= totele * cv_nloc ) cv_nod = cv_nod - 1
2716
Loop_X: do iloc = 1, x_nloc
2718
x_ndgln( ( ele - 1 ) * x_nloc + iloc ) = x_nod
2719
if ( x_nloc == 1 ) then
2720
x( x_nod ) = ( real( ele - 1 ) + 0.5 ) * dx
2722
x( x_nod ) = real( ele - 1 ) * dx + real( iloc - 1 ) * dx / real ( x_nloc - 1 )
2725
if( x_nloc /= 1 ) x_nod = x_nod - 1
2727
Loop_Mat: do iloc = 1, mat_nloc
2728
mat_nod = mat_nod + 1
2729
mat_ndgln( ( ele - 1 ) * mat_nloc + iloc ) = mat_nod
2732
end do Loop_Elements
2734
end if Conditional_NDIM
2736
write(357,*) 'Leaving allocating_global_nodes'
2740
end subroutine allocating_global_nodes
2743
! Initialising T and Told:
2744
subroutine initialise_scalar_fields( &
2745
problem, ndim, nphase, totele, domain_length, &
2746
x_nloc, cv_nloc, x_nonods, cv_nonods, &
2747
x_ndgln, cv_ndgln, &
2750
integer, intent( in ) :: problem, ndim, nphase, totele
2751
real, intent( in ) :: domain_length
2752
integer, intent( in ) ::x_nloc, cv_nloc, x_nonods, cv_nonods
2753
integer, dimension( totele * cv_nloc ) :: x_ndgln, cv_ndgln
2754
real, dimension( x_nonods ), intent( in ) :: x
2755
real, dimension( cv_nonods * nphase ), intent( inout ) :: told, t
2757
integer :: iloc, ele, x_nod, cv_nod
2760
write(357,*) 'In initialise_scalar_fields'
2762
Conditional_NDIM: if ( ndim == 1 ) then ! This may change for 2-3D
2764
dx = domain_length / real( totele )
2766
Loop_Element: do ele = 1, totele
2767
do iloc = 1, cv_nloc
2768
x_nod = x_ndgln( ( ele - 1 ) * x_nloc + iloc )
2769
cv_nod = cv_ndgln( ( ele - 1 ) * cv_nloc + iloc )
2771
Select case( problem )
2773
case( -1, 0 ); ! CV-Adv test-case
2774
if( ( x( x_nod ) > 0.199 ) .and. ( x( x_nod ) < 0.401 )) t( cv_nod ) = 1.
2775
if( problem == 0 ) t( cv_nod ) = exp( -(( x( x_nod ) - 0.6 ) / 0.2 )**2 )
2777
case( 1, 2 ); ! BL problem
2778
if( x( x_nod ) < dx ) t( cv_nod ) = 0.5
2781
write(357,*) 'It is necessary to set up -1 < problem < 1'
2790
end if Conditional_NDIM
2792
write(357,*) 'Leaving initialise_scalar_fields'
2795
end subroutine initialise_scalar_fields
2798
integer function Combination( n, r )
2799
! This function performs the combinatorial:
2800
! C(n,r) = n! / ( (n-r)! r! )
2804
Combination = Permut( n ) / ( Permut( n - r ) * Permut( r ))
2807
end function Combination
2809
integer function Permut( n )
2810
! This function performs probabilistic permutation:
2824
end module input_var
added
removed
multiphase/src/multi_IO.F90
