c subroutine napart(tn,xo2,xh2o,xco2,xh,xh2,xo1,xo3,xne, + xo2p,xop,xn2,xnat,ncol,iden,zsb,zst,npart,names,fout, + iday,slt,glat,glon,spval,idebug) c c Partition Na, returning columns of up to 12 components plus airglow c (max npart=13) in fout(ncol,npart). The desired components are given c by calling routine in char*8 names(npart), with the following valid c names (case sensitive). (Components are returned in fout in same order c as that given in names). c c 'NaS ','NaO ','NaO3 ','NaO2 ','NaOH ', c 'NaCO3 ','NaHCO3 ','NaS+ ','NaN2+ ','NaCO2+ ', c 'NaH2O+ ','NaO+ ','NaEMIS ' c c ixnas ,ixnao ,ixnao3 ,ixnao2 ,ixnaoh ,ixnaco3,ixnahco3,ixnasp, c ixnan2p,ixnaco2p,ixnah2op,ixnaop ,ixnaemis c c These have the following dependencies: c TN, O2, H2O, CO2, H, H2, O1, O3, NE, O2+ [NO+], O+, N2 NAT c (since NO+ is not available on histories, use NO+ = (Ne - O+ - O2+)) c If iden=0, input species are MMR, otherwise they are in CM3 c (with exception of xne, xo2p, and xop, which are always cm3) c c NAT (xnat(ncol)) is total sodium, a prognostic from the model. c parameter(mxpart=13) c c Input: c real + tn(ncol) ,xo2(ncol) ,xh2o(ncol) ,xco2(ncol) ,xh(ncol), + xh2(ncol) ,xo1(ncol) ,xo3(ncol) ,xne(ncol) ,xo2p(ncol), + xop(ncol) ,xn2(ncol) ,xnat(ncol) c c Output: c real fout(ncol,npart),xxnat(ncol) real xnas(ncol) ,xnao(ncol) ,xnao3(ncol) ,xnao2(ncol), + xnaoh(ncol) ,xnaco3(ncol) ,xnahco3(ncol),xnasp(ncol), + xnan2p(ncol) ,xnaco2p(ncol),xnah2op(ncol),xnaop(ncol), + xnaemis(ncol) ! ! Externals: real,external :: getsza c c Work: c real s1(ncol),s2(ncol),s3(ncol),s4(ncol),s5(ncol),s6(ncol), + s7(ncol),s8(ncol),s9(ncol) real r1(ncol),r2(ncol),r3(ncol),r4(ncol),r5(ncol),r6(ncol), + r7(ncol),r8(ncol),r9(ncol) real an(ncol),bn(ncol),cn(ncol),dn(ncol),en(ncol),fn(ncol), + gn(ncol),hn(ncol),nn(ncol),aa(ncol),bb(ncol),cc(ncol), + dd(ncol),rt(ncol),denom(ncol),exps(ncol),bbdd(ncol) real rkn1(ncol), rkn2(ncol), rkn3a(ncol), rkn3b(ncol), + rkn4(ncol), rkn5(ncol), rkn6(ncol), rkn7(ncol), + rkn8(ncol), rkn9(ncol), rkn11(ncol), rkn13(ncol), + rkn15(ncol), rkn16(ncol), rkn17(ncol), rkn18(ncol), + rkn19(ncol), rkn20(ncol), rkn21(ncol), rkn22(ncol), + rkn24a(ncol), rkn24b(ncol), rkn25(ncol), rkn27a(ncol), + rkn27b(ncol), rkn28(ncol), rkn29(ncol), rkn30(ncol) real xjn10, xjn12, xjn14, xjn23, xjn26 real rmass(3) character*8 names(npart),na_names(mxpart) character*3 xunits data na_names / + 'NaS ','NaO ','NaO3 ','NaO2 ','NaOH ', + 'NaCO3 ','NaHCO3 ','NaS+ ','NaN2+ ','NaCO2+ ', + 'NaH2O+ ','NaO+ ','NaEMIS '/ data rmass/32.,16.,28./, p0/5.e-4/, boltz/1.38e-16/ data rna,rmo3,rmh2o,rmh,rmco2,rmo1,rmh2 + /23.,48., 18., 1., 44., 16., 2./ c xunits = "mmr" if (iden.gt.0) xunits = "cm3" c if (idebug.gt.0) then write(6,"(/72('-'))") write(6,"(/'napart: ncol=',i2,' npart=',i2,' iday=',i3, + ' slt=',f4.1,' glat=',f7.2,' glon=',f8.2, + ' napart inputs follow:')") ncol,npart,iday,slt,glat, + glon write(6,"('tn=',/(6e12.4))") tn write(6,"('xo2 (',a,') =',/(6e12.4))") xunits,xo2 write(6,"('xh2o (',a,') =',/(6e12.4))") xunits,xh2o write(6,"('xco2 (',a,') =',/(6e12.4))") xunits,xco2 write(6,"('xh (',a,') =',/(6e12.4))") xunits,xh write(6,"('xh2 (',a,') =',/(6e12.4))") xunits,xh2 write(6,"('xo1 (',a,') =',/(6e12.4))") xunits,xo1 write(6,"('xo3 (',a,') =',/(6e12.4))") xunits,xo3 write(6,"('xne (cm3) =',/(6e12.4))") xne write(6,"('xo2p (cm3) =',/(6e12.4))") xo2p write(6,"('xop (cm3) =',/(6e12.4))") xop write(6,"('xn2 (',a,') =',/(6e12.4))") xunits,xn2 write(6,"('xnat (',a,') =',/(6e12.4))") xunits,xnat endif if (npart.gt.mxpart) then write(6,"(/'>>> WARNING napart: bad npart=',i3, + ' (must be <= ',i2,') -- returning'/)") + npart,mxpart return endif if (ncol.le.0) then write(6,"(/'>>> WARNING napart: bad ncol=',i3, + ' -- fout not defined')") return endif c c Insure positive o3 and h2o: c do k=1,ncol nbad_xo3 = 0 if (xo3(k).lt.0.) then xo3(k) = 1.e-20 nbad_xo3 = nbad_xo3+1 endif if (idebug.gt.0.and.nbad_xo3.gt.0) then write(6,"('NOTE NAPART: nbad_xo3 (number of xo3 < 0)=',i2)") + nbad_xo3 endif c nbad_xh2o = 0 if (xh2o(k).lt.0.) then xh2o(k) = 1.e-20 nbad_xh2o = nbad_xh2o+1 endif if (idebug.gt.0.and.nbad_xh2o.gt.0) then write(6,"('NOTE NAPART: nbad_xh2o (number of xh2o < 0)=',i2)") + nbad_xh2o endif enddo c c Use solar zenith angle za to determine visibility of sun for c photochem reactions: c za = getsza(iday,slt,glat,glon) day = 0. if (za.le.90.) day = 1. xjn10 = 4.e-3*day xjn12 = 1.e-3*day xjn14 = 2.e-5*day xjn23 = 0. xjn26 = 1.e-4*day if (idebug.gt.0) then write(6,"(/'napart: za=',f7.2,' xjn10=',e12.4,' xjn12=', + e12.4,' xjn14=',e12.4,/' xjn23=',e12.4,' xjn26=',e12.4)") + za,xjn10,xjn12,xjn14,xjn23,xjn26 endif c c Sodium reactions: c do k=1,ncol rkn1(k) = 1.1e-9*exp(-116./tn(k)) rkn2(k) = 2.2e-10*sqrt(tn(k)/200.) rkn3a(k) =1.1e-9*exp(-568./tn(k)) rkn3b(k)= 3.2e-10*exp(-550./tn(k)) rkn4(k)= 4.4e-10*exp(-507./tn(k)) rkn5(k)= 5.3e-30*(200./tn(k)) rkn6(k)= 1.3e-27*(200./tn(k)) rkn7(k)= 5.0e-30*(200./tn(k))**1.22 rkn8(k)= 5.0e-10*exp(-940./tn(k)) rkn9(k)= 1.e-9*exp(-1000./tn(k)) rkn11(k)= 4.e-11*exp(-550./tn(k)) rkn13(k)= 1.9e-28*(200./tn(k)) rkn15(k)= 1.4e-9 rkn16(k)= 1.0e-9 rkn17(k)= 2.6e-30*(200./tn(k))**2.19 rkn18(k)= 1.e-6*(200./tn(k))**0.5 rkn19(k)= 2.5e-10*sqrt(tn(k)/200.) rkn20(k)= 5.0e-10*exp(-1200./tn(k)) rkn21(k)= 1.0e-9*exp(-1400./tn(k)) rkn22(k)= 1.1e-11*exp(-910./tn(k)) rkn24a(k)=1.1e-9*exp(-1100./tn(k)) rkn24b(k)=1.1e-9*exp(-1400./tn(k)) rkn25(k)=3.e-10*exp(-668./tn(k)) rkn27a(k) = 8.e-10 rkn27b(k) = 8.e-10 rkn28(k) = 6.e-10 rkn29(k) = 8.e-10 rkn30(k)= 1.e-6*(200./tn(k))**0.5 enddo if (idebug.gt.0) then write(6,"(/'napart reactions:')") write(6,"('rkn1=',/(6e12.4))") rkn1 write(6,"('rkn2=',/(6e12.4))") rkn2 write(6,"('rkn3a=',/(6e12.4))") rkn3a write(6,"('rkn3b=',/(6e12.4))") rkn3b write(6,"('rkn4=',/(6e12.4))") rkn4 write(6,"('rkn5=',/(6e12.4))") rkn5 write(6,"('rkn6=',/(6e12.4))") rkn6 write(6,"('rkn7=',/(6e12.4))") rkn7 write(6,"('rkn8=',/(6e12.4))") rkn8 write(6,"('rkn9=',/(6e12.4))") rkn9 write(6,"('rkn11=',/(6e12.4))") rkn11 write(6,"('rkn13=',/(6e12.4))") rkn13 write(6,"('rkn15=',/(6e12.4))") rkn15 write(6,"('rkn16=',/(6e12.4))") rkn16 write(6,"('rkn17=',/(6e12.4))") rkn17 write(6,"('rkn18=',/(6e12.4))") rkn18 write(6,"('rkn19=',/(6e12.4))") rkn19 write(6,"('rkn20=',/(6e12.4))") rkn20 write(6,"('rkn21=',/(6e12.4))") rkn21 write(6,"('rkn22=',/(6e12.4))") rkn22 write(6,"('rkn24a=',/(6e12.4))") rkn24a write(6,"('rkn24b=',/(6e12.4))") rkn24b write(6,"('rkn25=',/(6e12.4))") rkn25 write(6,"('rkn27a=',/(6e12.4))") rkn27a write(6,"('rkn27b=',/(6e12.4))") rkn27b write(6,"('rkn28=',/(6e12.4))") rkn28 write(6,"('rkn29=',/(6e12.4))") rkn29 write(6,"('rkn30=',/(6e12.4))") rkn30 endif c c Convert input species from mmr to cm3 if necessary: c s1 = p0*exp(-zp) / barm*kT c if (iden.le.0) then ! input species in mmr -- convert to cm3 ds = (zst-zsb)/(ncol-1) exps(1) = exp(-zsb) expds = exp(-ds) do k=2,ncol exps(k) = exps(k-1)*expds enddo if (idebug.gt.0) write(6,"('zsb=',f6.2,' zst=',f6.2,' ds=',f6.2, + ' expds=',f6.2,' exps=',/(6e12.4))") zsb,zst,ds,expds,exps do k=1,ncol s1(k) = p0*exps(k)/(boltz*tn(k)* + (xo2(k)/rmass(1) + xo1(k)/rmass(2) + xn2(k)/rmass(3))) s2(k) = xn2(k)/rmass(3)*s1(k) s3(k) = xo2(k)/rmass(1)*s1(k) s4(k) = xo1(k)/rmo1*s1(k) s5(k) = xco2(k)/rmco2*s1(k) s6(k) = xh(k)/rmh*s1(k) s7(k) = xh2o(k)/rmh2o*s1(k) s8(k) = xo3(k)/rmo3*s1(k) s9(k) = xh2(k)/rmh2*s1(k) xxnat(k)=xnat(k)*s1(k)/rna enddo else ! input species already in cm3 s2(:) = xn2(:) s3(:) = xo2(:) s4(:) = xo1(:) s5(:) = xco2(:) s6(:) = xh(:) s7(:) = xh2o(:) s8(:) = xo3(:) s9(:) = xh2(:) xxnat(:) = xnat(:) endif if (idebug.gt.0) then write(6,"(/'napart: input number densities:')") if (iden.eq.0) + write(6,"('s1 (p0*exp(-z)/(barm*kt))=',/(6e12.4))") s1 write(6,"('s2 (n2)=',/(6e12.4))") s2 write(6,"('s3 (o2)=',/(6e12.4))") s3 write(6,"('s4 (o1)=',/(6e12.4))") s4 write(6,"('s5 (co2)=',/(6e12.4))") s5 write(6,"('s6 (h)=',/(6e12.4))") s6 write(6,"('s7 (h2o)=',/(6e12.4))") s7 write(6,"('s8 (o3)=',/(6e12.4))") s8 write(6,"('s9 (h2)=',/(6e12.4))") s9 write(6,"('xxnat (nat)=',/(6e12.4))") xxnat endif do k=1,ncol an(k) = rkn5(k)*s3(k)*s2(k)/(rkn19(k)*s4(k)+xjn26) bn(k) = rkn6(k)*s5(k)*s2(k)/(rkn20(k)*s4(k)+ + rkn21(k)*s6(k)) cn(k) = (rkn4(k)*s7(k)+rkn24a(k)*s9(k)+rkn21(k) + *s6(k)*bn(k))/(rkn11(k)*s6(k)+xjn12 + +rkn13(k)*s2(k)*s5(k)) denom(k) = rkn8(k)*s4(k)+rkn9(k)*s6(k)+xjn10 dn(k) = (rkn3a(k)*s8(k)+rkn20(k)*s4(k)*bn(k))/denom(k) en(k) = rkn7(k)*s3(k)*s2(k)/denom(k) fn(k) = rkn13(k)*s5(k)*s2(k)*cn(k)/(rkn22(k)*s6(k)) denom(k) = rkn2(k)*s4(k)+(rkn3a(k)+rkn3b(k))*s8(k)+rkn4(k)* + s7(k)+rkn5(k)*s3(k)*s2(k)+rkn6(k)*s5(k)*s2(k)+ + (rkn24a(k)+rkn24b(k))*s9(k)+rkn25(k)*s6(k) gn(k) = (rkn1(k)*s8(k)+rkn8(k)*s4(k)*en(k))/denom(k) hn(k) = (xjn26*an(k)+rkn8(k)*s4(k)*dn(k))/denom(k) nn(k) = gn(k)/(1.-hn(k)) r1(k) = nn(k) r2(k) = an(k)*nn(k) r3(k) = bn(k)*nn(k) r4(k) = cn(k)*nn(k) r5(k) = dn(k)*nn(k)+en(k) r6(k) = fn(k)*nn(k) aa(k)= xjn14+rkn15(k)*xo2p(k)+rkn16(k)* + (xne(k)-xop(k)-xo2p(k)) ! f(i,nnopk) bb(k) = rkn29(k)*s4(k)*rkn28(k)*s4(k)/(rkn29(k) + *s4(k)+rkn30(k)*xne(k)) cc(k) = rkn17(k)*s2(k)*s2(k) dd(k) = rkn27a(k)*s5(k) + rkn27b(k)*s7(k) + rkn28(k)*s4(k) + + rkn30(k)*xne(k) r7(k) = aa(k)/(cc(k)*(1.-bb(k)/dd(k))) r9(k) = cc(k)/dd(k) r8(k) = (1.+rkn28(k)*s4(k)/(rkn29(k)*s4(k)+ + rkn30(k)*xne(k))+(rkn27a(k)*s5(k)+ + rkn27b(k)*s7(k))/(rkn30(k)*xne(k))) + *r7(k)*r9(k) rt(k)=1./(1.+r1(k)+r2(k)+r3(k)+r4(k)+r5(k)+r6(k)+r7(k)+r8(k)) enddo if (idebug.gt.0) then write(6,"(/'napart an,bn,cn,dn,en,fn,gn,hn,nn=')") write(6,"('an=',/(6e12.4))") an write(6,"('bn=',/(6e12.4))") bn write(6,"('cn=',/(6e12.4))") cn write(6,"('dn=',/(6e12.4))") dn write(6,"('en=',/(6e12.4))") en write(6,"('fn=',/(6e12.4))") fn write(6,"('gn=',/(6e12.4))") gn write(6,"('hn=',/(6e12.4))") hn write(6,"('nn=',/(6e12.4))") nn endif if (idebug.gt.0) then write(6,"(/'napart aa,bb,cc,dd:')") write(6,"('aa=',/(6e12.4))") aa write(6,"('bb=',/(6e12.4))") bb write(6,"('cc=',/(6e12.4))") cc write(6,"('dd=',/(6e12.4))") dd do k=1,ncol bbdd(k) = bb(k)/dd(k) enddo write(6,"('bb/dd=',/(6e12.4))") bbdd endif if (idebug.gt.0) then write(6,"(/'napart r1-9,t:')") write(6,"('r1=',/(6e12.4))") r1 write(6,"('r2=',/(6e12.4))") r2 write(6,"('r3=',/(6e12.4))") r3 write(6,"('r4=',/(6e12.4))") r4 write(6,"('r5=',/(6e12.4))") r5 write(6,"('r6=',/(6e12.4))") r6 write(6,"('r7=',/(6e12.4))") r7 write(6,"('r8=',/(6e12.4))") r8 write(6,"('r9=',/(6e12.4))") r9 write(6,"('rt=',/(6e12.4))") rt write(6,"('xxnat=',/(6e12.4))") xxnat endif c c Neutrals output: c do k=1,ncol xnas(k)=rt(k)*xxnat(k) xnao(k)=r1(k)*xnas(k) xnao3(k)=r2(k)*xnas(k) xnaco3(k)=r3(k)*xnas(k) xnaoh(k)=r4(k)*xnas(k) xnao2(k)=r5(k)*xnas(k) xnahco3(k)=r6(k)*xnas(k) c c Ions output: c xnan2p(k)=r9(k)*r7(k)*xnas(k) xnaop(k) =rkn28(k)*s4(k)/(rkn29(k)*s4(k)+rkn30(k) + *xne(k))*xnan2p(k) xnasp(k) = (aa(k)*xnas(k)+rkn29(k)*s4(k)*xnaop(k))/cc(k) xnaco2p(k)=rkn27a(k)*s5(k)/(rkn30(k)*xne(k))*xnan2p(k) xnah2op(k)=rkn27b(k)*s7(k)/(rkn30(k)*xne(k))*xnan2p(k) c c Sodium airglow (height integrate and divide by 1.e6 to get Rayleighs) c (better to multiply by 1.e-6) c xnaemis(k) = rkn2(k)*xnao(k)*s4(k)*0.2 enddo if (idebug.gt.0) then write(6,"(/'napart neutral output (7 species):')") write(6,"('xnas=',/(6e12.4))") xnas write(6,"('xnao=',/(6e12.4))") xnao write(6,"('xnao3=',/(6e12.4))") xnao3 write(6,"('xnaco3=',/(6e12.4))") xnaco3 write(6,"('xnaoh=',/(6e12.4))") xnaoh write(6,"('xnao2=',/(6e12.4))") xnao2 write(6,"('xnahco3=',/(6e12.4))") xnahco3 c write(6,"(/'napart ion output (5 species):')") write(6,"('xnan2p=',/(6e12.4))") xnan2p write(6,"('xnaop=',/(6e12.4))") xnaop write(6,"('xnasp=',/(6e12.4))") xnasp write(6,"('xnaco2p=',/(6e12.4))") xnaco2p write(6,"('xnah2op=',/(6e12.4))") xnah2op endif c c Define partition(s) in fout: c do ip=1,npart ixp = ixfindc(na_names,mxpart,names(ip)) if (ixp.le.0) then write(6,"('>>> WARNING napart: unknown sodium name=',a,/ + ' Returning spval=',e9.2,/ + ' Valid Na partition names = ',/(8a8))") + names(ip),spval,na_names fout(:,ip) = spval goto 100 endif c c + 'NaS ','NaO ','NaO3 ','NaO2 ','NaOH ', c + 'NaCO3 ','NaHCO3 ','NaS+ ','NaN2+ ','NaCO2+ ', c + 'NaH2O+ ','NaO+ ','NaEMIS '/ c if (names(ip).eq."NaS ") fout(:,ip) = xnas(:) if (names(ip).eq."NaO ") fout(:,ip) = xnao(:) if (names(ip).eq."NaO3 ") fout(:,ip) = xnao3(:) if (names(ip).eq."NaO2 ") fout(:,ip) = xnao2(:) if (names(ip).eq."NaOH ") fout(:,ip) = xnaoh(:) if (names(ip).eq."NaCO3 ") fout(:,ip) = xnaco3(:) if (names(ip).eq."NaHCO3 ") fout(:,ip) = xnahco3(:) if (names(ip).eq."NaS+ ") fout(:,ip) = xnasp(:) if (names(ip).eq."NaN2+ ") fout(:,ip) = xnan2p(:) if (names(ip).eq."NaCO2+ ") fout(:,ip) = xnaco2p(:) if (names(ip).eq."NaH2O+ ") fout(:,ip) = xnah2op(:) if (names(ip).eq."NaO+ ") fout(:,ip) = xnaop(:) if (names(ip).eq."NaEMIS ") fout(:,ip) = xnaemis(:) 100 continue enddo if (idebug.gt.0) write(6,"(72('-')/)") return end