From 25c5d0a80f1894d58bd4dc65d02ce92fd50803de Mon Sep 17 00:00:00 2001 From: fangjian Date: Thu, 2 Nov 2023 17:22:04 +0000 Subject: [PATCH] clean nondimen and COMB in fludyna --- src/commvar.F90 | 4 +- src/fludyna.F90 | 509 ++++++++++++++++++++++++---------------------- src/solver.F90 | 16 +- src/thermchem.F90 | 4 +- 4 files changed, 285 insertions(+), 248 deletions(-) diff --git a/src/commvar.F90 b/src/commvar.F90 index 56fc63f..ef15db1 100644 --- a/src/commvar.F90 +++ b/src/commvar.F90 @@ -90,7 +90,7 @@ module commvar integer,allocatable :: imon(:,:) real(8) :: time,deltat real(8) :: uinf,vinf,winf,roinf,pinf,tinf - real(8) :: ref_t,reynolds,mach,rgas,gamma,prandtl + real(8) :: ref_t,reynolds,mach,rgas,cp,cv,gamma,prandtl real(8),allocatable :: schmidt(:) real(8) :: const1,const2,const3,const4,const5,const6,const7 real(8) :: tempconst,tempconst1 @@ -109,6 +109,8 @@ module commvar !| reynolds | Reynolds number. | !| mach | Mach number. | !| rgas | gas constant, p=ro*rgas*T. | + !| cp | specific heat at constant pressure. | + !| cv | specific heat at constant volume. | !| gamma | specific heat ratio. | !| prandtl | Prandtl number. | !| schmidt | the Schmidt number array for each scalar. | diff --git a/src/fludyna.F90 b/src/fludyna.F90 index 9cd726d..41e3e59 100644 --- a/src/fludyna.F90 +++ b/src/fludyna.F90 @@ -50,19 +50,8 @@ function thermal_scar(density,pressure,temperature,species) result(vout) real(8) :: vout real(8),intent(in) ,optional :: density,pressure,temperature,species(:) ! -#ifdef COMB - - if(present(density).and.present(temperature)) then - vout = density*temperature*rgcmix(species) - elseif(present(density).and.present(pressure)) then - vout = pressure/density/rgcmix(species) - elseif(present(temperature).and.present(pressure)) then - vout = pressure/temperature/rgcmix(species) - else - stop ' !! unable to use dimensional EoS @ thermal !!' - endif + real(8) :: rloc -#else if(nondimen) then ! @@ -78,25 +67,29 @@ function thermal_scar(density,pressure,temperature,species) result(vout) ! else ! +#ifdef COMB + rloc=rgcmix(species) +#else + rloc=rgas +#endif + ! if(present(density).and.present(temperature)) then - vout = density*temperature*rgas + vout = density*temperature*rloc elseif(present(density).and.present(pressure)) then - vout = pressure/(density*rgas) + vout = pressure/density/rloc elseif(present(temperature).and.present(pressure)) then - vout = pressure/(temperature*rgas) + vout = pressure/temperature/rloc else stop ' !! unable to use dimensional EoS @ thermal !!' endif ! endif -#endif - ! end function thermal_scar ! function thermal_1d(density,pressure,temperature,species,dim) result(vout) ! - use commvar,only : const2 + use commvar,only : const2,rgas ! ! arguments integer,intent(in) :: dim @@ -104,6 +97,8 @@ function thermal_1d(density,pressure,temperature,species,dim) result(vout) real(8),intent(in),optional :: density(:),pressure(:), & temperature(:),species(:,:) ! + real(8) :: rloc(dim) + ! if(nondimen) then ! if(present(density) .and. present(temperature)) then @@ -118,12 +113,18 @@ function thermal_1d(density,pressure,temperature,species,dim) result(vout) ! else ! +#ifdef COMB + rloc=rgcmix(species,dim) +#else + rloc=rgas +#endif + if(present(density).and.present(temperature)) then - vout = density*temperature*rgcmix(species,dim) + vout = density*temperature*rloc elseif(present(density).and.present(pressure)) then - vout = pressure/density/rgcmix(species,dim) + vout = pressure/density/rloc elseif(present(temperature).and.present(pressure)) then - vout = pressure/temperature/rgcmix(species,dim) + vout = pressure/temperature/rloc else stop ' !! unable to use dimensional EoS @ thermal !!' endif @@ -134,7 +135,7 @@ end function thermal_1d ! function thermal_3d(density,pressure,temperature,species,dim) result(vout) ! - use commvar,only : const2 + use commvar,only : const2,rgas ! ! arguments integer,intent(in) :: dim(3) @@ -142,6 +143,8 @@ function thermal_3d(density,pressure,temperature,species,dim) result(vout) real(8),intent(in),optional :: density(:,:,:),pressure(:,:,:), & temperature(:,:,:),species(:,:,:,:) ! + real(8) :: rloc(dim(1),dim(2),dim(3)) + ! if(nondimen) then ! if(present(density) .and. present(temperature)) then @@ -151,19 +154,24 @@ function thermal_3d(density,pressure,temperature,species,dim) result(vout) elseif(present(temperature) .and. present(pressure)) then vout=pressure/temperature*const2 else - stop ' !! unable to get thermal variable @ thermal !!' + stop ' !! unable to get thermal variable @ thermal_3d !!' endif ! else ! +#ifdef COMB + rloc=rgcmix(species,dim) +#else + rloc=rgas +#endif if(present(density).and.present(temperature)) then - vout = density*temperature*rgcmix(species,dim) + vout = density*temperature*rloc elseif(present(density).and.present(pressure)) then - vout = pressure/density/rgcmix(species,dim) + vout = pressure/density/rloc elseif(present(temperature).and.present(pressure)) then - vout = pressure/temperature/rgcmix(species,dim) + vout = pressure/temperature/rloc else - stop ' !! unable to use dimensional EoS @ thermal !!' + stop ' !! unable to use dimensional EoS @ thermal_3d !!' endif ! endif @@ -304,7 +312,7 @@ end subroutine updateq subroutine fvar2q_sca(q,density,velocity,pressure,temperature, & species,tke,omega) ! - use commvar, only: numq,ndims,num_species,const1,const6 + use commvar, only: numq,ndims,num_species,const1,const6,cv ! real(8),intent(in) :: density,velocity(:) real(8),intent(in),optional :: pressure,temperature,species(:),tke,omega @@ -312,43 +320,43 @@ subroutine fvar2q_sca(q,density,velocity,pressure,temperature, & ! ! local data integer :: jspec,j - real(8) :: var1,var2 + real(8) :: var1,var2,cotem ! q(1)=density q(2)=density*velocity(1) q(3)=density*velocity(2) q(4)=density*velocity(3) ! - if(nondimen) then - ! - if(present(temperature)) then - q(5)=density*( temperature*const1 + 0.5d0*(velocity(1)**2 + & - velocity(2)**2 + & - velocity(3)**2) ) - elseif(present(pressure)) then - q(5)=pressure*const6+0.5d0*density*( velocity(1)**2 + & - velocity(2)**2 + & - velocity(3)**2 ) - else - print*,' !! pressure or temperature required for energy calculation !!' - stop ' !! error @ fvar2q' - endif +#ifdef COMB + if(.not.present(species)) & + stop ' !! error @ fvar2q - species not set for energy calculation !!' ! + if(present(temperature)) then + var1=0.5d0*sum(velocity(:)*velocity(:)) + call cpeval(tmp=temperature,spc=species,ke=var1,eng=var2) + q(5)=var2*density else - ! - if(.not.present(species)) & - stop ' !! error @ fvar2q - species not set for energy calculation !!' - ! - if(present(temperature)) then - var1=0.5d0*sum(velocity(:)*velocity(:)) - call cpeval(tmp=temperature,spc=species,ke=var1,eng=var2) - q(5)=var2*density - else - print*,' !! temperature required for energy calculation !!' - stop ' !! error @ fvar2q' - endif - ! - endif !nondimen + print*,' !! temperature required for energy calculation !!' + stop ' !! error @ fvar2q' + endif + ! +#else + if(nondimen) then + cotem=const1 + else + cotem=cv + endif + ! + var1=0.5d0*sum(velocity(:)*velocity(:)) + if(present(temperature)) then + q(5)=density*(temperature*cotem+var1) + elseif(present(pressure)) then + q(5)=pressure*const6+density*var1 + else + print*,' !! pressure or temperature required for energy calculation !!' + stop ' !! error @ fvar2q' + endif +#endif ! if(num_species>0) then ! @@ -369,7 +377,7 @@ end subroutine fvar2q_sca ! subroutine fvar2q_1da(q,density,velocity,pressure,temperature,species) ! - use commvar, only: numq,ndims,num_species,const1,const6 + use commvar, only: numq,ndims,num_species,const1,const6,cv ! real(8),intent(in) :: density(:),velocity(:,:) real(8),intent(in),optional :: pressure(:),temperature(:), & @@ -378,49 +386,55 @@ subroutine fvar2q_1da(q,density,velocity,pressure,temperature,species) ! ! local data integer :: jspec,j - real(8) :: var1,var2 + real(8) :: var1,var2,cotem ! q(:,1)=density(:) q(:,2)=density(:)*velocity(:,1) q(:,3)=density(:)*velocity(:,2) q(:,4)=density(:)*velocity(:,3) ! - if(nondimen) then +#ifdef COMB + + if(.not.present(species)) & + stop ' !! error @ fvar2q - species not set for energy calculation !!' + ! + if(present(temperature)) then ! - if(present(temperature)) then - q(:,5)=density(:)*( temperature(:)*const1 + & - 0.5d0*(velocity(:,1)**2 + & - velocity(:,2)**2 + & - velocity(:,3)**2) ) - elseif(present(pressure)) then - q(:,5)=pressure(:)*const6+0.5d0*density(:)*( & - velocity(:,1)**2 + & - velocity(:,2)**2 + & - velocity(:,3)**2 ) - else - print*,' !! pressure or temperature required !!' - stop ' !! error @ fvar2q' - endif + do j=1,size(density) + var1=0.5d0*sum(velocity(j,:)*velocity(j,:)) + call cpeval(tmp=temperature(j),spc=species(j,:),ke=var1,eng=var2) + q(j,5)=var2*density(j) + enddo ! else - ! - if(.not.present(species)) & - stop ' !! error @ fvar2q - species not set for energy calculation !!' - ! - if(present(temperature)) then - ! - do j=1,size(density) - var1=0.5d0*sum(velocity(j,:)*velocity(j,:)) - call cpeval(tmp=temperature(j),spc=species(j,:),ke=var1,eng=var2) - q(j,5)=var2*density(j) - enddo - ! - else - print*,' !! temperature required for energy calculation !!' - stop ' !! error @ fvar2q' - endif - ! - endif !nondimen + print*,' !! temperature required for energy calculation !!' + stop ' !! error @ fvar2q' + endif + +#else + + if(nondimen) then + cotem=const1 + else + cotem=cv + endif + ! + if(present(temperature)) then + q(:,5)=density(:)*( temperature(:)*cotem + & + 0.5d0*(velocity(:,1)**2 + & + velocity(:,2)**2 + & + velocity(:,3)**2) ) + elseif(present(pressure)) then + q(:,5)=pressure(:)*const6+0.5d0*density(:)*( & + velocity(:,1)**2 + & + velocity(:,2)**2 + & + velocity(:,3)**2 ) + else + print*,' !! pressure or temperature required for energy calculation !!' + stop ' !! error @ fvar2q' + endif + +#endif ! if(num_species>0) then ! @@ -434,7 +448,7 @@ end subroutine fvar2q_1da ! subroutine fvar2q_3da(q,density,velocity,pressure,temperature,species,tke,omega) ! - use commvar, only: numq,ndims,num_species,const1,const6 + use commvar, only: numq,ndims,num_species,const1,const6,cv ! real(8),intent(in) :: density(:,:,:),velocity(:,:,:,:) real(8),intent(in),optional :: pressure(:,:,:),temperature(:,:,:), & @@ -444,56 +458,62 @@ subroutine fvar2q_3da(q,density,velocity,pressure,temperature,species,tke,omega) ! ! local data integer :: jspec,i,j,k - real(8) :: var1,var2 + real(8) :: var1,var2,cotem ! q(:,:,:,1)=density(:,:,:) q(:,:,:,2)=density(:,:,:)*velocity(:,:,:,1) q(:,:,:,3)=density(:,:,:)*velocity(:,:,:,2) q(:,:,:,4)=density(:,:,:)*velocity(:,:,:,3) ! - if(nondimen) then - ! - if(present(temperature)) then - q(:,:,:,5)=density(:,:,:)*( temperature(:,:,:)*const1 + & - 0.5d0*(velocity(:,:,:,1)**2 + & - velocity(:,:,:,2)**2 + & - velocity(:,:,:,3)**2) ) - elseif(present(pressure)) then - q(:,:,:,5)=pressure(:,:,:)*const6+0.5d0*density(:,:,:)*( & - velocity(:,:,:,1)**2 + & - velocity(:,:,:,2)**2 + & - velocity(:,:,:,3)**2 ) - else - print*,' !! pressure or temperature required !!' - stop ' !! error @ fvar2q' - endif +#ifdef COMB + + if(.not.present(species)) & + stop ' !! error @ fvar2q - species not set for energy calculation !!' + ! + if(present(temperature)) then + ! + do i=1,size(density,1) + do j=1,size(density,2) + do k=1,size(density,3) + ! + var1=0.5d0*sum(velocity(i,j,k,:)*velocity(i,j,k,:)) + call cpeval(tmp=temperature(i,j,k),spc=species(i,j,k,:), & + ke=var1,eng=var2) + q(i,j,k,5)=var2*density(i,j,k) + ! + enddo + enddo + enddo ! else - ! - if(.not.present(species)) & - stop ' !! error @ fvar2q - species not set for energy calculation !!' - ! - if(present(temperature)) then - ! - do i=1,size(density,1) - do j=1,size(density,2) - do k=1,size(density,3) - ! - var1=0.5d0*sum(velocity(i,j,k,:)*velocity(i,j,k,:)) - call cpeval(tmp=temperature(i,j,k),spc=species(i,j,k,:), & - ke=var1,eng=var2) - q(i,j,k,5)=var2*density(i,j,k) - ! - enddo - enddo - enddo - ! - else - print*,' !! temperature required for energy calculation !!' - stop ' !! error @ fvar2q' + print*,' !! temperature required for energy calculation !!' + stop ' !! error @ fvar2q' endif - ! - endif !nondimen + +#else + + if(nondimen) then + cotem=const1 + else + cotem=cv + endif + ! + if(present(temperature)) then + q(:,:,:,5)=density(:,:,:)*( temperature(:,:,:)*cotem + & + 0.5d0*(velocity(:,:,:,1)**2 + & + velocity(:,:,:,2)**2 + & + velocity(:,:,:,3)**2) ) + elseif(present(pressure)) then + q(:,:,:,5)=pressure(:,:,:)*const6+0.5d0*density(:,:,:)*( & + velocity(:,:,:,1)**2 + & + velocity(:,:,:,2)**2 + & + velocity(:,:,:,3)**2 ) + else + print*,' !! pressure or temperature required !!' + stop ' !! error @ fvar2q' + endif + +#endif ! if(num_species>0) then ! @@ -553,55 +573,56 @@ subroutine q2fvar_3da(q,density,velocity,pressure,temperature,species,tke,omega) ! endif ! - if(nondimen) then - ! - if(present(pressure) .or. present(temperature)) then - pressure(:,:,:) =( q(:,:,:,5)-0.5d0*density(:,:,:)*( & - velocity(:,:,:,1)**2+ & - velocity(:,:,:,2)**2+ & - velocity(:,:,:,3)**2) )/const6 - endif - ! - if(present(temperature)) then - dim(1)=size(q,1) - dim(2)=size(q,2) - dim(3)=size(q,3) - ! - temperature=thermal(pressure=pressure,density=density,dim=dim) - endif +#ifdef COMB + + if(.not.present(species)) & + stop ' !! error @ q2fvar - species not set for T calculation !!' + ! + if(present(pressure) .or. present(temperature)) then + ! + do i=1,size(q,1) + do j=1,size(q,2) + do k=1,size(q,3) + ! + var1=q(i,j,k,5)/density(i,j,k) & + -0.5d0*sum(velocity(i,j,k,:)*velocity(i,j,k,:)) + ! temperature(i,j,k)=tinf + call temperature_calc(tmp=temperature(i,j,k),den=density(i,j,k), & + spc=species(i,j,k,:),eint=var1) + enddo + enddo + enddo ! - else + endif + ! + if(present(pressure)) then + dim(1)=size(q,1) + dim(2)=size(q,2) + dim(3)=size(q,3) ! - if(.not.present(species)) & - stop ' !! error @ q2fvar - species not set for T calculation !!' - ! - if(present(pressure) .or. present(temperature)) then - ! - do i=1,size(q,1) - do j=1,size(q,2) - do k=1,size(q,3) - ! - var1=q(i,j,k,5)/density(i,j,k) & - -0.5d0*sum(velocity(i,j,k,:)*velocity(i,j,k,:)) - ! temperature(i,j,k)=tinf - call temperature_calc(tmp=temperature(i,j,k),den=density(i,j,k), & - spc=species(i,j,k,:),eint=var1) - enddo - enddo - enddo - ! - endif + pressure=thermal(temperature=temperature,density=density,species=species,dim=dim) + + endif + ! +#else + + if(present(pressure) .or. present(temperature)) then + pressure(:,:,:) =( q(:,:,:,5)-0.5d0*density(:,:,:)*( & + velocity(:,:,:,1)**2+ & + velocity(:,:,:,2)**2+ & + velocity(:,:,:,3)**2) )/const6 + endif + ! + if(present(temperature)) then + dim(1)=size(q,1) + dim(2)=size(q,2) + dim(3)=size(q,3) ! - if(present(pressure)) then - dim(1)=size(q,1) - dim(2)=size(q,2) - dim(3)=size(q,3) - ! - pressure=thermal(temperature=temperature,density=density, & - species=species,dim=dim) - endif - ! - endif !nondimen + temperature=thermal(pressure=pressure,density=density,dim=dim) + endif + ! + +#endif ! if(present(tke)) then tke(:,:,:)=q(:,:,:,5+num_species+1)/density @@ -644,45 +665,45 @@ subroutine q2fvar_1da(q,density,velocity,pressure,temperature,species,tke,omega) ! endif ! - if(nondimen) then +#ifdef COMB + + if(.not.present(species)) & + stop ' !! error @ q2fvar - species not set for T calculation !!' ! - if(present(pressure) .or. present(temperature)) then - pressure(:) =( q(:,5)-0.5d0*density(:)*( & - velocity(:,1)**2+ & - velocity(:,2)**2+ & - velocity(:,3)**2) )/const6 - endif + if(present(pressure) .or. present(temperature)) then ! - if(present(temperature)) then - dim=size(q,1) + do j=1,size(density) ! - temperature=thermal(pressure=pressure,density=density,dim=dim) - endif - ! - else + var1=q(j,5)/density(j)-0.5d0*sum(velocity(j,:)*velocity(j,:)) + call temperature_calc(tmp=temperature(j),den=density(j), & + spc=species(j,:),eint=var1) + enddo ! - if(.not.present(species)) & - stop ' !! error @ q2fvar - species not set for T calculation !!' - ! - if(present(pressure) .or. present(temperature)) then - ! - do j=1,size(density) - ! - var1=q(j,5)/density(j)-0.5d0*sum(velocity(j,:)*velocity(j,:)) - call temperature_calc(tmp=temperature(j),den=density(j), & - spc=species(j,:),eint=var1) - enddo - ! - endif + endif + ! + if(present(pressure)) then + dim=size(q,1) ! - if(present(pressure)) then - dim=size(q,1) - ! - pressure=thermal(temperature=temperature,density=density, & - species=species,dim=dim) - endif + pressure=thermal(temperature=temperature,density=density, & + species=species,dim=dim) + endif + +#else + + if(present(pressure) .or. present(temperature)) then + pressure(:) =( q(:,5)-0.5d0*density(:)*( & + velocity(:,1)**2+ & + velocity(:,2)**2+ & + velocity(:,3)**2) )/const6 + endif + ! + if(present(temperature)) then + dim=size(q,1) ! - endif !nondimen + temperature=thermal(pressure=pressure,density=density,dim=dim) + endif + +#endif ! if(present(tke)) then tke(:)=q(:,5+num_species+1)/density @@ -723,35 +744,33 @@ subroutine q2fvar_sca(q,density,velocity,pressure,temperature,species,tke,omega) ! endif ! - if(nondimen) then - ! - if(present(pressure) .or. present(temperature)) then - pressure =( q(5)-0.5d0*density*(velocity(1)**2+velocity(2)**2+ & - velocity(3)**2) )/const6 - endif - ! - if(present(temperature)) then - temperature=thermal(pressure=pressure,density=density) - endif - ! - else - ! - if(.not.present(species)) & - stop ' !! error @ q2fvar - species not set for T calculation !!' - ! - if(present(pressure) .or. present(temperature)) then - ! - var1=q(5)/density-0.5d0*sum(velocity(:)*velocity(:)) - call temperature_calc(tmp=temperature,den=density, & - spc=species(:),eint=var1) - endif +#ifdef COMB + if(.not.present(species)) & + stop ' !! error @ q2fvar - species not set for T calculation !!' ! - if(present(pressure)) then - pressure=thermal(temperature=temperature,density=density, & - species=species) - endif + if(present(pressure) .or. present(temperature)) then ! + var1=q(5)/density-0.5d0*sum(velocity(:)*velocity(:)) + call temperature_calc(tmp=temperature,den=density, & + spc=species(:),eint=var1) endif + ! + if(present(pressure)) then + pressure=thermal(temperature=temperature,density=density, & + species=species) + endif +#else + + if(present(pressure) .or. present(temperature)) then + pressure =( q(5)-0.5d0*density*(velocity(1)**2+velocity(2)**2+ & + velocity(3)**2) )/const6 + endif + ! + if(present(temperature)) then + temperature=thermal(pressure=pressure,density=density) + endif + +#endif ! if(present(tke)) then tke=q(5+num_species+1)/density @@ -780,15 +799,15 @@ real(8) function miucal(temper) ! tempconst=110.4d0/ref_t ! tempconst1=1.d0+tempconst ! - real(8) :: s,tmpers,tmper0,miu0,tnond + real(8) :: s,tmpers,tmper0,miu0,tnondim ! if(nondimen) then miucal=temper*sqrt(temper)*tempconst1/(temper+tempconst) else ! for air, kg/(m s) - tnond=temper/273.15d0 + tnondim=temper/273.15d0 ! - miucal=1.716d-5*tnond*sqrt(tnond)*(273.15d0+110.4d0)/(temper+110.4d0) + miucal=1.716d-5*tnondim*sqrt(tnondim)*(273.15d0+110.4d0)/(temper+110.4d0) ! ! tmpers=110.d0 ! tmper0=273.125d0 @@ -813,7 +832,7 @@ end function miucal !+-------------------------------------------------------------------+ real(8) function sos(tmp,spc) ! - use commvar, only: nondimen,mach,num_species + use commvar, only: nondimen,mach,num_species,gamma,rgas use thermchem, only: aceval ! real(8),intent(in) :: tmp @@ -825,11 +844,15 @@ real(8) function sos(tmp,spc) if(nondimen) then sos=sqrt(tmp)/mach else +#ifdef COMB if(.not. present(spc)) then stop ' !! error, species needed to calculate speed of sound @ function sos' else call aceval(tmp,spc,sos) endif +#else + sos=sqrt(gamma*rgas*tmp) +#endif endif ! return diff --git a/src/solver.F90 b/src/solver.F90 index e70679d..fccf7ad 100644 --- a/src/solver.F90 +++ b/src/solver.F90 @@ -27,7 +27,7 @@ module solver !+-------------------------------------------------------------------+ subroutine refcal ! - use commvar, only : numq,num_species,prandtl,gamma,rgas,ia,ja,ka, & + use commvar, only : numq,num_species,prandtl,gamma,rgas,cp,cv,ia,ja,ka, & uinf,vinf,winf,roinf,pinf,tinf,const1,const2, & const3,const4,const5,const6,const7,tempconst, & tempconst1,reynolds,ref_t,mach, & @@ -77,7 +77,6 @@ subroutine refcal endif ! prandtl=0.72d0 - rgas=287.1d0 ! if(nondimen) then ! @@ -104,6 +103,9 @@ subroutine refcal ! else ! + rgas=287.1d0 + cp =gamma/(gamma-1.d0)*rgas + cv = rgas/(gamma-1.d0) uinf=1.d0 vinf=0.d0 winf=0.d0 @@ -2329,7 +2331,8 @@ subroutine diffrsdcal6(timerept) use commvar, only : im,jm,km,numq,npdci,npdcj,npdck,difschm, & conschm,ndims,num_species,num_modequ, & reynolds,prandtl,const5,is,ie,js,je,ks,ke, & - turbmode,nondimen,schmidt,nstep,deltat,flowtype + turbmode,nondimen,schmidt,nstep,deltat, & + cp,flowtype use commarray, only : vel,tmp,spc,dvel,dtmp,dspc,dxi,x,jacob,qrhs, & rho,vor,omg,tke,miut,dtke,domg,res12 use commfunc, only : ddfc @@ -2409,6 +2412,7 @@ subroutine diffrsdcal6(timerept) else ! #ifdef COMB + call enthpy(tmp(i,j,k),hispec(:)) call convertxiyi(spc(i,j,k,:),xi(:),'Y2X') mw=sum(wmolar(:)*xi(:)) @@ -2424,6 +2428,8 @@ subroutine diffrsdcal6(timerept) spc=spc(i,j,k,:),rhodi=dispec(:,1)) ! print*,' ** miu=',miu,'cp=',cpe,'kamma=',kama,'rhodi=',dispec(:,1) end select +#else + miu=miucal(tmp(i,j,k)) #endif ! endif @@ -2484,7 +2490,11 @@ subroutine diffrsdcal6(timerept) if(nondimen) then hcc=(miu/prandtl)/const5 else +#ifdef COMB hcc=kama +#else + hcc=cp*miu/prandtl +#endif endif ! detk=0.d0 diff --git a/src/thermchem.F90 b/src/thermchem.F90 index 78a52cb..230a846 100644 --- a/src/thermchem.F90 +++ b/src/thermchem.F90 @@ -1311,10 +1311,11 @@ function rgcmix_3d(spc,dim) result(vout) real(8),intent(in) :: spc(:,:,:,:) real(8) :: vout(dim(1),dim(2),dim(3)) ! -#ifdef COMB !local integer :: i,j,k ! +#ifdef COMB + ! do i=1,dim(1) do j=1,dim(2) do k=1,dim(3) @@ -1322,6 +1323,7 @@ function rgcmix_3d(spc,dim) result(vout) enddo enddo enddo + ! #endif ! end function rgcmix_3d