diff --git a/bbb/bbb.v b/bbb/bbb.v index 34309f35..cd83efa7 100755 --- a/bbb/bbb.v +++ b/bbb/bbb.v @@ -91,6 +91,7 @@ tibg real [eV] /1.e-20/ +input #backgrd ion eng sor to limit te~tebg iteb integer /2/ +input #exponent of (tebg*ev/te)**iteb for bkg sor temin real [eV] /0.03/ +input #min value of te allow; if less, reset to temin2 real [eV] /0.03/ +input #soft floor with te=sqrt[te**2+(temin2*ev)**2] +tgmin real [eV] /0.03/ # min value of tg allowed pwrbkg_c real [W/m**3] /1.e3/ +input #const background factor in pwrebkg express pwribkg_c real [W/m**3] /1.e3/ +input #const background factor in pwribkg express cfwjdotelim real /1./ +input #factor scaling reduction of wjdote if te_i to cne_sgvi cne_sgvi real [1/m**3] /1.e18/ +input #ne for _i if ifxnsgi=1 ctsor real /1./ +input #Coef for eng. src. in Ti eq. 0.5*mi*up**2*psor @@ -324,6 +343,7 @@ vboost real /1./ +input #previously scaled eqp; no longer in use cvgp real /1./ +input #Coef for v.Grad(p) ion/elec eng. terms cvgpg real /1./ +input #Coef for v.Grad(pg) gas eng. terms cfvgpx(1:nispmx) real /nispmx*1./ +input #Coefs for x components of v.grad(p) in ti-eq +cftiexclg real /1./ +input #Coef for including atom gas in Ti eq. cfvgpy(1:nispmx) real /nispmx*1./ +input #Coefs for y components of v.grad(p) in ti-eq cfvgpgx(1:ngspmx) real /ngspmx*0./ +input # Coefs for v*grad(pg) term (poloidal component) for Tg if isupgon=0. cfvgpgy(1:ngspmx) real /ngspmx*0./ +input # Coefs for v*grad(pg) term (poloidal component) for Tg if isupgon=0. @@ -469,6 +489,7 @@ istgcore(ngspmx) integer /ngspmx*1/ +input #switch for neutral-density core B.C. #=0, set tg(ixcore,0,igsp)=ti(ixcore,0) #=1, set fixed temp tgcore(igsp) #if > 1, set zero grad; tg(,0,)=tg(,1,) +cftgticore(ngspmx) real /ngspmx*1./ +input #.. curcore(1:nispmx) real [A] /0.,30*0./ +input #value of current from core if isnicore=0 lzcore(1:nispmx) real [kg/ms] /nispmx*0./ +input #tor. ang. mom dens core bdry; phi eqn lzflux(1:nispmx) real [kg/s**2]/nispmx*0./ +input #tor. ang. mom dens flux core bdry; up eqn @@ -496,6 +517,8 @@ tepltl real [eV] /2./ +input #left plate Te B.C. if ibctepl=0 tipltl real [eV] /2./ +input #left plate Ti B.C. if ibctipl=0 tepltr real [eV] /2./ +input #right plate Te B.C. if ibcteplr=0 tipltr real [eV] /2./ +input #right plate Ti B.C. if ibctiplr=0 +cftgtipltl(ngspmx) real /ngspmx*1./ +input #left plate Tg B.C. if ibctgpl=? +cftgtipltr(ngspmx) real /ngspmx*1./ +input #right plate Tg B.C. if ibctgpr=? tbmin real [eV] /.1/ +input #min. wall & pf temp for extrap. b.c.(isextrt..) nbmin real [m**-3] /1.e17/ +input #min. wall & pf den for extrap. b.c.(isextrn..) ngbmin real [m**-3] /1.e10/ +input #min. core gas den for extrap. b.c.(isextrngc) @@ -525,6 +548,8 @@ istgpfc(ngspmx) integer/ngspmx*0/ +input #switch for PF BC on Tg(,0,igsp) # =2, set Tg scale length to lytg(1, # =3, eng flux = 2Tg*Maxw-flux # >2, report error in input +cftgtiwc(ngspmx) real /ngspmx*1./ +input #wall Tg B.C. if istgwc=? +cftgtipfc(ngspmx) real /ngspmx*1./ +input #pfc Tg B.C. if istgpfc=? tewalli(0:nx+1) _real [eV] +input #/(nx+2)*0./ #inner wall Te for istepfc=1.; = tedge if not set tiwalli(0:nx+1) _real [eV] +input #/(nx+2)*0./ @@ -871,12 +896,14 @@ recycmrb_use(0:ny+1,ngspmx,nxptmx) _real +maybeinput #outer plt mom-recycl coeff recycmlb(0:ny+1,ngspmx,nxptmx) _real +maybeinput #total inner plt mom recycling coeff recycmrb(0:ny+1,ngspmx,nxptmx) _real +maybeinput #total outer plt mom recycling coeff recyce real /0./ +input #energy recycling/Rp for inertial gas +recycwe real /0./ +input #energy recycling/Rp for inertial gas on walls and prfs recycl real /1./ +input #recycling coef. at a limiter (ix_lim) recycml real /0.1/ +input #momentum recycling/Rp for gas at limtr recycc(ngspmx) real /6*1./ +input #core recycling coeff. if isnicore=3 albedoc(ngspmx) real /6*1./ +input #core neut albedo for isngcore=0 albedolb(ngspmx,nxptmx) _real /1./ +input #albedo at inner plate if ndatlb=0 albedorb(ngspmx,nxptmx) _real /1./ +input #albedo at outer plate if ndatrb=0 +cfalbedo real /2./ +input #coef. for Tg Eq. due to albedo ndatlb(ngspmx,nxptmx) _integer /0/ +maybeinput #number of recycp data pts on inner plt ndatrb(ngspmx,nxptmx) _integer /0/ +maybeinput #number of recycp data pts on outer plt ydatlb(ngspmx,50,nxptmx) _real [m] /0./ +maybeinput #inner data pt location from sep. @@ -1701,6 +1728,7 @@ volpsor(0:nx+1,0:ny+1,1:nisp) _real [1/s]+maybeinput #current src into ions in c volmsor(0:nx+1,0:ny+1,1:nisp) _real [kg m/s**2] +maybeinput #up mom src in cell ix,iy voljcsor(0:nx+1,0:ny+1) _real [A] +maybeinput #uniform core-region curr sor. in ix,iy volpsorg(0:nx+1,0:ny+1,1:ngsp) _real [1/s]+maybeinput #curr source for gas in cell ix,iy +pwrsorg(0:nx+1,0:ny+1,1:ngsp) _real [1/s]+maybeinput #power source for gas in cell ix,iy; only for manufactured solution verification pondpot(0:nx+1,0:ny+1) _real [V] /0./ +maybeinput #elec ponderomotive potential psgov_use(0:nx+1,0:ny+1,1:ngsp) _real [1/m**3 s]+maybeinput #user-specified gas source jcvsor real [A] /0./ +maybeinput #total core-region current for voljcsor @@ -2076,6 +2104,7 @@ cfvli(nisp) _real #/nisp*0./+input #scal fac for individ ion rate nuvl l_parloss real [m] /1.e20/ +input #parall length for nuvl loss rate eqp(0:nx+1,0:ny+1) _real [1/m**3]#Te,i equipart. fact; needs *(Te-Ti)*vol eqpg(0:nx+1,0:ny+1,ngsp) _real [1/m**3]#Tg,i equipart. fact; needs *(Tg-Ti)*vol + #..: modified to incorporate the separation of Tg and Ti. engcoolm(0:nx+1,0:ny+1) _real [J/s] #cool rate ion/atoms by mols if ishymol=1 eeli(0:nx+1,0:ny+1) _real [J] #electron energy loss per ionization pradhyd(0:nx+1,0:ny+1) _real [W/m**3] /0./#power radiated by hydrogen diff --git a/bbb/boundary.m b/bbb/boundary.m index a83b2be6..3585e70b 100755 --- a/bbb/boundary.m +++ b/bbb/boundary.m @@ -130,7 +130,7 @@ c_mpi Use(MpiVars) #module defined in com/mpivarsmod.F.in if (isupgon(1) .eq. 1 .and. zi(ifld) .eq. 0.0) then if (isixcore(ix)==1) then # ix is part of the core boundary: if (isngcore(1) .eq. 0) then - t0 = max(tg(ix,0,1),temin*ev) + t0 = max(tg(ix,0,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(ifld)) ) nharmave = 2.*(ni(ix,0,ifld)*ni(ix,1,ifld)) / . (ni(ix,0,ifld)+ni(ix,1,ifld)) @@ -162,7 +162,7 @@ call xerrab("*** isngcore=2 option unvailable ***") endif else # ix is not part of the core boundary if (iscpli(ix) .eq. 1) call wsmodi(1) - t0 = max(tg(ix,0,1),temin*ev) + t0 = max(tg(ix,0,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(1)) ) fng_chem = 0. do ii = 1, ngsp #chem sputt of hydrogen - strange = 0 @@ -614,12 +614,12 @@ ccc yldot(iv2) = -nurlxi*(feiy(ix,0) - . 0.5*(ti(ix,1) + ti(ix,0))/ . (gyf(ix,0)*lytix(1,ix)) )/(temp0*ev) elseif (istipfc .eq. 4) then # sheath-like condition - t0 = max(tg(ix,0,1), temin*ev) + t0 = max(tg(ix,0,1), tgmin*ev) fngyw=0.25*sqrt(8*t0/(pi*mg(1)))*ng(ix,1,1)*sy(ix,0) yldot(iv2) = -nurlxi*( feiy(ix,0) - bceiw*fniy(ix,0,1)* . ti(ix,0) - - . bcenw*fniy(ix,0,iigsp)*tg(ix,0,1) + - . cgengw*2.*tg(ix,0,iigsp)*fngyw ) / + . cftiexclg*( bcenw*fniy(ix,0,iigsp)*tg(ix,0,1) - + . cgengw*2.*tg(ix,0,iigsp)*fngyw ) ) / . (vpnorm*ennorm*sy(ix,0)) endif endif # end if-test for core and p.f. boundaries @@ -636,9 +636,9 @@ ccc yldot(iv2) = -nurlxi*(feiy(ix,0) - if (isngonxy(ix,0,igsp) .eq. 1) then # ends just before 275 continue statem. if (iscpli(ix) .eq. 1) call wsmodi(igsp) iv = idxg(ix,0,igsp) - t0 = max(cdifg(igsp)*tg(ix,0,igsp),temin*ev) + t0 = max(cdifg(igsp)*tg(ix,0,igsp),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) - t1 = engbsr * max(tg(ix,0,1),temin*ev) + t1 = engbsr * max(tg(ix,0,1),tgmin*ev) c ... prepare impurity ion flux for possible recycling zflux = 0. if (igsp .gt. nhgsp) then @@ -683,7 +683,7 @@ ccc yldot(iv2) = -nurlxi*(feiy(ix,0) - ### wall BCs. calc diff neut sputt fluxes fng_chem = 0. do igsp2 = 1, 1+ishymol #include hydrogen neut sputt - t0 = max(tg(ix,0,igsp2),temin*ev) + t0 = max(tg(ix,0,igsp2),tgmin*ev) flx_incid=ng(ix,1,igsp2)*0.25*sqrt(8*t0/(pi*mg(igsp2))) if (isch_sput(igsp).ne.0) then call sputchem (isch_sput(igsp), t0/ev, tvwalli(ix), @@ -772,10 +772,24 @@ call sputchem (isch_sput(igsp), eng_sput, tvwalli(ix), iv = idxtg(ix,0,igsp) if (isixcore(ix)==1) then if (istgcore(igsp) == 0) then #set to ti - yldot(iv)=nurlxg*(ti(ix,0)-tg(ix,0,igsp))/(temp0*ev) + yldot(iv)=nurlxg*(ti(ix,0)*cftgticore(igsp) + . -tg(ix,0,igsp))/(temp0*ev) elseif (istgcore(igsp) == 1) then # set to tgcore yldot(iv)=nurlxg*(tgcore(igsp)*ev-tg(ix,0,igsp))/ . (temp0*ev) + elseif (istgcore(igsp) == 2) then + #if (isupgon(igsp)==1) then + t0 = max(tg(ix,0,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + nharmave = 2.*(ng(ix,0,igsp)*ng(ix,1,igsp)) / + . (ng(ix,0,igsp)+ng(ix,1,igsp)) + fng_alb = (1-albedoc(igsp))* + . nharmave*vyn*sy(ix,0) + yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0) + . /(vpnorm*ennorm*sy(ix,0)) + #else + # write(*,*) 'Istgcore',igsp,'==1 is not available' + #endif else # all others, set to zero y-gradient yldot(iv)=nurlxg*(tg(ix,1,igsp)-tg(ix,0,igsp))/ . (temp0*ev) @@ -794,12 +808,36 @@ call sputchem (isch_sput(igsp), eng_sput, tvwalli(ix), . 0.5*(tg(ix,1,igsp) + tg(ix,0,igsp))/ . (gyf(ix,0)*lytg(1,igsp)) )/(temp0*ev) elseif (istgpfc(igsp) == 3) #Maxwell thermal flux to wall - t0 = max(cdifg(igsp)*tg(ix,1,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ix,1,igsp), tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = -nurlxg*( fegy(ix,0,igsp) + 2*cgengmw* . ng(ix,1,igsp)*vyn*t0*sy(ix,0) )/ . (sy(ix,0)*vpnorm*ennorm) - elseif (istgpfc(igsp) > 3) + elseif (istgpfc(igsp) == 4) + t0 = max(tg(ix,0,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + nharmave = 2.*(ng(ix,0,igsp)*ng(ix,1,igsp)) / + . (ng(ix,0,igsp)+ng(ix,1,igsp)) + fng_alb = (1-albedoi(ix,1))* + . nharmave*vyn*sy(ix,0) + fng_chem = 0. #..double check + yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 + . - 2.*fng_chem*t0) + . /(vpnorm*ennorm*sy(ix,0)) + fniy_recy = 0. + if (matwalli(ix) .gt. 0) then + if (recycwit(ix,igsp,1) .gt. 0) then + fniy_recy = recycwit(ix,igsp,1)*fac2sp*fniy(ix,0,1) + if (isrefluxclip==1) fniy_recy=min(fniy_recy,0.) + yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 + . - 2.*fng_chem*t0 + . + fniy_recy*(1.-cfdiss) + . *cfalbedo*recycwe + . *ti(ix,0) ) + . /(vpnorm*ennorm*sy(ix,0)) + endif + endif + elseif (istgpfc(igsp) > 4) call xerrab("***Input error: invalid istgpfc ***") endif @@ -1097,7 +1135,7 @@ call xerrab("**INPUT ERROR: only iphibcc=1,2,3 available") c --- typically hydrogen only, so DIVIMP chem sputt not used here if (isupgon(1) .eq. 1 .and. zi(ifld) .eq. 0.0) then if (iscplo(ix) .eq. 1) call wsmodo(1) - t0 = max(tg(ix,ny+1,1),temin*ev) + t0 = max(tg(ix,ny+1,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(ifld)) ) fng_chem = 0. do ii = 1, ngsp @@ -1307,12 +1345,12 @@ c Force corner cells (0,ny+1) and (nx+1,ny+1) to relax to an average . 0.5*(ti(ix,ny) + ti(ix,ny+1))/ . (gyf(ix,ny)*lytix(2,ix)) )/(temp0*ev) elseif (istiwc .eq. 4) then #heat flux ~bcei*ti*elec_flux - t0 = max(tg(ix,ny+1,1), temin*ev) + t0 = max(tg(ix,ny+1,1), tgmin*ev) fngyw=0.25*sqrt(8*t0/(pi*mg(1)))*ng(ix,ny,1)*sy(ix,ny) yldot(iv2) = nurlxi*( feiy(ix,ny) - bceiw*ne(ix,ny)* . vey(ix,ny)*sy(ix,ny)*ti(ix,ny+1) - - . bcenw*fniy(ix,ny,iigsp)*tg(ix,ny+1,1) - - . cgengw*2.*tg(ix,ny+1,1)*fngyw ) / + . cftiexclg*( bcenw*fniy(ix,ny,iigsp)*tg(ix,ny+1,1) + + . cgengw*2.*tg(ix,ny+1,1)*fngyw ) ) / . (vpnorm*ennorm*sy(ix,ny)) endif endif @@ -1328,7 +1366,7 @@ ccc Now do the diffusive neutral density (ng) and Tg equations if(isngonxy(ix,ny+1,igsp) .eq. 1) then #skip do-loop if isngon=0 if (iscplo(ix) .eq. 1) call wsmodo(igsp) iv = idxg(ix,ny+1,igsp) - t0 = max(cdifg(igsp)*tg(ix,ny+1,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ix,ny+1,igsp), tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) if (igsp .gt. nhgsp) then zflux = 0. @@ -1339,7 +1377,7 @@ ccc Now do the diffusive neutral density (ng) and Tg equations c ... prepare chemical sputtering info fng_chem = 0. do igsp2 = 1, 1+ishymol #include hydrogen neut sputt - t0 = max(tg(ix,ny+1,igsp2), temin*ev) + t0 = max(tg(ix,ny+1,igsp2), tgmin*ev) flx_incid=ng(ix,ny+1,igsp2)*0.25*sqrt(8*t0/(pi*mg(igsp2))) if (isch_sput(igsp).ne.0) then call sputchem (isch_sput(igsp), t0/ev, tvwallo(ix), @@ -1434,12 +1472,35 @@ call sputchem (isch_sput(igsp), eng_sput, tvwallo(ix), . 0.5*(tg(ix,ny,igsp) + tg(ix,ny+1,igsp))/ . (gyf(ix,ny)*lytg(2,igsp)) )/(temp0*ev) elseif (istgwc(igsp) == 3) #Maxwell thermal flux to wall - t0 = max(cdifg(igsp)*tg(ix,ny,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ix,ny,igsp), tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = nurlxg*( fegy(ix,ny,igsp) - 2*cgengmw* . ng(ix,ny,igsp)*vyn*t0*sy(ix,ny) )/ . (sy(ix,ny)*vpnorm*ennorm) - elseif (istgwc(igsp) > 3) + elseif (istgwc(igsp) == 4) + t0 = max(tg(ix,ny+1,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + nharmave = 2.*(ng(ix,ny,igsp)*ng(ix,ny+1,igsp)) / + , (ng(ix,ny,igsp)+ng(ix,ny+1,igsp)) + fng_alb = (1-albedoo(ix,1))*nharmave*vyn*sy(ix,ny) + fng_chem = 0. + yldot(iv) = nurlxg*(fegy(ix,ny,igsp) - cfalbedo*fng_alb*t0 + . + 2.*fng_chem*t0) + . /(vpnorm*ennorm*sy(ix,ny)) + fniy_recy = 0. + if (matwallo(ix) .gt. 0) then + if (recycwot(ix,igsp) .gt. 0.) then + fniy_recy = recycwot(ix,igsp)*fac2sp*fniy(ix,ny,1) + if (isrefluxclip==1) fniy_recy=max(fniy_recy,0.) + yldot(iv)=nurlxg*(fegy(ix,ny,igsp) - cfalbedo*fng_alb*t0 + . + 2.*fng_chem*t0 + . + fniy_recy*(1.-cfdiss) + . *cfalbedo*recycwe + . *ti(ix,ny) ) + . /(vpnorm*ennorm*sy(ix,ny)) + endif + endif + elseif (istgwc(igsp) > 4) call xerrab("***Input error: invalid istgwc ***") endif endif @@ -1502,6 +1563,12 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 . ( ng(ixlb(jx)+1,ny+1,igsp) - ng(ixlb(jx),ny+1,igsp) ) . / n0g(igsp) endif + if(istgonxy(ixlb(jx),ny+1,igsp)==1) then + yldot(idxtg(ixlb(jx),ny+1,igsp)) = nurlxg* + . ( 0.5*(tg(ixlb(jx)+1,ny+1,igsp)+tg(ixlb(jx),ny,igsp)) + . - tg(ixlb(jx),ny+1,igsp) ) + . / (temp0*ev) + endif enddo enddo endif @@ -1543,6 +1610,12 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 . ( ng(ixrb(jx),ny+1,igsp) - ng(ixrb(jx)+1,ny+1,igsp) ) . / n0g(igsp) endif + if(istgonxy(ixrb(jx)+1,ny+1,igsp)==1) then + yldot(idxtg(ixrb(jx)+1,ny+1,igsp)) = nurlxg* + . ( 0.5*(tg(ixrb(jx),ny+1,igsp)+tg(ixrb(jx)+1,ny,igsp)) + . - tg(ixrb(jx)+1,ny+1,igsp) ) + . / (temp0*ev) + endif enddo enddo endif # end of if-test on xcnearrb @@ -1641,7 +1714,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 if(isfixlb(1).eq.2) yldot(iv) = nurlxg * . (ng(1,iy,igsp) - ng(0,iy,igsp))/n0g(igsp) if(isfixlb(1).eq.2 .and. yylb(iy,1).gt.rlimiter) then - t1 = engbsr * max(tg(1,iy,igsp),temin*ev) + t1 = engbsr * max(tg(1,iy,igsp),tgmin*ev) vxn = 0.25 * sqrt( 8*t1/(pi*mg(igsp)) ) flux_inc = fac2sp*fnix(0,iy,1) if (ishymol.eq.1 .and. igsp.eq.2) then @@ -1722,7 +1795,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 iv1 = idxn(ixt,iy,ifld) if (isupgon(1)==1 .and. zi(ifld)==0.0) then ## neutrals if (recylb(iy,1,jx) .gt. 0.) then # recycling - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mi(ifld)) ) yldot(iv1) = -nurlxg * . (fnix(ixt,iy,ifld) + recylb(iy,1,jx)*fnix(ixt,iy,1) - @@ -1731,7 +1804,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 . fngxslb(iy,1,jx) ) / (vpnorm*n0(ifld)*sx(ixt,iy)) elseif (recylb(iy,1,jx) <= 0. .and. . recylb(iy,1,jx) >= -1.) then # recylb is albedo - t0 = max(tg(ixt,iy,1),temin*ev) + t0 = max(tg(ixt,iy,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(1)) ) yldot(iv1) = -nurlxg * ( fnix(ixt,iy,ifld) + . (1+recylb(iy,1,jx))*ni(ixt,iy,ifld)*vyn*sx(ixt,iy) )/ @@ -1793,7 +1866,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 yldot(iv2) = nurlxu*(up(ixt1,iy,ifld) - . up(ixt,iy,ifld))/vpnorm else # neutral thermal flux to wall if recycm < -10.1 - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vxn = cgmompl*0.25*sqrt( 8*t0/(pi*mi(ifld)) ) vparn = up(ixt,iy,ifld) yldot(iv2) = -nurlxu*( fmix(ixt1,iy,ifld) + vparn*vxn* @@ -1827,9 +1900,15 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0 endif # end if-test on isupss endif # end if-test on isupgon endif # end if-test on isuponxy - kfeix = kfeix - cfvcsx(ifld)*0.5*sx(ixt,iy) + if (zi(ifld)==0.0) then + kfeix = kfeix - cftiexclg*cfvcsx(ifld)*0.5*sx(ixt,iy) + . *visx(ixt1,iy,ifld)*gx(ixt1,iy) + . *( up(ixt1,iy,ifld)**2 - up(ixt,iy,ifld)**2 ) + else + kfeix = kfeix - cfvcsx(ifld)*0.5*sx(ixt,iy) . *visx(ixt1,iy,ifld)*gx(ixt1,iy) . *( up(ixt1,iy,ifld)**2 - up(ixt,iy,ifld)**2 ) + endif enddo # end do-loop on ifld c Next, the potential equation -- @@ -1897,7 +1976,7 @@ cc fqp(ixt,iy)=-( abs((fqp(ixt,iy)*fqpsatlb(iy,jx)))**exjbdry/ # temperatures (ETM 7 Jan 2014) if(t0 < temin) f_cgpld = 0. if(t0 > 0.3) f_cgpld = 1. - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vxn = f_cgpld * 0.25 * sqrt(8*t0/(pi*mg(1))) c Do the electron temp Eqn ----------------------------------- @@ -1955,20 +2034,21 @@ cc if (recylb(iy,1,jx) .gt. 0.) then cc . recyce*upi(ixt,iy,1)**2)/(recylb(iy,1,jx)*ti(ixt,iy)) cc endif cc if (recyce .le. 0) bcen = 0. # gets back to old case - yldot(iv2) = -nurlxi*(totfeixl(iy,jx) - . -totfnix*ti(ixt,iy)*bcil(iy,jx) - . -cfneut*fnix(ixt,iy,iigsp)*tg(ixt,iy,1)*bcen + + yldot(iv2) = -nurlxi*(totfeixl(iy,jx) + . -totfnix*ti(ixt,iy)*bcil(iy,jx)+ + . cftiexclg*( -cfneut*fnix(ixt,iy,iigsp)*tg(ixt,iy,1)*bcen . +(cgengpl*2.*tg(ixt,iy,1) - cgpld*eion*ev)* . ng(ixt1,iy,1)*vxn*sx(ixt,iy) . -cmneut*fnix(ixt,iy,1)*recycp(1)* . cmntgpl*(ti(ixt,iy)-eidisspl*ev) - . )/(vpnorm*ennorm*sx(ixt,iy)) + . ) )/(vpnorm*ennorm*sx(ixt,iy)) else - yldot(iv2) = -nurlxi*(totfeixl(iy,jx) - . -totfnix*bcil(iy,jx)*ti(ixt,iy) - . -cmneut*fnix(ixt,iy,1)*recycp(1)* + yldot(iv2) = -nurlxi*(totfeixl(iy,jx) + . -totfnix*bcil(iy,jx)*ti(ixt,iy)+ + . cftiexclg*( -cmneut*fnix(ixt,iy,1)*recycp(1)* . cmntgpl*(ti(ixt,iy)-eidisspl*ev) - . )/(vpnorm*ennorm*sx(ixt,iy)) + . ) )/(vpnorm*ennorm*sx(ixt,iy)) endif #end loop for isupgon(1) elseif (ibctipl .eq. 0) then @@ -1993,7 +2073,7 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case flux_inc = 0.5*( fnix(ixt,iy,1) + fngx(ixt,iy,1) ) endif endif - t0 = max(tg(ixt1,iy,igsp), temin*ev) + t0 = max(tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = -nurlxg * ( fngx(ixt,iy,igsp) - . fngxlb_use(iy,igsp,jx) - @@ -2002,7 +2082,7 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case . / (vpnorm*n0g(igsp)*sx(ixt,iy)) elseif (recylb(iy,igsp,jx) <= 0. .and. . recylb(iy,igsp,jx) >= -1.) then # recylb is albedo - t0 = max(tg(ixt,iy,igsp), temin*ev) + t0 = max(tg(ixt,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = -nurlxg*( fngx(ixt,iy,igsp) + . (1+recylb(iy,igsp,jx))*ng(ixt,iy,igsp)*vxn*sx(ixt,iy) ) @@ -2065,7 +2145,7 @@ call sputchem (isch_sput(igsp),eng_sput,tvplatlb(iy,jx), endif if (isph_sput(igsp) .eq. 3) then #add chem sput from h neuts do igsp2 = 1, 1+ishymol #hydrogen neut fluxes only - t0p = max(tg(ixt1,iy,igsp2),temin*ev) + t0p = max(tg(ixt1,iy,igsp2),tgmin*ev) flx_incid = ng(ixt,iy,igsp2)*.25* . sqrt(8*t0p/(pi*mg(igsp2))) call sputchem (isch_sput(igsp),t0p/ev,tvplatlb(iy,jx), @@ -2077,7 +2157,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatlb(iy,jx), endif if (sputtlb(iy,igsp,jx) .ge. 0. .or. . abs(sputflxlb(iy,igsp,jx)).gt. 0.) then - t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) zflux = - sputtlb(iy,igsp,jx) * hflux - . sputflxlb(iy,igsp,jx) - @@ -2087,7 +2167,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatlb(iy,jx), yldot(iv) = -nurlxg * (fngx(ixt,iy,igsp) - zflux) / . (n0(igsp) * vpnorm * sx(ixt,iy)) elseif (sputtlb(iy,igsp,jx).ge.-9.9) then # neg. sputtlb ==> albedo - t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = -nurlxg*( fngx(ixt,iy,igsp) - . (1+sputtlb(iy,igsp,jx))*ng(ixt,iy,igsp)*vxn*sx(ixt,iy) ) @@ -2115,11 +2195,43 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatlb(iy,jx), elseif (istglb(igsp) == 2) #placeholder for gradient BC call xerrab("**INPUT ERROR: istglb=2 grad opt not implemented") elseif (istglb(igsp) == 3) #Maxwell thermal flux to wall - t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = -nurlxg*( fegx(ixt,iy,igsp) + 2*cgengmpl* . ng(ixt1,iy,igsp)*vxn*t0*sx(ixt,iy) )/ . (sx(ixt,iy)*vpnorm*ennorm) + elseif (istglb(igsp) == 4) + if (isupgon(igsp)==1) then + if (recylb(iy,igsp,jx) .gt. 0.) then + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) + fng_alb=(1-alblb(iy,igsp,jx))*ng(ixt1,iy,igsp) + . *vxn*sx(ixt,iy) + yldot(iv) = -nurlxg*( fegx(ixt,iy,igsp) + . +cfalbedo*fng_alb*t0 + . +recylb(iy,igsp,jx)*(1.-cfdiss) + . *fnix(ixt,iy,1) + . *recyce*cfalbedo + . *( kappal(iy,jx)*zi(1)*te(ixt,iy) + . +ti(ixt,iy) ) )/ + . (vpnorm*ennorm*sx(ixt,iy)) + elseif (recylb(iy,igsp,jx) <= 0. .and. + . recylb(iy,igsp,jx) >= -1.) then # recylb is albedo + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + fng_alb = (1+recylb(iy,igsp,jx))*ng(ixt1,iy,igsp) + . *vyn*sx(ixt,iy) + yldot(iv) = -nurlxg*( fegx(ixt,iy,igsp) + . +cfalbedo*fng_alb*t0 )/ + . (vpnorm*ennorm*sx(ixt,iy)) + elseif (recylb(iy,igsp,jx) < -1.) then #..half Maxwellian + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + yldot(iv) = -nurlxg*( fegx(ixt,iy,igsp) + . +cfalbedo*fnix(ixt,iy,iigsp)*t0 )/ + . (vpnorm*ennorm*sx(ixt,iy)) + endif + endif else call xerrab("**INPUT ERROR: istglb set to unknown option") endif @@ -2258,7 +2370,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0 if(isfixrb(1).eq.2) yldot(iv) = nurlxg * . (ng(nx,iy,igsp) - ng(nx+1,iy,igsp))/n0g(igsp) if(isfixrb(1).eq.2 .and. yyrb(iy,1).gt.rlimiter) then - t1 = engbsr * max(tg(nx,iy,1),temin*ev) + t1 = engbsr * max(tg(nx,iy,1),tgmin*ev) vxn = 0.25 * sqrt( 8*t1/(pi*mg(igsp)) ) flux_inc = fac2sp*fnix(nx,iy,1) if (ishymol.eq.1 .and. igsp.eq.2) then @@ -2339,7 +2451,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0 iv1 = idxn(ixt,iy,ifld) if (isupgon(1)==1 .and. zi(ifld)==0.0) then ## neutrals if (recyrb(iy,1,jx) .gt. 0.) then # recycling - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mi(ifld)) ) yldot(iv1) = nurlxg * . (fnix(ixt1,iy,ifld) + recyrb(iy,1,jx)*fnix(ixt1,iy,1) + @@ -2348,7 +2460,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0 . - fngxsrb(iy,1,jx) ) / (vpnorm*n0(ifld)*sx(ixt1,iy)) elseif (recyrb(iy,1,jx) <= 0. .and. . recyrb(iy,1,jx) >= -1.) then # recyrb is albedo - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(1)) ) yldot(iv1) = nurlxg * ( fnix(ixt1,iy,ifld) - . (1+recyrb(iy,1,jx))*ni(ixt,iy,ifld)*vyn*sx(ixt1,iy) )/ @@ -2414,7 +2526,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0 yldot(iv2) = nurlxu*(up(ixt2,iy,ifld) - . up(ixt1,iy,ifld))/vpnorm else #neutral thermal flux to wall if recycm < -10.1 - t0 = max(tg(ixt,iy,1),temin*ev) + t0 = max(tg(ixt,iy,1),tgmin*ev) vxn = cgmompl*0.25*sqrt( 8*t0/(pi*mi(ifld))) c... if up > 0, leave unchanged; if up<0, big reduction cc vparn = up(ixt,iy,ifld)*( 1./(1. + @@ -2456,10 +2568,15 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0 yldot(iv) =nurlxu*(up(ixt1,iy,ifld)-up(ixt,iy,ifld))/vpnorm endif # end if-test on isupgon endif # end if-test on isupon - - kfeix = kfeix - cfvcsx(ifld)*0.5*sx(ixt1,iy) + if (zi(ifld)==0.0) then + kfeix = kfeix - cftiexclg*cfvcsx(ifld)*0.5*sx(ixt1,iy) . *visx(ixt1,iy,ifld)*gx(ixt1,iy) . *( up(ixt1,iy,ifld)**2 - up(ixt2,iy,ifld)**2 ) + else + kfeix = kfeix - cfvcsx(ifld)*0.5*sx(ixt1,iy) + . *visx(ixt1,iy,ifld)*gx(ixt1,iy) + . *( up(ixt1,iy,ifld)**2 - up(ixt2,iy,ifld)**2 ) + endif enddo # end do-loop on ifld c Next, the potential equation -- @@ -2527,7 +2644,7 @@ cccTDR if ((isnewpot==1) .and. ((iy==1) .or. (iy==ny))) then # temperatures (ETM 7 Jan 2014) if(t0 < temin) f_cgpld = 0. if(t0 > 0.3) f_cgpld = 1. - t0 = max(tg(ixt1,iy,1),temin*ev) + t0 = max(tg(ixt1,iy,1),tgmin*ev) vxn = f_cgpld * 0.25 * sqrt(8*t0/(pi*mg(1))) c Do the electron temp Eqn ----------------------------------- @@ -2592,20 +2709,21 @@ cc if (recyrb(iy,1,jx) .gt. 0.) then cc bcen = 0. cc endif cc if (recyce .le. 0) bcen = 0. # gets back to old case - yldot(iv2) = nurlxi*(totfeixr(iy,jx) - . -totfnix*bcir(iy,jx)*ti(ixt,iy) - . -cfneut*fnix(ixt1,iy,iigsp)*bcen*tg(ixt,iy,1) + + yldot(iv2) = nurlxi*(totfeixr(iy,jx) + . -totfnix*bcir(iy,jx)*ti(ixt,iy)+ + . cftiexclg*( -cfneut*fnix(ixt1,iy,iigsp)*bcen*tg(ixt,iy,1) . -(cgengpl*2.*tg(ixt,iy,1) - cgpld*eion*ev)* . ng(ixt1,iy,1)*vxn*sx(ixt1,iy) . -cmneut*fnix(ixt1,iy,1)*recycp(1)* . cmntgpr*(ti(ixt,iy)-eidisspr*ev) - . )/(vpnorm*ennorm*sx(ixt1,iy)) + . ) )/(vpnorm*ennorm*sx(ixt1,iy)) else - yldot(iv2) = nurlxi*(totfeixr(iy,jx) - . -totfnix*bcir(iy,jx)*ti(ixt,iy) - . -cmneut*fnix(ixt1,iy,1)*recycp(1)* + yldot(iv2) = nurlxi*(totfeixr(iy,jx) + . -totfnix*bcir(iy,jx)*ti(ixt,iy)+ + . cftiexclg*( -cmneut*fnix(ixt1,iy,1)*recycp(1)* . cmntgpr*(ti(ixt,iy)-eidisspr*ev) - . )/(vpnorm*ennorm*sx(ixt1,iy)) + . ) )/(vpnorm*ennorm*sx(ixt1,iy)) endif # end test for isupgon(1)==1 elseif (ibctipr .eq. 0) then @@ -2630,7 +2748,7 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case flux_inc = 0.5*( fnix(ixt1,iy,1) + fngx(ixt1,iy,1) ) endif endif - t0 = max(tg(ixt1,iy,igsp), temin*ev) + t0 = max(tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = nurlxg * ( fngx(ixt1,iy,igsp) + . fngxrb_use(iy,igsp,jx) - @@ -2639,7 +2757,7 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case . / (vpnorm*n0g(igsp)*sx(ixt1,iy)) elseif (recyrb(iy,igsp,jx) <= 0. .and. . recyrb(iy,igsp,jx) >= -1.) then # recyrb is albedo - t0 = max(tg(ixt,iy,igsp), temin*ev) + t0 = max(tg(ixt,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = nurlxg*( fngx(ixt1,iy,igsp) - . (1+recyrb(iy,igsp,jx))*ng(ixt,iy,igsp)*vxn*sx(ixt1,iy) ) @@ -2702,7 +2820,7 @@ call sputchem (isch_sput(igsp),eng_sput,tvplatrb(iy,jx), endif if (isph_sput(igsp) .eq. 3) then # add chem sput from h. neut do igsp2 = 1, 1+ishymol - t0p = max(tg(ixt1,iy,igsp2),temin*ev) + t0p = max(tg(ixt1,iy,igsp2),tgmin*ev) flx_incid=ng(ixt,iy,igsp2)*.25* . sqrt(8*t0p/(pi*mg(igsp2))) call sputchem (isch_sput(igsp),t0p/ev,tvplatrb(iy,jx), @@ -2714,7 +2832,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatrb(iy,jx), endif if (sputtrb(iy,igsp,jx) .ge. 0. .or. . abs(sputflxrb(iy,igsp,jx)).gt.0.) then - t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) zflux = - sputtrb(iy,igsp,jx) * hflux - . sputflxrb(iy,igsp,jx) - @@ -2724,7 +2842,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatrb(iy,jx), yldot(iv) = nurlxg * (fngx(ixt1,iy,igsp) - zflux) / . (n0(igsp) * vpnorm * sx(ixt1,iy)) elseif (sputtrb(iy,igsp,jx).ge.-9.9) then # neg. sputtrb ==> albedo - t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), temin*ev) + t0 = max(cdifg(igsp)*tg(ixt1,iy,igsp), tgmin*ev) vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) yldot(iv) = nurlxg*( fngx(ixt1,iy,igsp) - . (1+sputtrb(iy,igsp,jx))*ng(ixt,iy,igsp)*vxn*sx(ixt1,iy) ) @@ -2757,6 +2875,38 @@ call xerrab("**INPUT ERROR: istgrb=2 grad opt not implemented") yldot(iv) = nurlxg*( fegx(ixt1,iy,igsp) - 2*cgengmpl* . ng(ixt1,iy,igsp)*vxn*t0*sx(ixt1,iy) )/ . (sx(ixt1,iy)*vpnorm*ennorm) + elseif (istgrb(igsp) == 4) + if (isupgon(igsp)==1) then + if (recyrb(iy,igsp,jx) .gt. 0.) then + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) + fng_alb = (1-albrb(iy,igsp,jx))*ng(ixt1,iy,igsp) + . *vxn*sx(ixt1,iy) + yldot(iv) = nurlxg * ( fegx(ixt1,iy,igsp) + . -cfalbedo*fng_alb*t0 + . +recyrb(iy,igsp,jx)*(1.-cfdiss) + . *fnix(ixt1,iy,1) + . *recyce*cfalbedo + . *( kappar(iy,jx)*zi(1)*te(ixt,iy) + . +ti(ixt,iy) ) )/ + . (vpnorm*ennorm*sx(ixt1,iy)) + elseif (recyrb(iy,igsp,jx) <= 0. .and. + . recyrb(iy,igsp,jx) >= -1.) then # recyrb is albedo + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + fng_alb = (1+recyrb(iy,igsp,jx))*ng(ixt1,iy,igsp) + . *vyn*sx(ixt1,iy) + yldot(iv) = nurlxg * ( fegx(ixt1,iy,igsp) + . -cfalbedo*fng_alb*t0 )/ + . (vpnorm*ennorm*sx(ixt1,iy)) + elseif (recyrb(iy,igsp,jx) < -1.) then #..half Maxwellian + t0 = max(tg(ixt1,iy,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + yldot(iv) = nurlxg * ( fegx(ixt1,iy,igsp) + . -cfalbedo*fnix(ixt1,iy,iigsp)*t0 )/ + . (vpnorm*ennorm*sx(ixt1,iy)) + endif + endif else call xerrab("**INPUT ERROR: istgrb set to unknown option") endif @@ -3411,6 +3561,9 @@ c This subroutine sets components of the array iseqalg(neq) to be unity if (isngonxy(ix,0,igsp).eq.1) then iseqalg(idxg(ix,0,igsp)) = 1 endif + if (istgonxy(ix,0,igsp).eq.1) then + iseqalg(idxtg(ix,0,igsp)) = 1 + endif 9 continue if (isphionxy(ix,0).eq.1) then @@ -3443,6 +3596,9 @@ c This subroutine sets components of the array iseqalg(neq) to be unity if (isngonxy(0,iy,igsp).eq.1) then iseqalg(idxg(0,iy,igsp)) = 1 endif + if (istgonxy(0,iy,igsp).eq.1) then + iseqalg(idxtg(0,iy,igsp)) = 1 + endif 13 continue do 15 ix = 1-ixmnbcl, nx if(isphionxy(ix,iy)==1) iseqalg(idxphi(ix,iy)) = 1 @@ -3473,6 +3629,9 @@ c This subroutine sets components of the array iseqalg(neq) to be unity if (isngonxy(nx+1,iy,igsp).eq.1) then iseqalg(idxg(nx+1,iy,igsp)) = 1 endif + if (istgonxy(nx+1,iy,igsp).eq.1) then + iseqalg(idxtg(nx+1,iy,igsp)) = 1 + endif 17 continue if (isphionxy(nx+1,iy).eq.1) then iseqalg(idxphi(nx+1,iy)) = 1 @@ -3502,6 +3661,9 @@ c This subroutine sets components of the array iseqalg(neq) to be unity if (isngonxy(ix,ny+1,igsp).eq.1) then iseqalg(idxg(ix,ny+1,igsp)) = 1 endif + if (istgonxy(ix,ny+1,igsp).eq.1) then + iseqalg(idxtg(ix,ny+1,igsp)) = 1 + endif 22 continue if (isphionxy(ix,ny+1).eq.1) then iseqalg(idxphi(ix,ny+1)) = 1 @@ -3556,6 +3718,10 @@ c This subroutine sets components of the array iseqalg(neq) to be unity iseqalg(idxg(nxc,iy,igsp)) = 1 iseqalg(idxg(nxc+1,iy,igsp)) = 1 endif + if (istgonxy(nxc,iy,igsp)*istgonxy(nxc+1,iy,igsp)==1) then + iseqalg(idxtg(nxc,iy,igsp)) = 1 + iseqalg(idxtg(nxc+1,iy,igsp)) = 1 + endif enddo if (isphionxy(nxc,iy)*isphionxy(nxc+1,iy)==1) then iseqalg(idxphi(nxc,iy)) = 1 @@ -3596,6 +3762,10 @@ c This subroutine sets components of the array iseqalg(neq) to be unity iseqalg(idxg(ix_lim,iy,igsp)) = 1 iseqalg(idxg(ix_lim+1,iy,igsp)) = 1 endif + if (istgonxy(ix_lim,iy,igsp)*istgonxy(ix_lim+1,iy,igsp)==1) then + iseqalg(idxtg(ix_lim,iy,igsp)) = 1 + iseqalg(idxtg(ix_lim+1,iy,igsp)) = 1 + endif enddo if (isphionxy(ix_lim,iy)*isphionxy(ix_lim+1,iy)==1) then iseqalg(idxphi(ix_lim,iy)) = 1 @@ -3631,6 +3801,10 @@ c This subroutine sets components of the array iseqalg(neq) to be unity . iseqalg(idxg(ixrb(1)+1,iy,igsp)) = 1 if (isngonxy(ixlb(2),iy,igsp)==1) . iseqalg(idxg(ixlb(2),iy,igsp)) = 1 + if (istgonxy(ixrb(1)+1,iy,igsp)==1) + . iseqalg(idxtg(ixrb(1)+1,iy,igsp)) = 1 + if (istgonxy(ixlb(2),iy,igsp)==1) + . iseqalg(idxtg(ixlb(2),iy,igsp)) = 1 enddo if (isphionxy(ixrb(1)+1,iy)==1) . iseqalg(idxphi(ixrb(1)+1,iy)) = 1 @@ -3892,7 +4066,7 @@ subroutine wsmodi(ig) if (jsor > 0) then do ixt = issori(isor), iesori(isor) if(ixt <= ixpt1(1) .or. ixt > ixpt2(1)) then - t0 = max(tg(ixt,0,ig),temin*ev) + t0 = max(tg(ixt,0,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here iwalli(isor) = iwalli(isor) + qe*(1 - albedoi(ixt,ig))* . ng(ixt,1,ig)*vyn*sy(ixt,0) @@ -3901,7 +4075,7 @@ subroutine wsmodi(ig) elseif (jsor < 0) then # signal for outer bdry;now change to jsor > 0 jsor = -jsor do ixt = issoro(isor), iesoro(isor) - t0 = max(tg(ixt,ny+1,ig),temin*ev) + t0 = max(tg(ixt,ny+1,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here iwalli(isor) = iwalli(isor) + qe*(1 - albedoo(ixt,ig))* . ng(ixt,ny,ig)*vyn*sy(ixt,ny) @@ -3911,7 +4085,7 @@ subroutine wsmodi(ig) if (abs(jsor) > 0) then do ixt = issori(jsor), iesori(jsor) if(ixt <= ixpt1(1) .or. ixt > ixpt2(1)) then - t0 = max(tg(ixt,0,ig),temin*ev) + t0 = max(tg(ixt,0,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here fngysi(ixt,ig)= iwalli(isor)*cplsori(jsor)*fwsori(ixt,jsor) ccc . - (1 - albedoi(ixt,ig))* @@ -3964,7 +4138,7 @@ subroutine wsmodo(ig) iwallo(isor) = 0. if (jsor > 0) then do ixt = issoro(isor), iesoro(isor) - t0 = max(tg(ixt,ny+1,ig),temin*ev) + t0 = max(tg(ixt,ny+1,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here iwallo(isor) = iwallo(isor) + qe*(1 - albedoo(ixt,ig))* . ng(ixt,ny,ig)*vyn*sy(ixt,ny) @@ -3972,7 +4146,7 @@ subroutine wsmodo(ig) elseif (jsor < 0) then # signal for inner bdry;now change to jsor > 0 jsor = -jsor do ixt = issori(isor), iesori(isor) - t0 = max(tg(ixt,0,ig),temin*ev) + t0 = max(tg(ixt,0,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here iwallo(isor) = iwallo(isor) + qe*(1 - albedoi(ixt,ig))* . ng(ixt,1,ig)*vyn*sy(ixt,0) @@ -3981,7 +4155,7 @@ subroutine wsmodo(ig) c... reinject current at coupled source jsor if (abs(jsor) > 0) then do ixt = issoro(jsor), iesoro(jsor) - t0 = max(tg(ixt,ny+1,ig),temin*ev) + t0 = max(tg(ixt,ny+1,ig),tgmin*ev) vyn = 0.25 * sqrt( 8*t0/(pi*mg(ig)) ) # mass only for scaling here fngyso(ixt,ig)=iwallo(isor)*cplsoro(jsor)*fwsoro(ixt,jsor) ccc . - (1 - albedoo(ixt,ig))* diff --git a/bbb/convert.m b/bbb/convert.m index 0b4e0532..1dd8f7ed 100755 --- a/bbb/convert.m +++ b/bbb/convert.m @@ -121,6 +121,7 @@ enddo do igsp = 1, ngsp if(istgonxy(ix,iy,igsp) .eq. 1) then + istgcon(igsp) = -1. iv = iv + 1 ntemp = ng(ix,iy,igsp) if(isflxvar == 0) ntemp=n0g(igsp) @@ -300,7 +301,7 @@ c_mpi integer status(MPI_STATUS_SIZE) if(isflxvar == 0) ntemp = n0g(igsp) tg(ix,iy,igsp) = yl(idxtg(ix,iy,igsp))*ennorm/ . (1.5*ntemp) - tg(ix,iy,igsp) = max(tg(ix,iy,igsp), temin*ev) + tg(ix,iy,igsp) = max(tg(ix,iy,igsp), tgmin*ev) endif 65 continue ntemp = nit(ix,iy) + cngtgx(1)*ng(ix,iy,1) @@ -514,6 +515,8 @@ subroutine convsr_aux (ixl, iyl) do 12 iy = js, je do 11 ix = is, ie pri(ix,iy,ifld) = ni(ix,iy,ifld) * ti(ix,iy) + if (ifld.eq.iigsp .and. istgon(1).eq.1) pri(ix,iy,ifld) = + . ni(ix,iy,ifld) * tg(ix,iy,1) if (zi(ifld).ne.0.) then pr(ix,iy) = pr(ix,iy) + pri(ix,iy,ifld) zeff(ix,iy)=zeff(ix,iy)+zi(ifld)**2*ni(ix,iy,ifld) diff --git a/bbb/geometry.m b/bbb/geometry.m index 4cb70d61..5d98807b 100755 --- a/bbb/geometry.m +++ b/bbb/geometry.m @@ -786,7 +786,8 @@ call set_indirect_address(0) c... For 1D poloidal flux-tube geometry (nysol=1), reset cell vertex c... locations near ix=ixpt2 on iy=1 surface if nxpt2msh or nxpt2psh > 1 - if (nxpt2msh+nxpt2psh > 0) call reset1dmeshpt + #if (nxpt2msh+nxpt2psh > 0) call reset1dmeshpt + if (nxpt2msh+nxpt2psh > 0 .or. nxpt1msh+nxpt1psh > 0) call reset1dmeshpt * Define guard cells around the edge of the mesh -- call guardc @@ -2889,10 +2890,14 @@ real xn(0:nx+1,0:ny+1), xvn(0:nx+1,0:ny+1), Use(Share) # nxomit,nxpt2msh,nxpt2psh,rxpt2msh,zxpt2psh c... local scalars - integer ix,ix2 + integer ix,ix2,ix1 real lenginc,theta - ix2 = nxomit+ixpt2(1) #consider outer flux-surface only + if (nxomit .ge. 0) then + ix2 = nxomit+ixpt2(1) #consider outer flux-surface only + else + ix2 = ixpt2(1) + endif c... First fix iy=1 cell vertices above X-point (ixpt2) theta = rxpt2msh/zxpt2msh #defines curvature of flux surface @@ -2912,6 +2917,32 @@ real xn(0:nx+1,0:ny+1), xvn(0:nx+1,0:ny+1), zm(ix+1,1,3) = zm(ix,1,4) enddo +#.. the inner X-pt +#.. +#.. + if (nxomit .ge. 0) then + ix1 = nxomit+ixpt1(1)+1 + else + ix1 = ixpt1(1)+1 + endif +c... First fix iy=1 cell vertices above X-point (ixpt2) + theta = rxpt1msh/zxpt1msh #defines curvature of flux surface + do ix = ix1, ix1+nxpt1msh-1 + rm(ix,1,4) = rm(ix,1,3) + zm(ix,1,4) = zm(ix,1,3) + zxpt1msh + rm(ix-1,1,3) = rm(ix,1,4) + zm(ix-1,1,3) = zm(ix,1,4) + enddo + +c... Second, fix iy=1 cell vertices below X-point + theta = rxpt1psh/zxpt1psh #defines curvature of flux surface + do ix = ix1-1, ix1-nxpt1psh, -1 + rm(ix,1,3) = rm(ix,1,4) + zm(ix,1,3) = zm(ix,1,4) - zxpt1psh + rm(ix+1,1,4) = rm(ix,1,3) + zm(ix+1,1,4) = zm(ix,1,3) + enddo + return end c*** End of subroutine reset1dmeshpt **** diff --git a/bbb/oderhs.m b/bbb/oderhs.m index c15bf2f2..a7191def 100755 --- a/bbb/oderhs.m +++ b/bbb/oderhs.m @@ -740,6 +740,9 @@ call xerrab('Only ismcnon=0 is validated with openmp. ixmp2 = ixpt1(1) + nxomit + 3*(ixpt2(1)-ixpt1(1))/4 endif +c set switches for adding Tg from a default Ti case + if ((istiexclg .eq. 1) .and. (istiinclg_test.ne.1)) cftiexclg = 0.0 + c Set switches for neutrals-related source terms in plasma equations c (MER 1996/10/28) c (IJ 2015/04/06) add ismcnon>=3 for external call to run_neutrals @@ -1969,7 +1972,7 @@ ccc if(isphion+isphiofft .eq. 1) call calc_currents psorgc(ix,iy,igsp) = -ng(ix,iy,igsp)*nuiz(ix,iy,igsp)*vol(ix,iy) + . psorbgg(ix,iy,igsp) psorc(ix,iy,ifld) = - psorgc(ix,iy,igsp) - psordis(ix,iy) = psorc(ix,iy,1) # changed below if ishymol=1 + psordis(ix,iy) = cfdiss*psorc(ix,iy,1) # overwritten below if ishymol=1 psorxrc(ix,iy,ifld) = -ni(ix,iy,ifld)*nurc(ix,iy,igsp)*vol(ix,iy) psorrgc(ix,iy,igsp) = -psorxrc(ix,iy,ifld) msor(ix,iy,ifld) = 0. @@ -2327,8 +2330,10 @@ call xerrab('*** nhgsp must exceed 1 for ishymol=1 ***') endif do iy = iys1, iyf6 do ix = ixs1, ixf6 - nuiz(ix,iy,2) = ne(ix,iy) * ( - . svdiss( te(ix,iy) ) + sigvi_floor ) + nuiz(ix,iy,2) = ne(ix,iy) * ( #.. Tom + . svdiss( te(ix,iy) ) + . + cfizmol*rsa(te(ix,iy),ne_sgvi,rtau(ix,iy),0) + . + sigvi_floor ) massfac = 16*mi(1)/(3*(mg(2)+mi(1))) nuix(ix,iy,2)= fnuizx*nuiz(ix,iy,2) + . massfac*( kelighi(2)*ni(ix,iy,1)+ @@ -2581,17 +2586,19 @@ ccc if (isupgon .eq. 1 .and. zi(ifld) .eq. 0.0) call neudif c ... Now include seic contribution from hydrogen atoms if isupgon=1 c ... Then "ion" species 2 (redundant as gas species 1) is hydr atom - if(isupgon(1)==1 .and. zi(2)<1.e-20 .and. istgon(1)==0) then + if(isupgon(1)==1 .and. zi(2)<1.e-20) then # .and. istgon(1)==0) then if(cfvgpx(2) > 0.) then do iy = j2, j5 do ix = i2, i5 ix1 = ixm1(ix,iy) ix2 = ixp1(ix,iy) iy1 = max(0,iy-1) - seic(ix,iy) = seic(ix,iy) + 0.5*cfvgpx(2)*( - , uuxg(ix, iy,1)*gpix(ix,iy,2) + + seic(ix,iy) = seic(ix,iy) + cftiexclg + . *0.5*cfvgpx(2)*( + . uuxg(ix, iy,1)*gpix(ix,iy,2) + . uuxg(ix1,iy,1)*gpix(ix1,iy,2) )*vol(ix,iy) - seic(ix,iy) = seic(ix,iy) + 0.5*cfvgpy(2)*( + seic(ix,iy) = seic(ix,iy) + cftiexclg + . *0.5*cfvgpy(2)*( . vyg(ix, iy,1)*gpiy(ix,iy,2) + . vyg(ix,iy1,1)*gpiy(ix,iy1,2) )*vol(ix,iy) enddo @@ -2607,7 +2614,7 @@ ccc if (isupgon .eq. 1 .and. zi(ifld) .eq. 0.0) call neudif . ave(gx(ix2,iy),gx(ix,iy))*tv . + up(ix1,iy,2)*rrv(ix1,iy)* . ave(gx(ix,iy),gx(ix1,iy))*t1 ) - seic(ix,iy) = seic(ix,iy) + t1*vol(ix,iy) + seic(ix,iy) = seic(ix,iy) + cftiexclg*t1*vol(ix,iy) enddo enddo endif #test on cfvgpx(2) > 0 or = 0 @@ -2632,7 +2639,7 @@ ccc if (isupgon .eq. 1 .and. zi(ifld) .eq. 0.0) call neudif do 937 ix = i1,i6 c ix1 = ixm1(ix,iy) - vtn = sqrt(max(tg(ix,iy,1),temin*ev)/mi(ifld)) + vtn = sqrt(max(tg(ix,iy,1),tgmin*ev)/mi(ifld)) qfl = flalfvgxa(ix)*nm(ix,iy,ifld)*vtn**2 if(isvisxn_old == 1) then lmfpn = 1./(sigcx * @@ -2667,7 +2674,7 @@ ccc if (isupgon .eq. 1 .and. zi(ifld) .eq. 0.0) call neudif tgupyface = 0.25*( tg(ix,iy,1)+ . tg(ix,iyp1,1)+tg(ix2,iy,1)+ . tg(ix3,iyp1,1) ) - vtn = sqrt(max(tgupyface,temin*ev)/mi(ifld)) + vtn = sqrt(max(tgupyface,tgmin*ev)/mi(ifld)) nmxface = 0.5*(nm(ix,iy,ifld)+nm(ix2,iy,ifld)) ngupyface = 0.25*( ni(ix,iy,ifld)+ . ni(ix,iyp1,ifld)+ni(ix2,iy,ifld)+ @@ -2997,11 +3004,11 @@ ccc iysptrx is the last closed flux surface (see S.R. nphygeo) c..1dn0802 c IJ 2016/10/10 add cfneutsor_ei multiplier to control fraction of neutral energy to add hcxi(ix,iy) = hcxi(ix,iy) - . + cfneut*cfneutsor_ei*kxn*( ng(ix ,iy,1)*ti(ix ,iy) + . + cftiexclg*cfneut*cfneutsor_ei*kxn*( ng(ix ,iy,1)*ti(ix ,iy) . +ng(ix1,iy,1)*ti(ix1,iy) ) / . (mi(1)*(nucx(ix,iy,1) + nucx(ix1,iy,1))) hcyi(ix,iy) = hcyi(ix,iy) - . + cfneut*cfneutsor_ei*kyn*( ngy0(ix,iy,1)*tiy0(ix,iy) + . + cftiexclg*cfneut*cfneutsor_ei*kyn*( ngy0(ix,iy,1)*tiy0(ix,iy) . +ngy1(ix,iy,1)*tiy1(ix,iy) ) / . (mi(1)*(nucx(ix,iy,1) + nucx(ix,iyp1,1))) endif @@ -3039,7 +3046,8 @@ ccc iysptrx is the last closed flux surface (see S.R. nphygeo) qshx = cshx * (tg(ix,iy,1)-tg(ix1,iy,1)) * gxf(ix,iy) hcxn(ix,iy) = cshx / . (1 + (abs(qshx/qflx))**flgamtg)**(1./flgamtg) - hcxi(ix,iy) = hcxi(ix,iy) + cfneut*cfneutsor_ei*hcxn(ix,iy) + hcxi(ix,iy) = hcxi(ix,iy) + + . cftiexclg*cfneut*cfneutsor_ei*hcxn(ix,iy) c Now for the radial flux limit - good for nonorthog grid too qfly = flalftgya(iy) * sqrt(tgavey/mi(iigsp)) * noavey * . tgavey @@ -3049,7 +3057,8 @@ ccc iysptrx is the last closed flux surface (see S.R. nphygeo) qshy = cshy * (tgy0(ix,iy1,1)-tgy1(ix,iy1,1))/dynog(ix,iy) hcyn(ix,iy) = cshy / . (1 + (abs(qshy/qfly))**flgamtg)**(1./flgamtg) - hcyi(ix,iy) = hcyi(ix,iy) + cfneut*cfneutsor_ei*hcyn(ix,iy) + hcyi(ix,iy) = hcyi(ix,iy) + + . cftiexclg*cfneut*cfneutsor_ei*hcyn(ix,iy) c 63 continue 62 continue @@ -3143,6 +3152,10 @@ ccc iysptrx is the last closed flux surface (see S.R. nphygeo) . cfhcygc(igsp)*noavey*kyg_use(ix,iy,igsp) enddo enddo + if (igsp.eq.1 .and. isupgon(igsp).eq.1) then + hcxg(:,:,igsp) = hcxn(:,:) + hcyg(:,:,igsp) = hcyn(:,:) + endif enddo endif @@ -3153,8 +3166,10 @@ ccc iysptrx is the last closed flux surface (see S.R. nphygeo) do iy = j1, j6 do ix = i1, i6 nhi_nha = ni(ix,iy,1)+ni(ix,iy,2) - eqpg(ix,iy,igsp) = cftgeqp*ng(ix,iy,igsp)*nhi_nha* - . keligig(igsp) +# eqpg(ix,iy,igsp) = cftgeqp*ng(ix,iy,igsp)*nhi_nha* +# . keligig(igsp) + eqpg(ix,iy,igsp) = cftgeqp*ng(ix,iy,igsp)* + . (ni(ix,iy,1)+cftiexclg*ni(ix,iy,2))*keligig(igsp) enddo enddo enddo @@ -3571,7 +3586,7 @@ c The convective component is already already added through uu(ix,iy). endif c... Now flux limit with flalfvgxy if ifld=2 if (ifld==2) then - t0 = max(tg(ix,iy,1),temin*ev) + t0 = max(tg(ix,iy,1),tgmin*ev) vtn = sqrt(t0/mg(1)) qfl = flalfvgxya(ix)*0.5*(sx(ix,iy)+sx(ix1,iy))*vtn**2* . nm(ix,iy,ifld) + cutlo @@ -3753,8 +3768,11 @@ c The convective component is already already added through uu(ix,iy). resmo(ix,iy,iigsp) = # TR resmo(ix,iy,ifld) #IJ 2016 . - cmneut * cmneutsor_mi * uesor_up(ix,iy,1) . -sx(ix,iy) * rrv(ix,iy) * - . cpgx*( ni(ix2,iy,iigsp)*ti(ix2,iy)- - . ni(ix,iy,iigsp)*ti(ix,iy) ) + . cpgx*( cftiexclg*(ni(ix2,iy,iigsp)*ti(ix2,iy)- + . ni(ix,iy,iigsp)*ti(ix,iy))+ + . (1.0-cftiexclg)* + . (ni(ix2,iy,iigsp)*tg(ix2,iy,1)- + . ni(ix,iy,iigsp)*tg(ix,iy,1)) ) . -cfupcx*0.25*volv(ix,iy)* . (nucx(ix,iy,1)+nucx(ix2,iy,1))* . (nm(ix,iy,iigsp)+nm(ix2,iy,iigsp))* @@ -3990,21 +4008,21 @@ c if (ifld .ne. iigsp) then do 726 iy = j4, j8 do 725 ix = i1, i5 floxi(ix,iy) = floxi(ix,iy) + - . cfcvti*2.5*cfneut*cfneutsor_ei*fnix(ix,iy,ifld) + . cftiexclg*cfcvti*2.5*cfneut*cfneutsor_ei*fnix(ix,iy,ifld) 725 continue # next correct for incoming neut pwr = 0 do jx = 1, nxpt #if at plate, sub (1-cfloxiplt)*neut-contrib if(ixmnbcl==1) then #real plate-need for parallel UEDGE iixt = ixlb(jx) #left plate if(fnix(iixt,iy,ifld) > 0.) then floxi(iixt,iy) = floxi(iixt,iy) - (1.-cfloxiplt)* - . cfcvti*2.5*cfneut*cfneutsor_ei*fnix(iixt,iy,ifld) + . cftiexclg*cfcvti*2.5*cfneut*cfneutsor_ei*fnix(iixt,iy,ifld) endif endif if(ixmxbcl==1) then #real plate-need for parallel UEDGE iixt = ixrb(jx) # right plate if(fnix(iixt,iy,ifld) < 0.) then floxi(iixt,iy) = floxi(iixt,iy) - (1.-cfloxiplt)* - . cfcvti*2.5*cfneut*cfneutsor_ei*fnix(iixt,iy,ifld) + . cftiexclg*cfcvti*2.5*cfneut*cfneutsor_ei*fnix(iixt,iy,ifld) endif floxi(ixrb(jx)+1,iy) = 0.0e0 #cosmetic endif @@ -4044,7 +4062,7 @@ c if (ifld .ne. iigsp) then do iy = j1, j5 do ix = i4, i8 floyi(ix,iy) = floyi(ix,iy) - . + cfneut * cfneutsor_ei * 2.5 * fniy(ix,iy,ifld) + . + cftiexclg*cfneut * cfneutsor_ei * 2.5 * fniy(ix,iy,ifld) enddo enddo c ... Make correction at walls to prevent recyc neutrals injecting pwr @@ -4052,12 +4070,12 @@ c if (ifld .ne. iigsp) then if (matwallo(ix) > 0 .and. recycwot(ix,1)>0.) then fniy_recy = max(recycwot(ix,1)*fac2sp*fniy(ix,ny,1), 0.) floyi(ix,ny) = floyi(ix,ny) + - . cfneut*cfneutsor_ei*2.5*(1.-cfloygwall)*fniy_recy + . cftiexclg*cfneut*cfneutsor_ei*2.5*(1.-cfloygwall)*fniy_recy endif if (matwalli(ix) > 0 .and. recycwit(ix,1,1)>0.) then fniy_recy = min(recycwit(ix,1,1)*fac2sp*fniy(ix,0,1), 0.) - floyi(ix,0) = floyi(ix,0) + - . cfneut*cfneutsor_ei*2.5*(1.-cfloygwall)*fniy_recy + floyi(ix,0) = floyi(ix,0) + + . cftiexclg*cfneut*cfneutsor_ei*2.5*(1.-cfloygwall)*fniy_recy endif enddo @@ -4183,7 +4201,7 @@ c if (ifld .ne. iigsp) then do 141 iy = j4, j8 do 142 ix = i1, i5 floxi(ix,iy) = floxi(ix,iy) + - . cfneut*cfneutsor_ei*cngtgx(1)*cfcvti*2.5*fngx(ix,iy,1) + . cftiexclg*cfneut*cfneutsor_ei*cngtgx(1)*cfcvti*2.5*fngx(ix,iy,1) 142 continue floxi(nx+1,iy) = 0.0e0 141 continue @@ -4193,7 +4211,7 @@ c if (ifld .ne. iigsp) then do 145 iy = j1, j5 do 144 ix = i4, i8 floyi(ix,iy) = floyi(ix,iy) - . + cfneut*cfneutsor_ei*cngtgy(1)*2.5*fngy(ix,iy,1) + . + cftiexclg*cfneut*cfneutsor_ei*cngtgy(1)*2.5*fngy(ix,iy,1) 144 continue 145 continue @@ -4321,7 +4339,7 @@ c if (ifld .ne. iigsp) then . (log(ti(ix2,iy)) + log(ti(ix,iy))) )* . ( (fcdif*kyi+kyi_use(ix,iy))*0.5* . (nit(ix2,iy)+nit(ix,iy)) - . + cfneut*cfneutsor_ei*0.25*(hcyn(ix ,iy)+hcyn(ix ,iy1) + . + cftiexclg*cfneut*cfneutsor_ei*0.25*(hcyn(ix ,iy)+hcyn(ix ,iy1) . +hcyn(ix2,iy)+hcyn(ix4,iy1)) ) * . ( grdnv/cos(angfx(ix,iy)) . - (log(ti(ix2,iy)) - log(ti(ix,iy)))* @@ -4331,8 +4349,8 @@ c if (ifld .ne. iigsp) then t1 = max(ti(ix2,iy),temin*ev) vttn = t0*sqrt( t0/mi(1) ) vttp = t1*sqrt( t1/mi(1) ) - qfl = flalftxy * 0.125 * sx(ix,iy) * (vttn + vttp) * - . (ni(ix,iy,1)+ng(ix,iy,1)+ni(ix2,iy,1)+ng(ix2,iy,1)) + qfl = flalftxy * (cftiexclg*0.125+(1.-cftiexclg)*0.25) * sx(ix,iy) * (vttn + vttp) * + . (ni(ix,iy,1)+cftiexclg*ng(ix,iy,1)+ni(ix2,iy,1)+cftiexclg*ng(ix2,iy,1)) feixy(ix,iy) = feixy(ix,iy) / . sqrt(1. + (feixy(ix,iy)/qfl)**2) @@ -4418,8 +4436,10 @@ c if (ifld .ne. iigsp) then seg_ue(ix,iy,jfld)=-( (fegx_ue(ix,iy,jfld)-fegx_ue(ix1,iy, jfld)) . + fluxfacy*(fegy_ue(ix,iy,jfld)-fegy_ue(ix, iy-1,jfld)) ) . *( (ni(ix,iy,jfld)*ti(ix,iy))/(ni(ix,iy,jfld)*ti(ix,iy)) ) - resei(ix,iy) = resei(ix,iy) + - . cmneutdiv*cmneutdiv_feg*seg_ue(ix,iy,jfld) + resei(ix,iy) = resei(ix,iy) + + . cftiexclg*cmneutdiv*cmneutdiv_feg*seg_ue(ix,iy,jfld) + reseg(ix,iy,1) = reseg(ix,iy,1) + + . cmneutdiv*cmneutdiv_feg*seg_ue(ix,iy,jfld) endif 309 continue 310 continue @@ -4543,11 +4563,36 @@ ccc call volave(nx, ny, j2, j5, i2, i5, ixp1(0,0), ixm1(0,0), c to the friction force between neutrals and ions t1 = 0.5*(up(ix,iy,1)+up(ix1,iy,1)) t2 = 0.5*(up(ix,iy,iigsp)+up(ix1,iy,iigsp)) + temp3 = cfnidhgy*0.25*(vy(ix,iy,iigsp)+vy(ix1,iy,iigsp)) + . *(vy(ix,iy,iigsp)+vy(ix1,iy,iigsp)) + temp4 = cfnidhg2*0.25*(v2(ix,iy,iigsp)+v2(ix1,iy,iigsp)) + . *(v2(ix,iy,iigsp)+v2(ix1,iy,iigsp)) + tv = cfticx*nucx(ix,iy,1)*ng(ix,iy,1)*vol(ix,iy) + t0 = 1.5*( tg(ix,iy,1)* (psor(ix,iy,1)+tv) + . -ti(ix,iy) * (psorrg(ix,iy,1)+tv) ) resei(ix,iy) = resei(ix,iy) + w0(ix,iy) - . + cfneut * cfneutsor_ei * cfnidh * 0.5*mi(1) * (t1-t2)*(t1-t2) * - . ( psor(ix,iy,1) + psorrg(ix,iy,1) - . + 2*cfticx*nucx(ix,iy,1)*ng(ix,iy,1)*vol(ix,iy) ) - . + cfneut * cfneutsor_ei * cnsor*eion*ev*psordis(ix,iy) + . + cfneut * cfneutsor_ei * cfnidh * 0.5*mi(1) * + . ( (t1-t2)*(t1-t2)+temp3+temp4 ) * + . ( psor(ix,iy,1) + cftiexclg*psorrg(ix,iy,1) + . + tv + cftiexclg * tv ) + . + (1.0-cftiexclg) * t0 + . + cftiexclg * cfneut * cfneutsor_ei * cnsor + . *( eion*ev+cfnidhdis* + . 0.5*mg(1)*(t2*t2+temp3+temp4) )*psordis(ix,iy) + . + cfnidh2* + . ( -mi(1)*t1*t2*(psor(ix,iy,1)+tv) + . +0.5*mi(1)*t1*t1* + . (psor(ix,iy,1)+psorrg(ix,iy,1)+2*tv) ) + reseg(ix,iy,1) = reseg(ix,iy,1) + pwrsorg(ix,iy,1) + . - t0+0.5*mg(1) * ( (t1-t2)*(t1-t2) + . +temp3+temp4 ) + . * (psorrg(ix,iy,1)+tv) + . + ( eion*ev + cfnidh*cfnidhdis* + . 0.5*mg(1)*(t2*t2+temp3+temp4) )*psordis(ix,iy) + . + cfnidh2* + . ( -mg(1)*t1*t2*(psorrg(ix,iy,1)+tv) + . +0.5*mg(1)*(t2*t2+temp3+temp4)* + . (psor(ix,iy,1)+psorrg(ix,iy,1)+2*tv) ) else resei(ix,iy) = resei(ix,iy) + w0(ix,iy) . + cfneut * cfneutsor_ei * ctsor*1.25e-1*mi(1)* @@ -4564,14 +4609,17 @@ ccc call volave(nx, ny, j2, j5, i2, i5, ixp1(0,0), ixm1(0,0), c ... If molecules are present as gas species 2, add ion/atom cooling - if(ishymol == 1) then - do iy = j2, j5 - do ix = i2, i5 - resei(ix,iy) = resei(ix,iy) - vol(ix,iy)*eqpg(ix,iy,2)* - . (ti(ix,iy)-tg(ix,iy,2)) - enddo - enddo - endif + # energy transfer between ions and molecueles due to + # ion/molecule elastic collisions have been moved in + # engbalg subroutine, so comment the following lines... +# if(ishymol == 1) then +# do iy = j2, j5 +# do ix = i2, i5 +# resei(ix,iy) = resei(ix,iy) - vol(ix,iy)*eqpg(ix,iy,2)* +# . (ti(ix,iy)-tg(ix,iy,2)) +# enddo +# enddo +# endif * -- Energy transfer to impurity neutrals at tg(,,igsp) if (ngsp >= 2) then # for now, specialized to igsp=2 only @@ -4838,7 +4886,12 @@ cc elseif (ishosor .ne. 0) wvh(ix,iy,ifld) = wvh(ix,iy,ifld) - . sin(thetacc)*cfvcsy(ifld)*cfvisy* . visy(ix,iy,ifld)*dupdx*dupdy - resei(ix,iy) = resei(ix,iy) + wvh(ix,iy,ifld)*vol(ix,iy) + if (zi(ifld)==0.0 .and. ifld.eq.iigsp) then + resei(ix,iy) = resei(ix,iy) + cftiexclg*wvh(ix,iy,ifld)*vol(ix,iy) + reseg(ix,iy,1) = reseg(ix,iy,1) + wvh(ix,iy,ifld)*vol(ix,iy) + else + resei(ix,iy) = resei(ix,iy) + wvh(ix,iy,ifld)*vol(ix,iy) + endif 155 continue # loop over up species ifld 156 continue 157 continue @@ -6017,7 +6070,8 @@ c If we are solving for the neutral parallel momentum (isupgon=1) Use(MCN_dim) # ngsp, ... Use(MCN_sources) # cfneut_sng, cfneutdiv_fng, ... mcfngx, mcfngy, ... Use(Interp) # ngs, tgs - + Use(Bfield) # rbfbt + * -- procedures -- real ave ave(t0,t1) = 2*t0*t1 / (cutlo+t0+t1) @@ -6152,8 +6206,8 @@ c write (*,*) "neudifpg" do 890 iy = j1, j5 do 889 ix = i4, i8 ngyface = 0.5*(ng(ix,iy,igsp)+ng(ix,iy+1,igsp)) - t0 = max(tg(ix,iy,igsp),temin*ev) - t1 = max(tg(ix,iy+1,igsp),temin*ev) + t0 = max(tg(ix,iy,igsp),tgmin*ev) + t1 = max(tg(ix,iy+1,igsp),tgmin*ev) vtn = sqrt( t0/mg(igsp) ) vtnp = sqrt( t1/mg(igsp) ) nu1 = nuix(ix,iy,igsp) + vtn/lgmax(igsp) @@ -6474,6 +6528,9 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg, do iy = j4, j6 do ix = i1, i6 uu(ix,iy,iigsp) = uug(ix,iy,igsp) + v2(ix,iy,iigsp) = ( uuxg(ix,iy,igsp) + . - up(ix,iy,iigsp)*rrv(ix,iy) ) + . /(rbfbt(ix,iy) + rbfbt(ixp1(ix,iy),iy))*2. enddo enddo endif @@ -7472,7 +7529,7 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg, enddo do igsp = 1, ngsp - if(istgon(igsp) == 1) then + if(istgon(igsp) == 1) then do iy = j2, j5 do ix = i2, i5 ix1 = ixm1(ix,iy) @@ -7481,8 +7538,8 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg, tv = (pg(ix2,iy,igsp) - pg(ix ,iy,igsp)) t1 = (pg(ix ,iy,igsp) - pg(ix1,iy,igsp)) segc(ix,iy,igsp) = 0.5*cvgpg*( - , uug(ix, iy,igsp)*ave(gx(ix2,iy),gx(ix, iy))*tv + - . uug(ix1,iy,igsp)*ave(gx(ix ,iy),gx(ix1,iy))*t1 )* + , uuxg(ix, iy,igsp)*ave(gx(ix2,iy),gx(ix, iy))*tv + + . uuxg(ix1,iy,igsp)*ave(gx(ix ,iy),gx(ix1,iy))*t1 )* . vol(ix,iy) t2 = cvgpg*0.5*( vyg(ix, iy,igsp)*dynog(ix, iy)* . (pgy1(ix, iy,igsp)-pgy0(ix, iy,igsp)) + @@ -7640,8 +7697,8 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp), ix4 = ixp1(ix,iy1) ix5 = ixm1(ix,iy+1) ix6 = ixp1(ix,iy+1) - t0 = max(tg(ix,iy,igsp),temin*ev) - t1 = max(tg(ix2,iy,igsp),temin*ev) + t0 = max(tg(ix,iy,igsp),tgmin*ev) + t1 = max(tg(ix2,iy,igsp),tgmin*ev) vtn = sqrt( t0/mg(igsp) ) vtnp = sqrt( t1/mg(igsp) ) nu1 = nuix(ix,iy,igsp) + vtn/lgmax(igsp) @@ -7673,8 +7730,8 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp), . - (log(tg(ix2,iy,igsp)) - log(tg(ix,iy,igsp)))* . gxf(ix,iy) )*sx(ix,iy) c... Flux limit with flalftxt even though hcys have parallel FL built in - t0 = max(tg(ix,iy,igsp),temin*ev) - t1 = max(tg(ix2,iy,igsp),temin*ev) + t0 = max(tg(ix,iy,igsp),tgmin*ev) + t1 = max(tg(ix2,iy,igsp),tgmin*ev) vttn = t0*sqrt( t0/mg(igsp) ) vttp = t1*sqrt( t1/mg(igsp) ) if(isfegxyqflave == 0) then @@ -7707,8 +7764,49 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp), . fegy(ix,iy,igsp)-fegy(ix, iy1,igsp) ) . + segc(ix,iy,igsp) reseg(ix,iy,igsp)= reseg(ix,iy,igsp) + vol(ix,iy)* - . eqpg(ix,iy,igsp)*(ti(ix,iy)-tg(ix,iy,igsp))+ - . cftgdiss(igsp)*psorg(ix,iy,igsp)*tg(ix,iy,igsp) + . eqpg(ix,iy,igsp)*(ti(ix,iy)-tg(ix,iy,igsp))#+ +# . cftgdiss(igsp)*psorg(ix,iy,igsp)*tg(ix,iy,igsp) + if (igsp.eq.1) then #..for D0, we should include D+ and D0 in Ti + seic(ix,iy) = seic(ix,iy)- vol(ix,iy)*(1.0-cftiexclg)* + . eqpg(ix,iy,igsp)* + . (ti(ix,iy)-tg(ix,iy,igsp)) + else + seic(ix,iy) = seic(ix,iy)- vol(ix,iy)* + . eqpg(ix,iy,igsp)* + . (ti(ix,iy)-tg(ix,iy,igsp)) + reseg(ix,iy,igsp) = reseg(ix,iy,igsp) + . + cftgeqp*ng(ix,iy,igsp)* + . (1.0-cftiexclg)*ng(ix,iy,1)*kelighg(igsp)* + . (tg(ix,iy,1)-tg(ix,iy,igsp))*vol(ix,iy) + reseg(ix,iy,1) = reseg(ix,iy,1) - cftgeqp*ng(ix,iy,igsp)* + . ng(ix,iy,1)*kelighg(igsp)* + . (tg(ix,iy,1)-tg(ix,iy,igsp))*vol(ix,iy) + if (ishymol.eq.1 .and. igsp.eq.2) then #..D2 dissociation + reseg(ix,iy,igsp) = + . reseg(ix,iy,igsp)+psorg(ix,iy,igsp) + . *1.5*tg(ix,iy,igsp) + t0 = cfnidhmol*0.25*(uuxg(ix,iy,igsp)+uuxg(ix1,iy,igsp)) + . *(uuxg(ix,iy,igsp)+uuxg(ix1,iy,igsp)) + t1 = cfnidhmol*0.25*(vyg(ix,iy,igsp)+vyg(ix1,iy,igsp)) + . *(vyg(ix,iy,igsp)+vyg(ix1,iy,igsp)) + t2 = 0. #.. molecule v in the tol direction, it seems to be assumed as 0 in neudifpg? + reseg(ix,iy,1) = reseg(ix,iy,1) + cfnidhdis* + . 0.5*mg(1)*(t0+t1+t2)*psordis(ix,iy) + seic(ix,iy) = seic(ix,iy) + cftiexclg*cfnidhdis* + . 0.5*mg(1)*(t0+t1+t2)*psordis(ix,iy) + t0 = cfnidhmol*0.25*(uuxg(ix,iy,igsp)+uuxg(ix1,iy,igsp)) + . *(uuxg(ix,iy,1)+uuxg(ix1,iy,1)) + t1 = cfnidhmol*0.25*(vyg(ix,iy,igsp)+vyg(ix1,iy,igsp)) + . *(vyg(ix,iy,1)+vyg(ix1,iy,1)) + t2 = 0. + reseg(ix,iy,1) = reseg(ix,iy,1) - cfnidhdis* + . mg(1)*(t0+t1+t2)*psordis(ix,iy) + seic(ix,iy) = seic(ix,iy) - cftiexclg*cfnidhdis* + . mg(1)*(t0+t1+t2)*psordis(ix,iy) + endif + endif + #..zml place holder for neutral-neutral collision, + #.. not included above? enddo enddo enddo @@ -7728,9 +7826,14 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp), do ifld = nhsp+1, nisp do iy = j2, j5 # iys,iyf limits dont seem to work(?) do ix = i2, i5 - resei(ix,iy) =resei(ix,iy) -cftiimpg*1.5*ni(ix,iy,ifld)* + # possible bugs here? resei is replaced with seic + # since resei is not defined before this subroutine called + # more needs to be done here.. e.g. for segc + seic(ix,iy) =seic(ix,iy) -cftiimpg*1.5*ni(ix,iy,ifld)* . (nucxi(ix,iy,ifld)+nueli(ix,iy,ifld))* . (ti(ix,iy) - tg(ix,iy,2))*vol(ix,iy) + #..zml place holder for things remain to be done for segc + #.. to include neutral-impurity ion collisions enddo enddo enddo @@ -7967,9 +8070,12 @@ ccc if (isngonxy(ix,iy,1) .eq. 1) nbidot = cngtgx(1)*yldot(idxg(ix,iy do 255 ifld = 1, nisp if (isnionxy(ix,iy,ifld) .eq. 1) then iv = idxn(ix,iy,ifld) - nbidot = nbidot + yldot(iv)*n0(ifld) + #.. separate ions and atoms + #nbidot = nbidot + yldot(iv)*n0(ifld) if (isupgon(1)==1 .and. zi(ifld)==0) then #neutral hyd nbgdot = yldot(iv)*n0(ifld) + else + nbidot = nbidot + yldot(iv)*n0(ifld) endif nbedot = nbedot + zi(ifld)*yldot(iv)*n0(ifld) endif @@ -8027,22 +8133,43 @@ ccc if (isngonxy(ix,iy,1) .eq. 1) nbidot = cngtgx(1)*yldot(idxg(ix,iy endif if(istionxy(ix,iy) == 1) then iv1 = idxti(ix,iy) - if(iseqalg(iv1) == 0) then - if(isupgon(1)==1) then #atom dens included in nbidot - yldot(iv1) = ( yldot(iv1)*nnorm - yl(iv1)*nbidot )/ - . (nit(ix,iy)+ ni(ix,iy,2)) - else #atom dens not included in nbidot + if (iseqalg(iv1) == 0) then + if(isupgon(1)==1) then + yldot(iv1) = ( yldot(iv1)*nnorm - + . yl(iv1)*(nbidot + cftiexclg*nbgdot) ) / + . (nit(ix,iy) + cftiexclg*ni(ix,iy,2) ) + else yldot(iv1) = ( yldot(iv1)*nnorm - yl(iv1)* . ( nbidot + cngtgx(1)*nbg2dot(1) ) ) / . (nit(ix,iy) + cngtgx(1)*ng(ix,iy,1)) endif endif endif - do igsp = 1, ngsp - if(istgonxy(ix,iy,igsp) == 1) then - iv = idxtg(ix,iy,igsp) - yldot(iv) = ( yldot(iv)*n0g(igsp) - - . yl(iv)*nbg2dot(igsp) )/ng(ix,iy,igsp) + do igsp = 1, ngsp + if (igsp == 1 .and. isupgon(1) == 1) then + if (istgonxy(ix,iy,igsp) == 1) then + iv = idxtg(ix,iy,igsp) + if (iseqalg(iv).eq.0) then +# yldot(iv) = ( yldot(iv)*n0g(igsp) - +# . yl(iv)*yldot(idxn(ix,iy,iigsp))*n0(ifld) ) +# . /ni(ix,iy,iigsp) + #..the above one causes problem when turning off ni(:,:,2) + yldot(iv) = ( yldot(iv)*n0g(igsp) - + . yl(iv)*nbgdot ) / ni(ix,iy,iigsp) + endif + endif + else + if (istgonxy(ix,iy,igsp) == 1) then + iv = idxtg(ix,iy,igsp) + if (iseqalg(iv) == 0) then +# yldot(iv) = ( yldot(iv)*n0g(igsp) - +# . yl(iv)*yldot(idxg(ix,iy,igsp))*n0g(igsp) ) +# . /ng(ix,iy,igsp) + #..the above one causes problem when turning off ng(:,:,2) + yldot(iv) = ( yldot(iv)*n0g(igsp) - + . yl(iv)*nbg2dot(igsp) ) / ng(ix,iy,igsp) + endif + endif endif enddo endif #end if-test on isflxvar @@ -8068,12 +8195,14 @@ ccc if (isngonxy(ix,iy,1) .eq. 1) nbidot = cngtgx(1)*yldot(idxg(ix,iy # z0pe,z0pi,r0pe,r0pi,zwpe,zwpi,rwpe,rwpi, # z0ni,r0ni,zwni,rwni,voljcsor,jcvsor, # ix_sjcsor, ix_ejcsor, iy_sjcsor, iy_ejcsor, - # thetarot,rcutmin,zcutmin,effvng + # thetarot,rcutmin,zcutmin,effvng, + # pwrsorg #.. manufactured solution Use(Phyvar) # ev Use(Bcond) # islimsor,rlimiter Use(Parallv) # nxg,nyg Use(Xpoint_indices) # ixpt1,ixpt2,iysptrx Use(Share) # nxomit + Use(UEpar) # ismanufactured #..zml manufactured solution * -- local scalars -- real effvni, effvup, effvpe, effvpi, effvjel, zc, rc, ivolcurt, @@ -8108,9 +8237,11 @@ call s2fill (nx+2, ny+2, 0., volmsor(0:nx+1,0:ny+1,ifld), 1, nx+2) mvolcurt = mvolcurt + mvolcur(ifld) enddo do igsp = 1, ngsp - call s2fill (nx+2, ny+2, 0., volpsor(0:nx+1,0:ny+1,igsp), 1, nx+2) - call s2fill (nx+2, ny+2, 0., volmsor(0:nx+1,0:ny+1,igsp), 1, nx+2) + call s2fill (nx+2, ny+2, 0., volpsorg(0:nx+1,0:ny+1,igsp), 1, nx+2) ivolcurgt = ivolcurgt + ivolcurg(igsp) + #..zml manufactured solution + if (ismanufactured(igsp) .ne. 1) + . call s2fill (nx+2, ny+2, 0., pwrsorg(0:nx+1,0:ny+1,igsp), 1, nx+2) enddo cccMER NOTE: generalize the following for multiple x-points cc Define index ranges for a localized ion-loss sink; crude & temporary diff --git a/com/com.v b/com/com.v index af44091d..8d737898 100755 --- a/com/com.v +++ b/com/com.v @@ -276,7 +276,13 @@ zxpt2msh real /0.01/ [m] # spacing iy=1 1D vertices reset above ixpt2 rxpt2msh real /0.0/ [m] # spacing iy=1 1D vertices reset above ixpt2 zxpt2psh real /0.01/ [m] # spacing iy=1 1D vertices reset below ixpt2 rxpt2psh real /0.0/ [m] # spacing iy=1 1D vertices reset below ixpt2 -ismpsym integer /0/ +gridgen # =1 re-constructs "guard" cells at midplane +nxpt1msh integer /0/ # number of 1D cells above ixpt1 reset/no-cross +nxpt1psh integer /0/ # number of 1D cells below ixpt1 reset/no-cross +zxpt1msh real /0.01/ [m] # spacing iy=1 1D vertices reset above ixpt1 +rxpt1msh real /0.0/ [m] # spacing iy=1 1D vertices reset above ixpt1 +zxpt1psh real /0.01/ [m] # spacing iy=1 1D vertices reset below ixpt1 +rxpt1psh real /0.0/ [m] # spacing iy=1 1D vertices reset below ixpt1 +ismpsym integer /0/ +girdgen # =1 re-constructs "guard" cells at midplane # of "dnbot" via up/down symmetry isudsym integer /0/ +input #=1 up-down symmetric setup (only down part is modeled) islimon integer /0/ +input diff --git a/pyscripts/__src__.py b/pyscripts/__src__.py index 3c2c2db5..a0be297b 100644 --- a/pyscripts/__src__.py +++ b/pyscripts/__src__.py @@ -1 +1 @@ -__src__ = '/home/meyer8/UEDGE' +__src__ = '/home/meyer8/UEDGE' \ No newline at end of file