! module denmodule use params,only: nlev implicit none real :: alam,elec ! formerly in /qiondc/ real :: alamm(nlev) ! formerly in /lamtrsf/ real :: ar(54) ! formerly in /qiondc/ real :: phoxo2,phoxno,phoxi,phoxt ! formerly in /piontr/ real,dimension(nlev) :: qhoxi,qhoxo2,qhoxno,qhoxt ! formerly /ionhtr/ logical,private :: debug=.false. contains !----------------------------------------------------------------------- subroutine denmod(iday,ut,dt) use flds_atmos,only: xno,xno2,xnn2,xnco2,xnh2o,xno2,xno3,xnh,xnno, | xnno2,tn,rho,xn2d,xnh2,xn4s use flds_ionzrt,only: qpro,qop2p,qop2d,qop4s use flds_modelz,only: zpht use flds_col,only: txno2 use flds_ratecoef,only: alp1,alp2,alp3,rk1,rk2,rk3,rk4,rk5,rk6, | rk7,rk8,rk9,rk10,rk11,rk12,rk13,rk14,rk15,rk16,rk17,rk18,rk19, | rk20,rk21,rk22,rk23,rk24,rk25,rk26,rk27,rk28 use flds_ionzrt,only: qnp,qnop,qo2p,qn2p,qop use flds_ionatm,only: xne,xnpi,xnni,xneee,xio2p,xinop,xar,te, | parcd,sped,shal,xmup,xmun,sigp,sign,xmolpi,xmolni,xiop,xinp, | xihp,xin2p,xiop2p,xiop2d,xiop4s,alamxx,alamxy,te,ti use opdif_solve,only: opdifsv use flds_heat,only: qjoul,qn use input,only: istep,nstep ! ! Args: integer,intent(in) :: iday real,intent(in) :: ut,dt ! ! Local: integer :: k,j,m integer,parameter :: kiy=36 real,parameter :: | bgaus = 0.4, | xme = 5.486e-4, | qchr = 1.602e-19, ! | colfac = 1.0 | colfac = 1.5 ! old model real :: xnt,o2s,hno3,qio2p,xposi,xnegi,sum,sumn,xmolp,xmoln, | wnegi,wposi,bgaus2,bg,ctts,we,vinn,tw,te12,ven,vinp,vei, | cn,ckoe,ckop,ckon,qchr2,vp,vn,ve,ckpe,ckpp,ckpn,ckhe,ckhp, | ckhn,airm,alam1,alam2,aa,bb,cc,dd,ee,ff,gg,hh,a0,a1,a2,a3,a4, | dum,vopo,vo2po2,vopn2,vopo2,vnopo,vnopn2,vnopo2,vo2po,vo2pn2, | tr,wt1,wt2,wt3,winop,wi,wiop,wio2p,tlamxx,rw,rzm,rz,emv, | efield,ajoul,ax,tlam real,dimension(nlev) :: pop4s,xlop4s real :: xmass(40) = | (/ 0., 0., 0., 0., 0., 30., 48., 66., 84., 58., | 76., 94., 74., 92.,110., 32., 64., 80., 50., 36., | 19., 37., 55., 73., 91.,109.,127.,145., 16., 32., | 48., 64., 46., 62., 60., 76., 62., 60., 44., 34. /) ! ! rlam was declared complex in old denmod, but dummy arg in ! sub quart was implicitly real (see solvers.F) complex :: rlam ! real :: rlam real :: h,rl,cdelta,g,root,delta ! output from quart real :: acof(5) ! coeff input to quart ! ! write(6,"('denmod: iday=',i4,' ut=',f8.2)") iday,ut qchr2 = qchr*qchr bgaus2 = bgaus*bgaus bg = bgaus*1.e-4 ctts=9.57946e+3*bgaus we = ctts/xme alamm = 0. ! whole array init (was in old glbmean.F) xar = 1.e-20 ! whole array init (not in old code) do k=1,kiy xnt=xno(k)+xno2(k)+xnn2(k) o2s = 1.e-20 hno3 = 1.e-20 qio2p = qo2p(k)+qn2p(k)+qop(k)+qnp(k) elec = xne(k) ! elec is module data above alam = alamm(k) ! alam and alamm are module data above call ioncomp(iday,ut,k,zpht(k),tn(k),xnt,xnn2(k),xnco2(k), | xnh2o(k),xno2(k),xno3(k),xno(k),xnno(k),xnno2(k),o2s,hno3, | qio2p,qnop(K),txno2(k),alp1(k),alp2(k)) ! ! qhoxi,qhoxo2,qhoxno,qhoxt are denmod module data above, ! phoxt, phoxo2, phoxno are denmod module data above, ! qpro is in module flds_ionzrt, and was defined by solheat ! qhoxi(k) = phoxt*qpro(k) qhoxo2(k) = phoxo2 qhoxno(k) = phoxno qhoxt(k) = phoxt ! ! xnpi, etc are in module flds_ionatm: xnpi(k) = ar(4) xnni(k) = ar(5) xneee(k) = ar(2) xio2p(k) = ar(16) xinop(k) = ar(6) xne(k) = ar(2) ! ! Note: this new alam will be reused by sub ioncomp in next k-iteration: alam = ar(3) alamm(k) = alam ! ! xposi,xnegi are local: xposi = ar(4) xnegi = ar(5) ! ! ar(54) is denmod module data, xar(60) is in module flds_ionatm: xar(1:42,k) = ar(1:42) ! For Dan Marsh's d-region, comment out xar(43:54,k) statement. ! (see re-init of xar=1.e-20 below) xar(43:54,k) = ar(43:54) ! ! xmass is local (was initialized w/ array constructor above): sum = 0. ! local do j=6,28 ! was do 46 sum = sum+xmass(j)*ar(j) enddo ! j=6,28 sumn = 0. do j=29,40 ! was do 47 sumn = sumn+xmass(j)*ar(j) enddo ! j=29,40 xmolp = sum/ar(4) xmoln = sumn/ar(5) wnegi=ctts/xmoln wposi=ctts/xmolp vinn=2.6e-09*xnt/sqrt(xmoln) tw=te(k) ! te is in flds_ionatm te12=sqrt(tw) ven=2.33e-11*xnn2(k)*tw*(1.-1.21e-4*tw)+1.82e-10*xno2(k)*te12* | (1.+ 3.6e-2*te12)+8.2e-10*xno(k)*te12 vinp=2.6e-09*xnt/sqrt(xmolp) vei=34.+4.18*alog10(tw**3/xne(k))/tw**1.5 cn=1./(bg*qchr) ckoe=cn*we/(vei+ven) ckop=cn*wposi/(vinp+vei) ckon=cn*wnegi/(vinn+vei) parcd(k)=xne(k)*1.e+6*qchr2*(ckoe+(1.+alam)*ckop+alam*ckon) vp=vinp vn=vinn ve=ven+vei ckpe=cn*we*ve/(ve*ve+we*we) ckpp=cn*wposi*vp/(vp*vp+wposi*wposi) ckpn=cn*wnegi*vn/(vn*vn+wnegi*wnegi) sped(k)=xne(k)*1.e+6*qchr2*(ckpe+(1.+alam)*ckpp+alam*ckpn) ckhe=cn*we*we/(ve*ve+we*we) ckhp=cn*wposi*wposi/(vp*vp+wposi*wposi) ckhn=cn*wnegi*wnegi/(vn*vn+wnegi*wnegi) shal(k)=xne(k)*1.e+6*qchr2*(ckhe-(1.+alam)*ckhp-alam*ckhn) airm = 28.9 xmup(k)=35./sqrt(11.5*airm*xmolp/(airm+xmolp)) xmun(k)=35./sqrt(11.5*airm*xmoln/(airm+xmoln)) sigp(k)=1.602*xnpi(k)*xmup(k)*2.7e+2/xnt sign(k)=1.602*xnni(k)*xmun(k)*2.7e+2/xnt xmolpi(k)=xmolp xmolni(k)=xmoln xiop(k) =1.e-20 xinp(k) =1.e-20 xihp(k) =1.e-20 xin2p(k) =1.e-20 xiop2p(k)=1.e-20 xiop2d(k)=1.e-20 xiop4s(k)=1.e-20 ! alam1 is local to denmod (note there is also a different alam1 in ioncomp) alam1=sped(k)*bgaus2/rho(k)*1.e-11 alam2=shal(k)*bgaus2/rho(k)*1.e-11 alamxx(k)=alam1 alamxy(k)=alam2 enddo ! k=1,kiy (with ioncomp call) ! call fprint8('denmod before pdedif:',kiy, ! | sped , shal , alamxx , alamxy , dum,dum,dum,dum, ! | 'sped','shal','alamxx','alamxy',' ',' ',' ',' ') do k=1,nlev ! was do 55 xiop2p(k)=qop2p(k)/((rk16(k)+rk17(k))*xnn2(k)+rk18(k)*xno(k) | +(rk19(k)+rk20(k))*xne(k)+rk23(k)*xno2(k)+rk21(k) | +rk22(k)) xiop2d(k)=(qop2d(k)+rk20(k)*xiop2p(k)*xne(k)+rk22(k)*xiop2p(k)) | /(rk24(k)*xnn2(k)+rk25(k)*xno(k)+rk26(k)*xne(k) | +rk27(k)*xno2(k)+rk28(k)) xinp(k)=(qnp(k)+rk10(k)*xiop4s(k)*xn2d(k)+rk17(k)*xiop2p(k)* | xnn2(k))/((rk6(k)+rk7(k))*xno2(k)+rk8(k)*xno(k)) xihp(k)=rk12(k)*xnh(k)*xiop4s(k)/(rk11(k)*xno(k)) enddo ! k=1,nlev ! 2/16/05: no diffs ! call fprint8('denmod xiop4s, et.al.:',kmx, ! | xiop4s , xiop2p , xiop2d , xinp , xihp ,dum,dum,dum, ! | 'xiop4s','xiop2p','xiop2d','xinp','xihp','','','') ! call fprint8('denmod xiop4s, et.al.:',nlev, ! | xiop4s , xiop2p , xiop2d , xinp , xihp , xnh , xno, qop2p, ! | 'xiop4s','xiop2p','xiop2d','xinp','xihp','xnh','xno','qop2p') ! ! 8/22/05 btf: added zero init to pop4s,xlop4s: pop4s(:) = 0. xlop4s(:) = 0. do m=1,kiy ! was do 5 pop4s(m)=qop4s(m)+(rk18(m)*xno(m)+rk19(m)*xne(m)+rk21(m))* | xiop2p(m)+(rk25(m)*xno(m)+rk26(m)*xne(m)+rk28(m))* | xiop2d(m)+rk8(m)*xno(m)*xinp(m)+rk11(m)*xno(m)*xihp(m) xlop4s(m)=rk1(m)*xno2(m)+rk2(m)*xnn2(m)+rk10(m)*xn2d(m)+ | rk12(m)*xnh(m)+rk13(m)*xnco2(m)+rk14(m)*xnh2(m)+ | rk15(m)*xnh2o(m) xiop4s(m)=pop4s(m)/xlop4s(m) xiop(m)=xiop4s(m) enddo ! m=1,kiy ! call fprint8('denmod before opdifsv:',nlev, ! | pop4s , qop4s , xno , xne , xiop2p , xiop2d , xinp , xihp, ! | 'pop4s','qop4s','xno','xne','xiop2p','xiop2d','xinp','xihp') ! call fprint8('denmod before opdifsv:',nlev, ! | xlop4s , xno2 , xnn2 , xn2d , xnh , xnco2 , xnh2 , xnh2o, ! | 'xlop4s','xno2','xnn2','xn2d','xnh','xnco2','xnh2','xnh2o') ! ! O+ diffusion solver (opdifsv.F) calls pdedif (pdedif.F). if (debug) write(6,"('denmod calling opdifsv..')") call opdifsv(ut,dt,colfac) if (debug) write(6,"('denmod after opdifsv..')") ! do k=kiy,nlev aa=qnop(k)+rk2(k)*xiop4s(k)*xnn2(k)+rk7(k)*xno2(k)*xinp(k) bb=qo2p(k)+rk1(k)*xno2(k)*xiop4s(k)+rk6(k)*xno2(k)*xinp(k) | +rk13(k)*xnco2(k)*xiop4s(k)+rk23(k)*xno2(k)*xiop2p(k) | +rk27(k)*xno2(k)*xiop2d(k) cc=rk4(k)*xn4s(k)+rk5(k)*xnno(k) dd=qn2p(k)+rk16(k)*xiop2p(k)*xnn2(k)+rk24(k)*xiop2d(k)*xnn2(k) ee=rk3(k)*xno(k)+rk9(k)*xno2(k) ff=xiop4s(k) gg=xinp(k) hh=rk9(k)*xno2(k) a4=alp1(k)*alp2(k)*alp3(k) a3=alp1(k)*(alp2(k)*ee+alp3(k)*cc)-alp1(k)*alp2(k)*alp3(k)* | (ff+gg) a2=alp1(k)*ee*cc-alp1(k)*(alp2(k)*ee+alp3(k)*cc)*(ff+gg)- | alp1(k)*alp2(k)*dd-alp1(k)*alp3(k)*bb-alp2(k)*alp3(k)*aa a1=-alp1(k)*(ee*(cc*(ff+gg)+bb)+dd*(cc+hh))-alp2(k)*(ee*(aa+dd)- | hh*dd)-alp3(k)*cc*(aa+bb) a0=-ee*cc*(aa+bb+dd) acof(1)=a0 acof(2)=a1 acof(3)=a2 acof(4)=a3 acof(5)=a4 ! ! Quartic solver (in this module): ! ! write(6,"('denmod before quart: acof=',5e12.4)") acof if (debug) write(6,"('denmod calling quart: k=',i4, | ' nlev=',i4)") k,nlev call quart(acof,root,delta,g,h,rl,cdelta,rlam) ! ! write(6,"('denmod after quart: k=',i3,' root,delta,g,h=', ! | 4e12.4)") k,root,delta,g,h ! write(6,"(' k=',i3,' rl,cdelta,rlam=', ! | 4e12.4)") k,rl,cdelta,rlam ! 4e12.4 because rlam is complex ! write(6,"('denmod after quart: k=',i3,' root=',e12.4)") k,root ! xne(k) = root xiop4s(k)=ff xiop(k)=xiop4s(k) xin2p(k)=dd/(ee+alp3(k)*xne(k)) xinp(k)=gg xio2p(k)=(bb+hh*xin2p(k))/(cc+alp2(k)*xne(k)) xinop(k)=(aa+(ee-hh)*xin2p(k)+cc*(bb+hh*xin2p(k))/ | (cc+alp2(k)*xne(k)))/(alp1(k)*xne(k)) xnt=xno(k)+xno2(k)+xnn2(k) xmolp=(xin2p(k)*28.+xio2p(k)*32.+xiop4s(k)*16.+xinop(k)*30.+ | xinp(k)*14.)/xne(k) wposi=ctts/xmolp wnegi=0. alam=0. qhoxi(k)=0. qhoxo2(k)=0. qhoxno(k)=0. xneee(k)=xne(k) xnpi(k)=xiop4s(k)+xio2p(k)+xinop(k)+xin2p(k)+xinp(k)+xihp(k)+ | xiop2p(k)+xiop2d(k) xnni(k)=1.e-20 we=1.7588028e+7*bgaus tr=(ti(k)+tn(k))*0.5 vopo =3.67e-11*xno(k)*sqrt(tr)*(1.-0.064*alog10(tr))**2*colfac vo2po2 =2.59e-11*xno2(k)*sqrt(tr)*(1.-0.073*alog10(tr))**2 vopn2 =6.82e-10*xnn2(k) vopo2 =6.64e-10*xno2(k) vnopo =2.44e-10*xno(k) vnopn2 =4.34e-10*xnn2(k) vnopo2 =4.27e-10*xno2(k) vo2po =2.31e-10*xno(k) vo2pn2 =4.13e-10*xnn2(k) wi=9.6489e+3*bgaus wiop=wi/16. wio2p=wi/32. winop=wi/30. wt1=(vopo+vopn2+vopo2)/wiop wt2=(vo2po2+vo2po+vo2pn2)/wio2p wt3=(vnopo+vnopo2+vnopn2)/winop tw=te(k) te12=sqrt(tw) ven=2.33e-11*xnn2(k)*tw*(1.-1.21e-4*tw)+1.82e-10*xno2(k)*te12* | (1.+ 3.6e-2*te12)+8.2e-10*xno(k)*te12 sped(k)=1.e+10*qchr/bgaus*(xiop4s(k)*wt1/(1.+wt1*wt1)+xio2p(k) | *wt2/(1.+wt2*wt2)+xinop(k)*wt3/(1.+wt3*wt3)+xne(k)*ven/we/ | (1.+(ven/we)**2)) shal(k)=1.e+10*qchr/bgaus*(xne(k)/(1.+(ven/we)**2)-xiop4s(k)/ | (1.+wt1*wt1)-xio2p(k)/(1.+wt2*wt2)-xinop(k)/(1.+wt3*wt3)) vinp=vopo+vopn2+vopo2+vo2po+vo2po2+vo2pn2+vnopo+vnopo2+ | vnopn2 vei=34.+4.18*alog10(tw**3/xne(k))/tw**1.5 vinn=0. cn=1./(bg*qchr) ckoe=cn*we/(vei+ven) ckop=cn*wposi/(vinp+vei) ckon=cn*wnegi/(vinn+vei) parcd(k)=xne(k)*1.e+6*qchr2*(ckoe+(1.+alam)*ckop+alam*ckon) qhoxt(k)=0. xmup(k)=0. xmun(k)=0. sigp(k)=0. sign(k)=0. xmolpi(k)=0. xmolni(k)=0. alam1=sped(k)*bgaus2/rho(k)*1.e-11 alam2=shal(k)*bgaus2/rho(k)*1.e-11 alamxx(k)=alam1 alamxy(k)=alam2 xar(2:54,k) = 1.e-20 enddo ! k=kiy,nlev (includes quartic) (was do 44) ! tlamxx = 0. do k=2,nlev rz=rho(k) rzm=rho(k-1) rw=0.5*(zpht(k)-zpht(k-1))*1.e+5 tlamxx=tlamxx+(alamxx(k)*rz+alamxx(k-1)*rzm)*rw enddo ! k=2,nlev ax=5.163e+11 tlam=tlamxx*ax ! Alternative values for ajoul (commented out in old code): 0. and 1.0e11 ajoul=3.0e+10 emv=sqrt(ajoul*(bgaus*1.e-3)**2/tlam) efield=emv/(bgaus*1.e-3) do k=1,nlev ! was do 82 qjoul(k)=alamxx(k)*efield*efield qn(k)=qn(k)+qjoul(k) ! qn was init to nighttime ionization in solheat enddo ! k=1,nlev ! if (istep==nstep) then ! call fprint8('Ionosphere fields (denmod):',nlev, ! | sped , shal , alamxx , alamxy , qjoul , qn, dum,dum, ! | 'sped','shal','alamxx','alamxy','qjoul','qn',' ',' ') ! ! Diffs w/ old code in xne after 1st timestep: ! call fprint8('Ionosphere fields (denmod):',nlev, ! | xne , xiop4s , xin2p , xinp , xio2p , xinop, xnpi,dum, ! | 'xne','xiop4s','xin2p','xinp','xio2p','xinop','xnpi',' ') ! endif end subroutine denmod !----------------------------------------------------------------------- subroutine ioncomp(iday,ut,kht,h,te,am,an2,co2,h2o,o2,o3,o1,ano1, | ano2,o2s,hno3,qio2p,qinop,o2col,alph1,alph3) ! ! Args: integer,intent(in) :: iday,kht real,intent(in) :: ut,h,te,am,an2,co2,h2o,o2,o3,o1,ano1,ano2,o2s, | hno3,qio2p,qinop,o2col,alph1,alph3 ! ! Local: logical,parameter :: iswtc = .true. ! print switch integer :: iabc,i real,parameter :: | r6 =1.e-9 ,r7 =1.e-9 ,r8 =1.e-9 ,r14=1.e-9 ,r15=1.e-9 , | r16=1.e-9 ,r22=1.e-9 ,r23=1.e-9 ,r24=1.e-9 ,r25=7.e-11 , | r41=1.5e-9 ,r42=1.e-9 ,r43=1.e-9 ,r47=2.8e-28 ,r44=1.e-9 , | r45=2.e-10 ,r46=1.4e-9 ,r48=3.e-10 ,rn2=1.e-31 ,rn3=1.5e-10 , | rn4=1.5e-10 ,rn5=2.5e-10 ,rn6=1.1e-10 ,rn7=6.e-10 ,rn8=5.5e-10 , | rn9=1.1e-11 ,rn10=2.5e-12,rn11=1.2e-10,rn12=4.3e-10,rn13=2.5e-10, | rn14=4.8e-11,rn15=4.e-10 ,rn16=1.4e-10,rn17=1.9e-10,rn18=8.e-10 , | rn19=4.e-10 ,rn20=2.e-10 ,rn21=1.5e-11,rn22=1.3e-10,rn23=8.e-10 , | rn26=1.e-9 ,rn29=1.2e-9 ,rn30=1.4e-9 ,rn31=3.e-10 ,rn32=2.e-10 , | rn33=5.8e-10,rn34=3.e-28 ,rn35=4.7e-29,rn37=1.e-9 ,rn38=1.6e-9 , | rn39=8.e-10 ,rn40=3.e-9 ,rn41=2.8e-9 ,r49=4.4e-10 real,parameter :: | s1=0., s2=0., s3 =0., s4 =0., s5 =0., s6=0., s7=0., s8=0.,s9=0., | s10=0.,s11=0., s12=0. ! real :: alam1,tnwew,alph2,rec1,rec2,rec3,alph4,alph5, | ao2,ao4,ao5,ao2h,ahhoh,qno,qo2,aneg,test,test1,test2,eqdfy,ab1, | aqdfy real :: r1,r2,r3,r5,r9,r10,r11,r13,r17,r18,r19,r21,r26,r27,r28, | r29,r30,r31,r36,r37,r38,r39,r32,r33,r34,r35,r4,r12,r20,r40,rn1, | rn24,rn25,rn27,rn28,rn36 real :: c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16, | c17,c18,c19,c20,c21,c22,c23,c24,c25,c26,c27,c28,c29,c30,c31, | c32,c33,c34,c35,c36,c37,c38,c39,c40,c41,c42,c43,c44 real :: ano,anon1,anoc1,anoh1,anohn,anohc,anoh2,anoh2n,anoh2c, | anoh3 real :: w1,w2,w3,w4,w5,w6,w7,w8 real :: sum1,sum2,sum3,sum4,sum5,sumt,sumn1,sumn2 real :: ga1,ga2,ga3,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12 real :: bl10,bl7,bl12,bl8,bl2,bla2,blb11,bl11,bla11,bl9,blb3,blc3, | bld3,bl3,bla3,bl4,bla4,blb5,blc5,bl5,bla5,blb6,bl6,blc6,bla6 real :: coe1,coe2,coe3 real :: al1,al2,al3,al4,al5,al6,al7,al8,al9,al10,al11,al12 real :: bo2,bo1,bo3,bco3,bno2,bno3,bo4,bco4,bo2no,bo2n2,bo1n2, | bo2h2o real :: aneg1,qit,rk5,alf1,alf2,alf3,xl1,xl2,xl3,gam1,p1,p2,bet1, | xl4,xl5,xl6,gam2,p3,p4,bet2,bet3 ! ! Print inputs: ! write(6,"('ioncomp: kht=',i3,' inputs 1-6 =',6e12.4)") ! | kht,h,te,am,an2,co2,h2o ! write(6,"('ioncomp: kht=',i3,' inputs 7-12 =',6e12.4)") ! | kht,o2,o3,o1,ano1,ano2,o2s ! write(6,"('ioncomp: kht=',i3,' inputs 13-18=',6e12.4)") ! | kht,hno3,qio2p,qinop,o2col,alph1,alph3 r1 = 2.e-31*(300./te)**4.4 r2 = 1.5e6/(te**5.4*exp(2450./te)) r3 = 7.e-30*(300./te)**3. r5 = 1.8e-28*(308./te)**4.7 r9 = 2.e-31*(300./te)**4.4 r10 = 1.5e6/(te**5.4*exp(2150./te)) r11 = 7.e-30*(300./te)**3. r13 = 1.e-27*(308./te)**4.7 r17 = 2.e-31*(300./te)**4.4 r18 = 1.5e6/(te**5.4*exp(1800./te)) r19 = 7.e-30*(300./te)**3. r21 = 1.e-27*(308./te)**4.7 r26 = 2.4e-27*(300./te)**4. r27 = 1.36e11/(te**5.*exp(8360./te)) r28 = 9.e-28*(300./te)**4. r29 = 6.9e11/(te**5.*exp(7670./te)) r30 = 9.e-28*(300./te)**4. r31 = 2.1e11/(te**5.*exp(6540./te)) r36 = 2.3e-27*(300./te)**4. r37 = 1.95e11/(te**5.*exp(11000./te)) r38 = 3.4e-27*(300./te)**4. r39 = 9.6e11/(te**5.*exp(17100./te)) r32 = 9.e-28*(300./te)**4. r33 = 1.26e11/(te**5.*exp(5830/te)) r34 = 9.e-28*(300./te)**4. r35 = 2.3e10/(te**5.*exp(5000./te)) r4 = 3.11e4/(te**4.*exp(4590./te)) r12 = 3.11e4/(te**4.*exp(4025./te)) r20 = 3.11e4/(te**4.*exp(3335./te)) r40 = 2.6e-30*(300./te)**3.2 rn1 = 1.4e-29*(300./te)*exp(-600./te) rn24 = 3.1e-31*(300./te)**2.5 rn25 = 6.3e-32*(300./te)**4.2 rn27 = 4.7e-31*(300./te)**1.9 rn28 = 1.8e-31*(300./te)**2.5 rn36 = 9.e-12*(te/300.)**1.5 ! alam = 0. sumt = 0. tnwew = te alph2 = 3.e-6 alph4 = 7.0e-8+2.e-25*(an2+o2+o1) qno = qinop qo2 = qio2p ! iabc = 0 100 continue ! was 30 alam1 = alam*alph4 rec1 = (alph1+alam1)*elec rec2 = (alph2+alam1)*elec rec3 = (alph3+alam1)*elec c19 = r49*ano1+(r40*o2+r47*h2o)*am+rec3 c20 = c19/(r48*o1) c21 = r41*h2o+r42*o3+r48*o1+rec2 c22 = r40*o2*am/c21 ao2 = qo2/(r48*o1*(c20-c22)) if (ao2 > 1.e7) ao2=1.e7 ao4 = c22*ao2 ao5 = r42*o3*ao4/(r43*h2o+rec2) c23 = h2o*(r41*ao4+r43*ao5+r47*am*ao2) ao2h = c23/((r44+r45)*h2o+rec2) ahhoh = r44*h2o*ao2h/(r46*h2o+rec2) c1 = (r1*an2+r3*co2+r5*h2o)*am c2 = r2*am+r6*co2+r8*h2o c3 = r4*am+r7*h2o c4 = r1*an2*am/(c2+rec2) ! iabc = iabc+1 ! c5 = (r3*co2*am+r6*co2*c4)/(c3+rec2) ano = (qno+r49*ano1*ao2)/(c1+rec1-r2*am*c4-r4*am*c5) anon1 = c4*ano anoc1 = c5*ano c6 = (r9*an2+r11*co2+r13*h2o)*am c7 = r10*am+r14*co2+r16*h2o c8 = r12*am+r15*h2o c9 = r9*an2*am/(c7+rec2) c10 = (r11*co2*am+r14*co2*c9)/(c8+rec2) c11 = (r5*am*ano+r7*anoc1+r8*anon1)*h2o anoh1 = c11/(c6+rec2-(r10*c9+r12*c10)*am) anohn = c9*anoh1 anohc = c10*anoh1 ! c12 = (r17*an2+r19*co2+r21*h2o)*am c13 = r18*am+r22*co2+r24*h2o c14 = r20*am+r23*h2o c15 = r17*an2*am/(c13+rec2) c16 = (r19*am+r22*c15)*co2/(c14+rec2) c17 = (r13*am*anoh1+r15*anohc+r16*anohn)*h2o anoh2 = c17/(c12+rec2-am*(r18*c15+r20*c16)) anoh2n = c15*anoh2 anoh2c = c16*anoh2 c18 = (r21*am*anoh2+r23*anoh2c+r24*anoh2n)*h2o anoh3 = c18/(r25*h2o+rec2) ! c24 = r34*am*h2o/(r35*am+rec2) c25 = (r34*h2o+r33-r35*c24)*am+rec2 c26 = r32*am*h2o/c25 c27 = (r32*h2o+r31-r33*c26)*am+rec2 c28 = r30*am*h2o/c27 c29 = (r30*h2o+r29-r31*c28)*am+rec2 c30 = r28*am*h2o/c29 c31 = (r28*h2o+r27-r29*c30)*am+rec2 c32 = r26*am*h2o/c31 c33 = r25*h2o*anoh3 c34 = r36*am*h2o c35 = (r26*h2o+r37-r27*c32)*am+rec2 c36 = r46*h2o*ahhoh c37 = r38*am*h2o c38 = r37*am c39 = (r36*h2o+r39)*am+rec2 c40 = r45*h2o*ao2h c41 = r39*am c42 = r38*am*h2o+rec2 c43 = c42*(c35*c36+c33*c38)+c35*c37*c40 c44 = c42*(c34*c38-c35*c39)+c35*c37*c41 w2 = -c43/c44 w1 = (c40+c41*w2)/c42 w3 = (c33+c34*w2)/c35 w4 = c32*w3 w5 = c30*w4 w6 = c28*w5 w7 = c26*w6 w8 = c24*w7 sum1 = anon1+anoc1+anoh1+anohn sum2 = anohc+anoh2+anoh2n+anoh2c+anoh3 sum3 = ao4+ao5+ao2h+ahhoh sum4 = sum1+sum2+sum3 sum5 = w1+w2+w3+w4+w5+w6+w7+w8 sumt = ano+ao2+sum4+sum5 alph5 = alph4*sumt ga1 = (rn3+rn4)*o1+rn7*o3+rn20*o2s+rn23*ano2+rn41*hno3 ga2 = (rn24*o2+rn25*an2+rn34*h2o+rn35*co2)*am g1 = ga1+ga2+s1+alph5 ga3 = rn17*o1+(rn18+rn19)*o3+rn29*ano2+rn40*hno3 g2 = ga3+(rn27*o2+rn28*an2)*am+s2+alph5 g3 = rn5*o1+rn8*co2+rn10*ano1+s3+alph5 g4 = rn6*o1+rn9*ano1+rn32*ano2+rn39*hno3+s4+alph5 g5 = rn11*o3+rn38*hno3+s5+alph5 g6 = s6+alph5 g7 = rn12*co2+rn13*ano1+rn15*o1+rn30*h2o+s7+alph5 g8 = rn14*ano1+rn16*o1+rn22*o3+s8+alph5 g9 = rn21*ano1+s9+alph5 g10 = rn26*o2+s10+alph5 g11 = rn37*o2+s11+alph5 g12 = rn31*o3+rn33*co2+s12+alph5 bl10 = rn25*an2*am/g10 bl7 = (rn24*o2*am+rn26*o2*bl10)/g7 bl12 = (rn34*h2o*am+rn30*h2o*bl7)/g12 bl8 = (rn35*am+rn12*bl7+rn33*bl12)*co2/g8 bl2 = rn4*o1/g2 bla2 = rn36*o3/g2 blb11 = rn28*an2*am bl11 = blb11*bl2/g11 bla11 = blb11*bla2/g11 bl9 = (rn13*bl7+rn14*bl8)*ano1/g9 blb3 = rn18*o3+rn27*o2*am blc3 = rn7*o3+rn15*o1*bl7+rn22*o3*bl8 bld3 = blb3*bl2+rn31*o3*bl12+rn37*o2*bl11 bl3 = (blc3+bld3)/g3 bla3 = (blb3*bla2+rn37*o2*bla11)/g3 bl4 = (rn8*co2*bl3+rn16*o1*bl8)/g4 bla4 = rn8*co2*bla3/g4 blb5 = rn23*ano2+rn9*ano1*bl4+rn10*ano1*bl3 blc5 = rn21*ano1*bl9+rn29*ano2*bl2 bl5 = (blb5+blc5)/g5 bla5 = (rn9*ano1*bla4+rn10*ano1*bla3+rn29*ano2*bla2)/g5 blb6 = hno3*(rn41+rn40*bl2+rn39*bl4+rn38*bl5) bl6 = (blb6+rn32*ano2*bl4+rn11*o3*bl5)/g6 blc6 = hno3*(rn40*bla2+rn39*bla4+rn38*bla5) bla6 = (blc6+rn32*ano2*bla4+rn11*o3*bla5)/g6 coe1 = (rn1*o2+rn2*an2)*o2+rn5*o1*bla3 coe2 = rn6*o1*bla4+rn19*o3*bla2 coe3 = rn5*o1*bl3+rn6*o1*bl4+rn19*o3*bl2 al1 = (coe1+coe2)/(g1-coe3) al2 = bl2*al1+bla2 al3 = bl3*al1+bla3 al4 = bl4*al1+bla4 al5 = bl5*al1+bla5 al6 = bl6*al1+bla6 al7 = bl7*al1 al8 = bl8*al1 al9 = bl9*al1 al10 = bl10*al1 al11 = bl11*al1+bla11 al12 = bl12*al1 sumn1 = al1+al2+al3+al4+al5+al6 sumn2 = al7+al8+al9+al10+al11+al12 alam = sumn1+sumn2 aneg = elec*(1.+alam) test = sumt-aneg test1 = abs(test) aqdfy=0.01 if (h < 10.) aqdfy=0.05 test2=test1-aqdfy*aneg if (test2 > 0.) then ! (was goto 41 from arithmetic if) ab1=elec*sumt/(1.+alam) ab1=abs(ab1) elec=sqrt(ab1) if (iabc < 100) goto 100 endif ! if (kht==1.or.mod(kht,10)==0) ! | write(6,"(' kht iabc elec alam sumt', ! | ' aneg test1 test2 ab1')") ! write(6,"(2i5,7e12.4))") kht,iabc,elec,alam,sumt,aneg,test1, ! | test2,ab1 bo2 = al1*elec bo1 = al2*elec bo3 = al3*elec bco3 = al4*elec bno2 = al5*elec bno3 = al6*elec bo4 = al7*elec bco4 = al8*elec bo2no = al9*elec bo2n2 = al10*elec bo1n2 = al11*elec bo2h2o = al12*elec aneg1 = elec*alam ! ! ar(54) is denmod module data ! Variables assigned to ar() are either args or local to ioncomp ar(1)=h ! height ar(2)=elec ! electrons ar(3)=alam ! lambda ar(4)=sumt ! total +ions ar(5)=aneg1 ! total -ions ar(6)=ano ! no+ ar(7)=anoh1 ! no+(h2o) ar(8)=anoh2 ! no+(h2o)2 ar(9)=anoh3 ! no+(h2o)3 ar(10)=anon1 ! no+(n2) ar(11)=anohn ! no+(h2o)(n2) ar(12)=anoh2n ! no+(h2o)2(n2) ar(13)=anoc1 ! no+(co2) ar(14)=anohc ! no+(h2o)(co2) ar(15)=anoh2c ! no+(h2O)2(co2) ar(16)=ao2 ! o2+ ar(17)=ao4 ! o4+ ar(18)=ao5 ! o5+ ar(19)=ao2h ! o2+(h2o) ar(20)=ahhoh ! h+(h2o)(oh) ar(21)=w1 ! h+(h2o) ar(22)=w2 ! h+(h2o)2 ar(23)=w3 ! h+(h2o)3 ar(24)=w4 ! h+(h2o)4 ar(25)=w5 ! h+(h2o)5 ar(26)=w6 ! h+(h2o)6 ar(27)=w7 ! h+(h2o)7 ar(28)=w8 ! h+(h2o)8 ar(29)=bo1 ! o- ar(30)=bo2 ! o2- ar(31)=bo3 ! o3- ar(32)=bo4 ! o4- ar(33)=bno2 ! no2- ar(34)=bno3 ! no3- ar(35)=bco3 ! co3- ar(36)=bco4 ! co4- ar(37)=bo2no ! o2-(no) ar(38)=bo2n2 ! o2-(n2) ar(39)=bo1n2 ! o-(n2) ar(40)=bo2h2O ! o2-(h2o) ar(41)=te ! temp ar(42)=am ! m ar(43)=co2 ! co2 ar(44)=h2o ! h2o ar(45)=o3 ! o3 ar(46)=o1 ! o ar(47)=o2s ! o2(singlet) ar(48)=ano1 ! no ar(49)=ano2 ! no2 ar(50)=qno ! q(no+) ar(51)=qo2 ! q(o2+) ar(52)=o2 ! o2 ar(53)=hno3 ! hno3 ar(54)=o2col ! o2 column ! write(6,"('ar=',/,(6e12.4))") ar qit=qio2p+qinop rk5=4.4e-10 alf1=r40*o2*am/(alph3*elec+r40*o2*am+rk5*ano1) alf2=r41/(r48*o1+alph2*elec+r41*h2o) alf3=(r44+r45)*h2o/((r44+r45)*h2o+alph2*elec) phoxo2=2.*(0.904*qit+r48*o1*(0.904*qit*alf1/(alph2*elec+r41*h2o | +(1.-alf1)*r48*o1)))*alf1*alf2*alf3/qit xl1=(r9*an2+r11*co2+r13*h2o)*am+alph2*elec xl2=r12*am+r15*h2o+alph2*elec xl3=r10*am+r14*co2+r16*h2o+alph2*elec gam1=r9*an2*am p1=gam1*ar(7) p2=r11*ar(7)*co2*am+r14*ar(11)*co2 bet1=1.-(alph2*elec/(xl1-am*(r10*gam1/xl3+r12*(r11*co2*am+gam1* | r14*co2)/xl3/xl2))) xl4=(r17*an2+r19*co2+r21*h2o)*am+alph2*elec xl5=r23*h2o+r20*am+alph2*elec xl6=r22*co2+r24*h2o+r18*am+alph2*elec gam2=r17*an2*am p3=gam2*ar(8) p4=r22*ar(12)*co2+r19*ar(8)*co2*am bet2=1.-(alph2*elec/(xl4-am*(r18*gam2/xl6+r20*(r19*co2*am+gam2* | r22*co2)/xl6/xl5))) bet3=r25/(r25*h2o+alph2*elec) phoxno=2.*bet1*bet2*bet3*(.096+.904*r49*ano1*(1.+alph1*r48*o1/ | (alph2*elec+r41*h2o+(1.-alph1)*r48*o1))/(alph2* | elec+r40*o2*am+r49*ano1)) phoxt=phoxo2+phoxno ! ! Print ar: if (iswtc) then ! if (kht==1) then ! write(6,"(/,'ioncomp: iday=',i4,' ut=',f8.2)") ! | iday,ut ! write(6,"(' k h elec alam sumt', ! | ' aneg1')") ! endif ! write(6,"(i4,f8.2,4e12.4)") kht,(ar(i),i=1,5) endif end subroutine ioncomp !----------------------------------------------------------------------- subroutine quart(a,root,delta,g,h,rl,cdelta,rlam) implicit none ! ! Args: real,intent(in) :: a(5) real,intent(out) :: root,delta,g,h,rl,cdelta complex,intent(out) :: rlam ! ! Local: complex :: cp real,parameter :: e=1.e-300 real :: a0,a1,a2,a3,a4,ri,rj,rk,ch,cg,p,q,pq ! A0=A(5) A1=A(4)/4. A2=A(3)/6. A3=A(2)/4. A4=A(1) G=A0**2*A3-3.*A0*A1*A2+2.*A1**3 H=A0*A2-A1**2 RI=A0*A4-4.*A1*A3+3.*A2**2 RJ=A0*(A2*A4-A3**2)-A1*(A1*A4-A3*A2)+A2*(A1*A3-A2**2) RK=A0**2*RI-3.*H**2 RL=12.*H**2-A0**2*RI CH=-RI/12. CG=RJ/4. DELTA=RI**3-27.*RJ**2 CDELTA=CG**2+4.*CH**3 CP=CDELTA+E CP=(.5*(CG+CP**.5)+E)**(1./3.) RLAM=-2.*REAL(CP) P=A0*RLAM+A1**2-A0*A2+E P=SQRT(P) Q=(2.*RLAM+A2)**2-A0*A4+E Q=SQRT(Q) PQ=2.*A1*RLAM+A1*A2-A0*A3+E P=SIGN(P,Q*PQ) ROOT=P-A1 ROOT=(ROOT+SQRT(ROOT**2-A0*(A2+2.*RLAM-Q)))/A0 end subroutine quart !----------------------------------------------------------------------- end module denmodule