主程序：1.首先设置子程序运行的条件，if(t.gt.t1inj) call chem
           2.chem.f子程序中:
               do i4=ifirst, ncells
               ...........
               if(issot.eq.1) call issot(spdg,........,dnuc(i4))
               .......
               enddo
           3.issot.f子程序为：
      tijk=temp
      if(tijk.lt.tcuts) return
      rtijk=1./tijk
      rgast=rgas*tijk
      davg=0.5*(drncr+dnuc)
      expe2=exp(-ee2*rtijk)
      expe3=exp(-ee3*rtijk)
      xmpn=0.0
      xnrn=0.0
      vovind=0.0
 twlvth=1.0/12.0
c +++
c +++ 1st, the formation of radical nuclei (surovikin eq. 1, p. 93)
c +++
      hcn=spd(kfuel)*ramummw
      surv0=hcn*1.e+13*exp(-ee1*rtijk)
      rnmass=pio6*rhorn*drn**3
      xnrn=spd(krn)/rnmass
      pnmass=pio6*rhopn*davg**3
      xmpn=spd(kpn)/pnmass
c +++

      fma1=pio6*drn**3*avogad*rhorn
      fm1star=0.5*fma1
      surg0=0.5 * pi*drn**2*sqrt(8.*rgast/(pi*fm1star))
      sqrtk=sqrt(rgast/(pi2*fma1))*exp(-ee4*rtijk)
      surk=pi*davg**2*xnrn*xmpn*sqrtk
      xnrn=xnrn+dtc*(surv0 + surf*xnrn - surg0*xnrn**2 - surk)
c +++
c +++ 2nd, the growth of radical nuclei (surovikin eq. 3, p. 95)

      fm2star=amw(kfuel)*fma1/(amw(kfuel)+fma1)
      rate=(ra+0.5*drn)**2*sqrt(8.*rgast/(pi*fm2star))*hcn*expe2
      drn=drn+dtc*(2.*xmc*xnc/(rhorn*drn**2)*rate)
c +++
c +++ 3rd, the rate of consumption of hydrocarbon (surovikin eq.5, p.96)

      sqrtrm=sqrt(rgast*rpi2mw)
      dspdhcp=hcn*expe3*sqrtrm*xmpn*pi*davg**2*amummw
      dspdhcrn=rate*xnrn*pi*amummw
      spd(kfuel)=spd(kfuel)-dtc*(dspdhcp+dspdhcrn)
      if(drn.lt.drncr) go to 10
c +++
c +++ 4th, the oxidation of soot particles (ref. haynes & wagner, 1981)

      hwka=fka*exp(-eka*rtijk)
      hwkb=fkb*exp(-ekb*rtijk)
      hwkt=fkt*exp(-ekt*rtijk)
      hwkz=fkz*exp(ekz*rtijk)
      po2=rgast*rmw(ko2)*spd(ko2)*1.0e-6
      hwx=po2/(po2+hwkt/hwkb)
      hww=12.*(hwx*hwka*po2/(1.+hwkz*po2)+hwkb*po2*(1.-hwx))
      zoh=0.1*spd(koh)*avogad*rmw(koh)*sqrt(rgast*rmw(koh)*r2pi)
      dnuc=dnuc-2.*dtc*rrhopn*(hww+xmc*zoh)
c +++
c +++ 5th, create the soot-forming species (surovikin eq. 10, p. 97)
c +++
      davgpr=0.5*(drncr+dnuc)
      xnrnpr=spd(kfict)*rrnmin
      surkpr=pi*davgpr**2*xnrnpr*xmpn*sqrtk
      dnprdt=surv0 + surf*xnrnpr - surg0*xnrnpr**2 - surkpr
      xnrnpr=xnrnpr+dtc*dnprdt
c +++
c +++ 6th, create soot particle nuclei.  update particle number
c +++ density xmpn, and in turn spd(i4,kpn)
c +++
      if(con0.eq.0.0) con0=spd(kfuel)
      vovind=vovcof*(spd(kfuel)/con0)
      dnprdt=max(dnprdt,0.0)
c +++
c +++ 7th, particles destroyed by particle + particle collisions, and
c +++ calculate growth rate of particle nuclei (surovikin eq.9, p.97)
c +++ (radical nucleus has reached the critical diameter)
c +++ m3* is the reduced mass and h0 is the rate for these collisions
c +++
      fm3star=pi*twlvth*rhopn*avogad*dnuc**3
      h0=0.5*sqrt(8.*pi*rgast/fm3star)*dnuc**2
      delpn=vovind*dnprdt-h0*xmpn**2
      xmpn=xmpn+dtc*delpn
      dspdrn=vovind*dnprdt*pio6*rhopn*drncr**3
      ddavg=(third/(xmpn*davg**2))
     &     *((frcarb*dspdhcp+dspdrn)*dtc/(rhopn*pio6)-delpn*dtc*davg**3)
c ddavg=max(ddavg,0.0)
      dnuc=dnuc+ddavg*2.0
在此语句下：编写write(95,*) t,dnuc 可以输出dnuc,但是t为0。
                如果在主程序中输出这两个变量时，dnuc,t均为0。
请问这个该如何是好呢？该怎么改呢？