! module tetimodule use params,only: nlev,mmx implicit none logical,private :: debug=.false. contains !----------------------------------------------------------------------- subroutine teti(iday,ut,dt) use flds_ionatm,only: xne,qse,te,ti,xio2p,xinop,xiop use flds_atmos,only: xno2,xnn2,xno,xnh,xnhe,tn,rho use flds_ionzrt,only: qop,qnop,qo2p,qn2p,qnp use flds_heat,only: xlen,qjoul,qn implicit none ! ! Args: integer,intent(in) :: iday real,intent(in) :: ut,dt ! ! Local: integer :: k,m1 real :: r,r2,r3,r4,epsi,a,b,clei ! do k=1,nlev ! was do 1 r=alog(xne(k)/(xno2(k)+xnn2(k)+0.1*xno(k)+xnh(k)+xnhe(k))) r2=r*r r3=r2*r r4=r3*r epsi=exp(-(12.75+6.941*r+1.166*r2+0.08034*r3+0.001996*r4)) qse(k)=(qop(k)+qnop(k)+qo2p(k)+qn2p(k)+qnp(k))*epsi enddo ! k=1,nlev m1 = nlev-mmx do k=1,m1 ! was do 2 te(k)=tn(k) ti(k)=tn(k) xlen(k)=1.e-20 enddo ! k=1,m1 ! ! te solver: call tesolv(ut,dt) ! do k=1,nlev ! was do 7 clei=((2.406e-7*xio2p(k)+4.8e-7*xiop(k)+2.566e-7*xinop(k))* | xne(k)/te(k)**1.5) a=clei b=(6.6e-14*xnn2(k)+5.8e-14*xno2(k)+0.21e-14*xno(k)* | sqrt(2.*tn(k)))*xiop(k)+(5.45e-14*xno2(k)+5.9e-14*xnn2(k)+ | 4.5e-14*xno(k))*xinop(k)+(5.8e-14*xnn2(k)+4.4e-14*xno(k)+ | 0.14e-14*xno2(k)*sqrt(tn(k)))*xio2p(k) ti(k)=(a*te(k)+b*tn(k)+qjoul(k)*rho(k)/1.602e-12)/(a+b) if (xlen(k) < 1.e-20) xlen(k)=1.e-20 qn(k)=qn(k)+xlen(k) enddo ! k=1,nlev end subroutine teti !----------------------------------------------------------------------- subroutine tesolv(ut,dt) use flds_modelz,only: zp use flds_ionatm,only: te use flds_atmos,only: tn implicit none ! ! Args: real,intent(in) :: ut,dt ! ! Local: real :: usoln(1,mmx),zpts(mmx) integer,parameter :: ikopf=30000, ikop=400 real :: ekcn(ikopf),tm,eps,ddt,tend integer :: iwcn(ikop),itype,nbug,index,nwkcn,niwcn,npde, | mblft(1),mbrht(1),m1,m,k real :: wkcn(ikopf) ! ! Common /syspar/ is retained from old code for pdedif solver: real :: uround integer :: min,mout common/syspar/ uround,min,mout ! uround = 7.e-15 ! /syspar/ for pdedif min = 5 ! /syspar/ for pdedif mout = 6 ! /syspar/ for pdedif tm = 0. eps=1.e-2 ddt=1.e-9 itype=2 nbug=0 index=1 nwkcn=ikopf niwcn=ikop npde=1 mblft(1)=1 mbrht(1)=2 tend = dt m1=nlev-mmx do k=1,mmx ! was do 1 m=m1+k zpts(k)=zp(m) usoln(1,k) = te(m) enddo ! k=1,mmx ! write(6,"('tesolv before pdedif: usoln=',/,(6e12.4))") usoln ! ! pde solver for te (pdedif.F): if (debug) write(6,"('tesolv calling pdedif..')") call pdedif(func,bnbdy,npde,mmx,tend,eps,zpts,itype,mblft, | mbrht,nbug,tm,ddt,index,usoln,wkcn,iwcn,nwkcn,niwcn) if (debug) write(6,"('tesolv after pdedif..')") ! do k=1,mmx ! was do 2 m=m1+k te(m)=amax1(usoln(1,k),tn(m)) enddo ! k=1,mmx ! diffs starting at k=57 (nstep=19 in solver when nbug=2) ! write(6,"('tesolv after pdedif: te=',/,(6e12.4))") te end subroutine tesolv !----------------------------------------------------------------------- subroutine func(npde,nzp2,zn,tm,um,uzn,an,bn,cn) use flds_ionatm,only: xne,ti,xio2p,xinop,xiop,qse use flds_atmos,only: xnn2,xno2,xno,sht,rho use flds_modelz,only: zpht use flds_atmos,only: tn use flds_heat,only: xlen use cons,only: boltz,erg ! ! User-provided subroutine for pdedif solver (see pdedif.F). ! implicit none ! ! Args: integer,intent(in) :: npde,nzp2 real,intent(in) :: zn(nzp2),um(npde,nzp2),uzn(npde,nzp2),tm real,intent(out) :: an(npde,nzp2),bn(npde,nzp2),cn(npde,nzp2) ! ! Local: integer :: m,m1,k,kk,mmxm1,mmx1,mmx2 real,dimension(mmx) :: conde,dconde,dsht,qseh real :: quench,a6300,tee,tnn,tii,sqte,xene,qd1,qd2,qd3, | phi100,vhf,clight,phihf,xmt,vp,ve,dzit,dziu,tetn,teti, | xleoe,xleh,xlehe,xleo2e,xlen2r,xleo2r,xlen2e,aza,clei, | xleo1d,xlen2v,dca,xleo3p,xleo2v,ce,rdglo,xlei,a1d ! quench = 2.3e-11 a1d=6.81e-3 ! local (not the one in flds_ratecoef) a6300=5.15e-3 m1=nlev-mmx do k=1,mmx ! was do 22 m=m1+k kk=k+1 tee=amax1(um(1,kk),100.) sqte=sqrt(tee) xene=xne(m) qd1=2.82e-17*sqte-3.41e-21*tee*sqte qd2=2.2e-16+7.92e-18*sqte qd3=3.4e-16 conde(k)=(7.5e+5*tee*tee*sqte/(1.+3.22e+4*tee*tee/xene* | (qd1*xnn2(m)+qd2*xno2(m)+qd3*xno(m))))*1.602e-12 enddo ! k=1,mmx ! ! Alternate values for phi100 (commented out in old code): 12., 0.12 phi100 = 0. vhf = 5.0e+7 clight=3.e+10 do k=1,mmx ! was do 76 m=m1+k kk=k+1 tee=amax1(um(1,kk),100.) sqte=sqrt(tee) xmt=xno(m)+xno2(m)+xnn2(m) vp=5.6e+4*sqrt(xne(m)) ve=5.4e-10*xmt*sqte+(34.+4.18*alog(tee**3/xne(m)))*xne(m)/ | (tee*sqte) phihf=phi100*(100.e+5/(zpht(m)*1.e+5))**2 qseh(k)=phihf/clight*(vp*vp/(vhf*vhf+ve*ve))*ve enddo ! k=1,mmx dzit=0.5 dziu=0.25 mmxm1=mmx-1 do k=2,mmxm1 ! was do 23 m=m1+k dconde(k)=(conde(k+1)-conde(k-1))/dzit dsht(k)=(sht(m+1)-sht(m-1))/dzit enddo ! k=2,mmxm1 dconde(1)=(conde(2)-conde(1))/dziu dsht(1)=(sht(m1+1)-sht(m1))/dziu dconde(mmx)=(conde(mmx)-conde(mmx-1))/dziu dsht(mmx)=(sht(nlev)-sht(nlev-1))/dziu mmx1=mmx+1 mmx2=mmx+2 do kk=2,mmx1 ! was do 2 k=kk-1 m=m1+k tee=amax1(um(1,kk),100.) tnn=tn(m) tii=ti(m) xene=xne(m) tetn=tee-tnn teti=tee-tii sqte=sqrt(tee) xlen2r=2.9e-14*xene*xnn2(m)*tetn/sqte xleo2r=6.9e-14*xene*xno2(m)*tetn/sqte xlen2e=1.77e-19*xene*xnn2(m)*(1.-1.21e-4*tee)*tee*tetn xleo2e=1.21e-18*xene*xno2(m)*(1.+3.6e-2*sqte)*sqte*tetn xleoe=7.9e-19*xene*xno(m)*(1.+5.7e-4*tee)*sqte*tetn ! xleh=9.63e-16*xene*xnh(m)*(1.-1.35e-4*tee)*sqte*tetn ! xlehe=2.46e-17*xene*xnhe(m)*sqte*tetn xleh=0. xlehe=0. if (tee >= 1000.) then if (tee <= 2000.) then aza=2.e-7*exp(-4605.2/tee) else aza=2.53e-6*sqte*exp(-17620./tee) endif else aza=5.71e-8*exp(-3352.6/tee) endif dca=exp(3200.*(1./tee-1./tnn)) xlen2v=1.3e-4*xnn2(m)*xene*(1.-dca)*aza xleo2v=3.125e-21*xene*xno2(m)*tetn*tee*tee xleo3p=3.4e-12*tetn*(1.-7.e-5*tee)/tnn*xno(m)*xene* | (150./tee+0.40) xleo1d=1.07e-10*xene*xno(m)*sqte*exp(-2.27e+4/tee)* | (0.406+0.357e-4*tee-(0.333+0.183e-4*tee)*exp(-1.37e+4/tee)- | (0.456+0.174e-4*tee)*exp(-2.97e+4/tee)) clei=((2.406e-7*xio2p(m)+4.8e-7*xiop(m)+2.566e-7*xinop(m))* | xene/tee**1.5) xlei=clei*teti rdglo=(xleo1d/1.96)*a6300/(a1d+(1.+quench*xnn2(m)/a1d)) xlen(m)=(xlen2r+xleo2r+xlen2e+xleo2e+xleoe+xlen2v+xleo2v+xleo3p+ | xleh+xlehe+(xleo1d-1.96*rdglo)+xlei)*erg/rho(m) ce=1.5*boltz*xene an(1,kk)=conde(k)/(sht(m)*sht(m)*ce) bn(1,kk)=(dconde(k)/(sht(m)*sht(m))-conde(k)*dsht(k)/ | (sht(m)**3))/ce cn(1,kk)=(qse(m)-(xlen2r+xleo2r+xlen2e+xleo2e+xleoe+xleh+xlehe | +xlen2v+xleo2v+xleo3p+xleo1d+xlei))*1.602e-12/ce | +qseh(k)/ce enddo ! kk=2,mmx1 (was do 2) an(1,1)=an(1,2) bn(1,1)=bn(1,2) cn(1,1)=cn(1,2) an(1,mmx2)=an(1,mmx1) bn(1,mmx2)=bn(1,mmx1) cn(1,mmx2)=cn(1,mmx1) ! write(6,"('teti func: an=',/,(6e12.4))") an ! write(6,"('teti func: bn=',/,(6e12.4))") bn ! write(6,"('teti func: cn=',/,(6e12.4))") cn end subroutine func !----------------------------------------------------------------------- subroutine bnbdy(npde,tm,zl,zr,al,ar,bl,br) use flds_atmos,only: tn,xno2,xnn2,xno,sht use flds_ionatm,only: xne use input,only: f107a use cons,only: erg ! ! User-provided subroutine for pdedif solver (see pdedif.F). ! implicit none ! ! Args: integer,intent(in) :: npde real,intent(in) :: tm,zl,zr real,intent(in) :: al(npde),ar(npde) real,intent(out) :: bl(npde),br(npde) ! ! Local: integer :: m,m1 real :: qtophe,conde,xene,tee,qd1,qd2,qd3,sqte ! integer :: ncalls=0 ! ! ncalls = ncalls+1 m1 = nlev-mmx m = m1+1 bl(1) = tn(m) ! ! qtophe is a tunable parameter: ! qtophe=(-4.47E+7*F107A*0.125-7.0E+8) ! Roble tune ! qtophe=(-4.47E+7*F107A*0.125-7.0E+8)*3.0 ! Roble tune ! ! qtophe=-4.47e+7*f107a*1.2 ! use this value for solar minimum qtophe=-4.47e+7*f107a*0.7 ! use this value for solar medium (old model) ! qtophe=-4.47e+7*f107a*0.1 ! use this value for solar maximum ! tee=amax1(ar(1),100.) xene=xne(nlev) sqte=sqrt(tee) qd1=2.82e-17*sqte-3.41e-21*tee*sqte qd2=2.2e-16+7.92e-18*sqte qd3=3.4e-16 conde=(7.5e+5*tee*tee*sqte/(1.+3.22e+4*tee*tee/xene* | (qd1*xnn2(nlev)+qd2*xno2(nlev)+qd3*xno(nlev))))*1.602e-12 br(1)=-(qtophe*erg*sht(nlev))/conde ! write(6,"('teti bnbdy: ncalls=',i4,' qtophe,erg,sht,conde=', ! | 4e12.4,' br=',e12.4)") ncalls,qtophe,erg,sht(nlev),conde,br(1) ! write(6,"('teti bnbdy: ncalls=',i4,' bl=',e12.4,' br=',e12.4)") ! | ncalls,bl(1),br(1) end subroutine bnbdy !----------------------------------------------------------------------- end module tetimodule