Introduction
This page provides an example of code blocks that may be used to perform a calculation with drift scalars for coal combustion.
Physical properties
Local variables to be added
The following local variables need to be defined for the examples in this section:
integer ivart, iel, ifac
integer ivar
integer f_id
integer icla
integer iscdri, keydri, iflid, nfld, keyccl
double precision xvart
double precision aa1, bb1, cc1, dd1
double precision aa2, bb2, cc2, dd2
double precision aa3, bb3, cc3, dd3
double precision aa4, bb4, cc4, dd4
double precision aa5, bb5, cc5, dd5
double precision aa6, bb6, cc6, dd6
double precision aa7, bb7, cc7, dd7
double precision visco_O2, visco_CO, visco_H2, visco_N2
double precision visco_SO2, visco_NH3, visco_CO2
character*80 fname
double precision, allocatable, dimension(:) :: visco
double precision, dimension(:), pointer :: cpro_rom1, cpro_rom2, cpro_diam2
double precision, dimension(:), pointer :: cpro_temp1, cpro_x2, cpro_x1
double precision, dimension(:), pointer :: cpro_ym1_3, cpro_ym1_5, cpro_ym1_7
double precision, dimension(:), pointer :: cpro_ym1_8
double precision, dimension(:), pointer :: cpro_ym1_9, cpro_ym1_11, cpro_ym1_12
double precision, dimension(:), pointer :: cpro_taup
double precision, dimension(:), pointer :: cpro_taupg
Initialization and finalization
The following initialization block needs to be added for the following examples:
allocate(visco(ncelet))
call field_get_val_s(iym1(3), cpro_ym1_3)
call field_get_val_s(iym1(5), cpro_ym1_5)
call field_get_val_s(iym1(7), cpro_ym1_7)
call field_get_val_s(iym1(8), cpro_ym1_8)
call field_get_val_s(iym1(9), cpro_ym1_9)
call field_get_val_s(iym1(11), cpro_ym1_11)
call field_get_val_s(iym1(12), cpro_ym1_12)
call field_get_key_id("drift_scalar_model", keydri)
call field_get_key_id("scalar_class", keyccl)
call field_get_n_fields(nfld)
In theory Fortran 95 deallocates locally-allocated arrays automatically, but deallocating arrays in a symmetric manner to their allocation is good practice, and it avoids using a different logic for C and Fortran.
Body
Here is the corresponding code:
call field_get_val_s(itemp1, cpro_temp1)
call field_get_val_s(irom1, cpro_rom1)
if (ntcabs.le.1) then
do iel = 1, ncel
visco(iel) = viscl0
cpro_rom1(iel) = ro0
enddo
do icla = 1, nclacp
call field_get_val_s(irom2(icla), cpro_rom2)
call field_get_val_s(idiam2(icla), cpro_diam2)
do iel = 1, ncel
cpro_rom2(iel) = rho20(icla)
cpro_diam2(iel) = diam20(icla)
enddo
enddo
endif
aa1 = 4.0495d-6
bb1 = 6.22d-8
cc1 = -2.3032d-11
dd1 = 4.4077d-15
aa2 = 9.9987d-6
bb2 = 5.1578d-8
cc2 = -1.8383d-11
dd2 = 3.33307d-15
aa3 = 2.894d-6
bb3 = 2.22508d-8
cc3 = -8.041d-12
dd3 = 1.4619d-15
aa4 = 4.3093d-6
bb4 = 5.0516d-8
cc4 = -1.7869d-11
dd4 = 3.2136d-15
aa5 = -1.9889d-6
bb5 = 5.365d-8
cc5 = -1.4286d-11
dd5 = 2.1639d-15
aa6 = -1.293d-6
bb6 = 4.1194d-8
cc6 = -1.772d-11
dd6 = 1.8699d-15
aa7 = 4.4822d-7
bb7 = 5.4327d-8
cc7 = -1.7581d-11
dd7 = 2.9979d-15
if (ntcabs.gt.1) then
do iel = 1, ncel
xvart = cpro_temp1(iel)
visco_o2 = aa1 + xvart*bb1 + cc1*xvart**2 + dd1*xvart**3
visco_co = aa2 + xvart*bb2 + cc2*xvart**2 + dd2*xvart**3
visco_h2 = aa3 + xvart*bb3 + cc3*xvart**2 + dd3*xvart**3
visco_n2 = aa4 + xvart*bb4 + cc4*xvart**2 + dd4*xvart**3
visco_so2 = aa5 + xvart*bb5 + cc5*xvart**2 + dd5*xvart**3
visco_nh3 = aa6 + xvart*bb6 + cc6*xvart**2 + dd6*xvart**3
visco_co2 = aa7 + xvart*bb7 + cc7*xvart**2 + dd7*xvart**3
visco(iel) = ( cpro_ym1_8(iel) * visco_o2 &
+ cpro_ym1_3(iel) * visco_co &
+ cpro_ym1_5(iel) * visco_h2 &
+ cpro_ym1_12(iel)* visco_n2 &
+ cpro_ym1_11(iel)* visco_so2 &
+ cpro_ym1_7(iel) * visco_nh3 &
+ cpro_ym1_9(iel) * visco_co2 )/ &
( cpro_ym1_8(iel) + cpro_ym1_3(iel) &
+ cpro_ym1_5(iel) + cpro_ym1_12(iel) &
+ cpro_ym1_11(iel)+ cpro_ym1_7(iel) &
+ cpro_ym1_9(iel))
enddo
endif
call field_get_val_s_by_name("x_c", cpro_x1)
do iflid = 0, nfld-1
call field_get_key_int(iflid, keyccl, icla)
call field_get_key_int(iflid, keydri, iscdri)
if (icla.le.-1.and.btest(iscdri, drift_scalar_add_drift_flux)) then
call field_get_name(iflid, fname)
call field_get_id('drift_tau_'//trim(fname), f_id)
call field_get_val_s(f_id, cpro_taupg)
do iel = 1, ncel
cpro_taupg(iel) = 0.d0
enddo
endif
enddo
do iflid = 0, nfld-1
call field_get_key_int(iflid, keyccl, icla)
call field_get_key_int(iflid, keydri, iscdri)
if (icla.ge.1.and.btest(iscdri, drift_scalar_add_drift_flux)) then
call field_get_val_s(irom2(icla), cpro_rom2)
call field_get_val_s(idiam2(icla), cpro_diam2)
call field_get_val_s(ix2(icla), cpro_x2)
call field_get_name(iflid, fname)
call field_get_id('drift_tau_'//trim(fname), f_id)
call field_get_val_s(f_id, cpro_taup)
do iel = 1, ncel
cpro_taup(iel) = cpro_x1(iel) * cpro_rom2(iel) &
* cpro_diam2(iel)**2 &
/ (18.d0*visco(iel))
enddo
do iel = 1, ncel
cpro_taupg(iel) = cpro_taupg(iel) &
- ( cpro_taup(iel) * cpro_x2(iel) )
enddo
endif
enddo