! module composition ! ! Chemistry and composition (formerly in compo.F) ! use params,only: nlev,zpb,dzp use flds_fields,only: edy ! (was also difk in old code) use cons,only: rmo2,rmo,rmn2 use input,only: istep,nstep implicit none real :: xno_lb, xno3_lb, ps1b,ps2b, expdz,expdzinv ! ! fs1, et.al. were formerly in /source/ (compo.F): ! Defined by sub chemistry below. real,dimension(nlev) :: exps,difk, | fs1,fs2,fs11,fs12,fs21,fs22,t0 real :: b(2,2),fb(2),phi(2,3) ! ! phno, et.al. were formerly in /photoo/ (compo.F): ! Defined by sub photo below. real :: phno,phch4,phh2,phco2,phco,phh2o,phhe,pharg,phn4s ! ! Selected species densities in mass-mixing ratios (psi's): ! (analogous to the f-array species of old code, e.g., f(i,nps)) real,dimension(nlev) :: | ps_o2 ,ps_ox ,ps_n2 ,ps_ar ,ps_he ,ps_noz ,ps_n4s, | ps_ch4,ps_h2 ,ps_co2 ,ps_co ,ps_hox ,ps_h2o ,ps_nat ! ! Selected species volume number densities (cm-3): ! (analogous to fno, fno2, etc in old code) real,dimension(nlev) :: | xn_o2 ,xn_o ,xn_n2 ,xn_ar ,xn_he ,xn_noz ,xn_n4s, | xn_ch4,xn_h2 ,xn_co2 ,xn_co ,xn_hox ,xn_h2o ,xn_nat ! real,parameter :: | rmo2inv=1./rmo2, | rmoinv =1./rmo, | rmn2inv=1./rmn2, | t00=273., | tau=1.86e+3, | dfact = 1. real :: stepinv real :: phij,phip,phie,fluxhox,sumnap real,dimension(nlev) :: rk12_comp, rk13_comp, rk14_comp, rk15_comp ! ! Function expo checks exp() arg and flushes the result to 0 if arg ! is too small, or to largest value allowed if arg is too big, avoiding ! an fpe underflow or overflow. This function is used only if EXPO ! is set (i.e., when fpe traps are set, see Makefile). Use of expo() ! does slow the code down considerably. ! real,external :: expo ! util.F contains !----------------------------------------------------------------------- subroutine chemistry(iday,ut,dt) use flds_atmos,only: xno,xno2,xnn2,xnhox,xnh,xnoh,xnho2,xnnoz, | xnno,xnno2,xnox,xno3,xn2d,xno21d,xnco2,xno21s,xno1d,xnh2,xnh2o, | xnch4,ch3,ch3o2,ch3ooh,ch3o,ch2o,hco,xnh2o2,xncl,xnco,xn4s, | xnclo,tn,xnhe,xnarg,rho,xnhtot use flds_ionzrt,only: qtin,qnspe,qnp,xjch3oo,xjch2oa,xjch2ob, | dinolya use flds_dissoc,only: pnoeuv,do3har,do2lya,do2src,do2euv,dch4t, | dh2o2t,dch4t,dco2t,dh2oeuv,dh2olya,dh2ot,dh2osrb,dh2osrc, | dno2t,dnoeuv,dnosrb,do2t,do3t,diop,do2srb,dio2p,do2hz use flds_ratecoef,only: rk1,rk2,rk3,rk4,rk5,rk6,rk7,rk8,rk9,rk10, | rk11,rk12,rk13,rk14,rk15,rk16,rk23,rk24,rk27,alp1,alp2,alp3, | rb1,rb2,rb3,rb4,rb5,rb6,rb7,rb8,rb9,rb10,rb11,rb12,rb13,rkm1, | rkm2a,rkm2b,rkm3,rkm4,rkm5a,rkm5b,rkm6,rkm7a,rkm7b,rkm8a,rkm9, | rkm10,rkm11,rkm13,rkm14,rkm15,rkm16,rkm17,rkm20,rkm21,rkm22, | rkm23,rkm24,rkm25,rkm26,rkm27,rkm28,rkm29,rkm30,rkm31,rkm32, | rkm33,rkm34,rkm35,rkm36,rkm37,rkm38,rkm39,rkm40,rkm41,rkm42, | adl,a1d,asg,gam1,gam2a,gam2b,gam3,gam4,gam5,gam6,gam7,gam8, | gam9,gam10,gam11,gam12,gam13,gam14,gam15,del1,del2,del3 use flds_ionatm,only: xne,xinop,xiop,xin2p,xio2p,xihp,xinp, | xiop2d,xiop2p,xiop4s,te use flds_prodloss,only: xln2d,pn2d,po21dl,xlo21dl,po1d,xlo1d, | pch3,xlch3,po21sg,xlo21sg,pch3o2,xlch3o2,pch3oo,xlch3oo, | pch3o,xlch3o,pch2o,xlch2o,phco,xlhco,ph2o2,xlh2o2,pch4, | xlch4,pco2,xlco2,pco,xlco,ph2,xlh2,ph2o,xlh2o,xloh,ph,xlh, | xlho2,pho2,poh,xlhoxt,phox,xlhox,xlno,pno,xlno2,pno2,pn4s, | xln4s,xlnoxt,xlnox,pnox,po,xlo,po3,xlo3,pna use flds_fields,only: o3or,hhoxr,ooxr,rnonoz,rno2no,ohhr,ho2hr use flds_sodium,only: xnat,pnas use flds_heat,only: xno21dc,xno21sc,xmetah use flds_modelz,only: zpht,gz use denmodule,only: qhoxo2,qhoxno use cons,only: erg,brn2d,brn4s,xo3lb,xn4slb use solarheat,only: xj762 use flds_hist,only: f_hist,nf_hist,addfhist ! use input,only: istep,nstep implicit none ! ! Args: integer,intent(in) :: iday real,intent(in) :: ut,dt ! ! Local: integer :: k real,parameter :: eps = 1. real,dimension(nlev) :: phoxic,ratio1,ratio2,ratio3,ratio4,ratio5, | ratio6,ratio7 real,dimension(nlev) :: | qn2d,qn4s,dn2, ! formerly in /chemtr/ | a4,a5,b1,b2,b3,c1,c2,c3,d1 ! formerly prodloss.h real :: xm,dum,wl1s,wl1d,ene1s,ene1d,as,zps,h1s, aa,bb,cc,bb2 character(len=80) :: char80 ! ! falpxxx were formerly in /branch/, but were used only in sub chemstry. ! Here they are local to sub chemistry: real,parameter :: | falp14s=0.2, falp12d=0.8, falp23p=0.15, | falp21d=0.85, falp34s=0.1, falp32d=0.9 ! ! Main levels loop: do k=1,nlev ! was do 1 xm=xno(k)+xno2(k)+xnn2(k) xnhox(k)=xnh(k)+xnoh(k)+xnho2(k) xnox(k)=xno(k)+xno3(k) xnnoz(k)=xnno(k)+xnno2(k) phoxic(k)=(qhoxo2(k)+qhoxno(k))*qtin(k) ! ! Metastable species n(2d): qn2d(k)=qnspe(k)*brn2d+pnoeuv(k)+qnp(k)*brn2d qn4s(k)=qnspe(k)*brn4s+pnoeuv(k)+qnp(k)*brn4s dn2(k)=qn2d(k)+qn4s(k) pn2d(k)=qn2d(k)+rk3(k)*xin2p(k)*xno(k)+alp1(k)*xinop(k)*xne(k)* | falp12d+alp3(k)*xin2p(k)*xne(k)*falp32d xln2d(k)=rb2(k)*xno2(k)+rb4(k)*xno(k)+rb5(k)* | xne(k)+rb6(k)*xnno(k)+rb7(k)+rk10(k)*xiop(k) xn2d(k)=pn2d(k)/xln2d(k) ! ! Metastable species o2(1sg): po21sg(k)=xj762*xno2(k)+rkm2a(k)*xno1d(k)*xno2(k) xlo21sg(k)=rkm13(k)*xnn2(k)+rkm14(k)*xnco2(k)+rkm15(k)*xno3(k)+ | rkm16(k)*xno(k)+rkm17(k)*xno2(k)+asg xno21s(k) = po21sg(k)/xlo21sg(k) ! ! Metastable species o2(1dl): po21dl(k)=do3har(k)*xno3(k)*0.9+(rkm13(k)*xnn2(k)+rkm14(k)* | xnco2(k)+rkm15(k)*xno3(k)+rkm16(k)*xno(k)+ | rkm17(k)*xno2(k))*xno21s(k) xlo21dl(k)=rkm9(k)*xno2(k)+rkm10(k)*xnn2(k)+rkm11(k)*xno(k)+adl xno21d(k)=po21dl(k)/xlo21dl(k) ! ! Metastable species o(1d): po1d(k)=do2src(k)*xno2(k)+do2lya(k)*xno2(k)+do3har(k)*xno3(k)+ | do2euv(k)*xno2(k)+ | rb2(k)*xn2d(k)*xno2(k)*0.1+alp2(k)*xio2p(k)* | xne(k)*0.85 xlo1d(k)= rkm1(k)*xnn2(k)+(rkm2a(k)+rkm2b(k))*xno2(k)+ | rkm3(k)*xnh2o(k)+rkm8a(k)*xno(k)+ | rkm4(k)*xnh2(k)+(rkm5a(k)+rkm5b(k))*xnch4(k)+rkm6(k)* | xnco2(k)+(rkm7a(k)+rkm7b(k))*xno3(k)+a1d xno1d(k)=po1d(k)/xlo1d(k) ! ! Methane oxidation: ch3,ch3o2,ch3ooh,ch3o,ch2o,cho ! ch3: pch3(k)=(gam1(k)*xnoh(k)+gam2a(k)*xno1d(k)+gam13(k)*xno(k)+ | dch4t(k))*xnch4(k) xlch3(k) = gam3(k)*xno2(k)*xm+gam4(k)*xno(k) ch3(k)=pch3(k)/xlch3(k) ! ! ch3o2: pch3o2(k)=gam3(k)*xno2(k)*xm*ch3(k) xlch3o2(k)=gam5(k)*xnno(k)+gam6(k)*xnho2(k)-(gam8(k)*xnoh(k) | *gam6(k)*xnho2(k))/(xjch3oo(k)+gam8(k)*xnoh(k)) ch3o2(k)=pch3o2(k)/xlch3o2(k) ! ! ch3ooh: pch3oo(k)=gam6(k)*xnho2(k)*ch3o2(k) xlch3oo(k)=xjch3oo(k)+gam8(k)*xnoh(k) ch3ooh(k)=pch3oo(k)/xlch3oo(k) ! ! ch3o: pch3o(k)=gam5(k)*xnno(k)*ch3o2(k)+xjch3oo(k)*ch3ooh(k) xlch3o(k) = gam9(k)*xno2(k) ch3o(k)=pch3o(k)/xlch3o(k) ! ! ch2o: pch2o(k)=gam2b(k)*xno1d(k)*xnch4(k)+gam4(k)*xno(k)*ch3(k)+ | gam9(k)*xno2(k)*ch3o(k) xlch2o(k)=xjch2oa(k)+xjch2ob(k)+gam10(k)*xnoh(k)+gam11(k)*xno(k) ch2o(k)=pch2o(k)/xlch2o(k) ! ! hco (cho?): phco(k)=(xjch2oa(k)+gam10(k)*xnoh(k)+gam11(k)*xno(k))*ch2o(k) xlhco(k)=gam12(k)*xno2(k) hco(k)=phco(k)/xlhco(k) ! ! h2o2: ph2o2(k)=rkm35(k)*xnho2(k)*xnho2(k) xlh2o2(k)=dh2o2t(k)+rkm27(k)*xno(k)+rkm32(k)*xnoh(k) xnh2o2(k)=ph2o2(k)/xlh2o2(k) ! ! ch4 production and loss: pch4(k)=0. xlch4(k)=dch4t(k)+gam1(k)*xnoh(k)+gam13(k)*xno(k)+ | (rkm5a(k)+rkm5b(k))*xno1d(k)+del3(k)*xncl(k) ! ! co2 production and loss: pco2(k)=rkm42(k)*xnco(k)*xnoh(k)+gam14(k)*xno(k)*xm*xnco(k) xlco2(k)=dco2t(k)+rk13(k)*xiop(k)+gam15(k)*xno1d(k) ! ! co production and loss: pco(k)=dco2t(k)*xnco2(k)+rk13(k)*xiop(k)*xnco2(k)+ | (gam1(k)*xnoh(k)+(rkm5a(k)+rkm5b(k))*xno1d(k)+ | gam13(k)*xno(k)+del3(k)*xncl(k))*xnch4(k)+ | gam15(k)*xno1d(k)*xnco2(k) xlco(k)=gam14(k)*xno(k)*xm+rkm42(k)*xnoh(k) ! ! h2o production and loss: ph2o(k)=rkm30(k)*xnoh(k)*xnoh(k)+rkm31(k)*xnoh(k)*xnho2(k) | +rkm32(k)*xnoh(k)*xnh2o2(k)+rkm33(k)*xnoh(k)*xnh2(k) | +rkm40(k)*xnh(k)*xnho2(k)+gam1(k)*xnch4(k)*xnoh(k) | +gam8(k)*ch3ooh(k)*xnoh(k)+gam10(k)*ch2o(k) | *xnoh(k)+del3(k)*xncl(k)*xnch4(k) xlh2o(k)=dh2ot(k)+rkm3(k)*xno1d(k)+rk15(k)*xiop(k) | +phoxic(k)*0.5 ! ! h2 production and loss: ph2(k)=(dh2olya(k)+dh2oeuv(k))*xnh2o(k) | +rkm38(k)*xnh(k)*xnho2(k) | +rkm41(k)*xnh(k)*xnh(k)*xm+gam2b(k)*xnch4(k)*xno1d(k) | +xjch2ob(k)*ch2o(k) xlh2(k)=rkm4(k)*xno1d(k)+rkm28(k)*xno(k)+rkm33(k)*xnoh(k) | +rk14(k)*xiop(k) ! ! ho2 production and loss: pho2(k)=rkm27(k)*xno(k)*xnh2o2(k)+rkm29(k)*xno3(k)*xnoh(k) | +rkm32(k)*xnoh(k)*xnh2o2(k)+rkm36(k)*xnh(k)*xno2(k)* | xm+gam9(k)*ch3o(k)*xno2(k)+gam12(k)*hco(k)* xno2(k) xlho2(k)=rkm26(k)*xno(k)+rkm31(k)*xnoh(k)+rkm34(k)*xno3(k) | +(rkm38(k)+rkm39(k)+rkm40(k))*xnh(k)+rb10(k)* | xnno(k)+gam6(k)*ch3o2(k)+2.*rkm35(k)*xnho2(k) ! ! oh production and loss ! 8/23/05 btf: change rkm33 to rkm3 in oh production (also in old code). poh(k)=2.*rkm3(k)*xnh2o(k)*xno1d(k)+rkm4(k)*xnh2(k)*xno1d(k) | +rkm26(k)*xno(k)*xnho2(k)+rkm27(k)*xno(k)*xnh2o2(k) | +rkm28(k)*xno(k)*xnh2(k)+rkm34(k)*xno3(k)*xnho2(k) | +rkm37(k)*xno3(k)*xnh(k)+2.*rkm39(k)*xnh(k)*xnho2(k) | +gam2a(k)*xnch4(k)*xno1d(k)+gam11(k)*ch2o(k)*xno(k) | +gam13(k)*xnch4(k)*xno(k)+xjch3oo(k)*ch3ooh(k) | +rb10(k)*xnno(k)*xnho2(k)+2.*dh2o2t(k)*xnh2o2(k) | +dh2ot(k)*xnh2o(k) xloh(k)=rkm25(k)*xno(k)+rkm29(k)*xno3(k)+rkm31(k)*xnho2(k) | +rkm32(k)*xnh2o2(k)+rkm33(k)*xnh2(k) | +rkm42(k)*xnco(k)+gam1(k)*xnch4(k)+gam8(k)* | ch3ooh(k)+gam10(k)*ch2o(k)+rb8(k)*xn4s(k) | +2.*rkm30(k)*xnoh(k) ! ! h production and loss ph(k)=rkm4(k)*xnh2(k)*xno1d(k)+rkm25(k)*xno(k)*xnoh(k) | +rkm28(k)*xno(k)*xnh2(k)+rkm33(k)*xnoh(k)*xnh2(k) | +rkm42(k)*xnco(k)*xnoh(k)+gam4(k)*ch3(k)*xno(k) | +rb8(k)*xn4s(k)*xnoh(k)+rk11(k)*xihp(k)*xno(k) | +rk14(k)*xiop(k)*xnh2(k)+dch4t(k)*xnch4(k) | +xjch2oa(k)*ch2o(k) xlh(k)=rkm36(k)*xno2(k)*xm+rkm37(k)*xno3(k)+(rkm38(k)+ | rkm39(k)+rkm40(k))*xnho2(k)+rk12(k)*xiop(k) | +2.*rkm41(k)*xm*xnh(k) ! ! hox production and loss phox(k)=2.*rkm27(k)*xno(k)*xnh2o2(k)+2.*rkm3(k)*xno1d(k)* | xnh2o(k)+2.*rkm4(k)*xno1d(k)*xnh2(k)+2.*rkm28(k)* | xno(k)*xnh2(k)+gam9(k)*ch3o(k)*xno2(k)+gam12(k)* | hco(k)*xno2(k)+gam2a(k)*xnch4(k)*xno1d(k)+ | gam11(k)*ch2o(k)*xno(k)+gam13(k)*xnch4(k)*xno(k) | +gam4(k)*ch3(k)*xno(k)+rk11(k)*xihp(k)*xno(k) | +rk14(k)*xiop(k)*xnh2(k)+xjch3oo(k)*ch3ooh(k) | +2.*dh2o2t(k)*xnh2o2(k)+(dh2osrc(k)+dh2osrb(k))* | xnh2o(k)+dch4t(k)*xnch4(k)+xjch2oa(k)*ch2o(k) | +phoxic(k)*xnh2o(k) xlhoxt(k)=2.*rkm31(k)*xnoh(k)*xnho2(k)+2.*rkm38(k)*xnh(k)* | xnho2(k)+2.*rkm40(k)*xnh(k)*xnho2(k)+2.*rkm35(k) | *xnho2(k)*xnho2(k)+2.*rkm30(k)*xnoh(k)*xnoh(k) | +2.*rkm41(k)*xm*xnh(k)*xnh(k)+gam6(k)*ch3o2(k)* | xnho2(k)+gam1(k)*xnch4(k)*xnoh(k)+gam8(k)* | ch3ooh(k)*xnoh(k)+gam10(k)*ch2o(k)*xnoh(k) | +rk12(k)*xiop(k)*xnh(k) ! ! hox partitioning: ratio1=oh/h, ratio2=ho2/h, ratio3=h/hox ! ratio1(k)=(rkm36(k)*xno2(k)*xm+rkm37(k)*xno3(k))/(rkm25(k)* | xno(k)+rkm33(k)*xnh2(k)+rkm42(k)*xnco(k)+rb8(k)* | xn4s(k)) ratio2(k)=(rkm36(k)*xno2(k)*xm+rkm29(k)*xno3(k)*ratio1(k))/ | (rkm26(k)*xno(k)+rkm34(k)*xno3(k)+rb10(k)*xnno(k)) ratio3(k)=1./(1.+ratio1(k)+ratio2(k)) xlhox(k)=(2.*rkm31(k)*ratio1(k)*ratio2(k)+2.*(rkm38(k)+ | rkm40(k))*ratio2(k)+2.*rkm35(k)*ratio2(k)*ratio2(k)+ | 2.*rkm30(k)*ratio1(k)*ratio1(k)+2.*rkm41(k)*xm)* | ratio3(k)*ratio3(k)*xnhox(k)+(gam6(k)*ch3o2(k)* | ratio2(k)+(gam1(k)*xnch4(k)+gam8(k)*ch3ooh(k)+ | gam10(k)*ch2o(k))*ratio1(k)+rk12(k)*xiop(k))*ratio3(k) ! ! n(4s) production and loss: pn4s(k)=qn4s(k)+rk2(k)*xiop(k)*xnn2(k)+rk6(k)*xinp(k) | *xno2(k)+rk8(k)*xinp(k)*xno(k)+alp1(k)*xinop(k) | *xne(k)*(1.-falp12d)+2.*alp3(k)*xin2p(k)*xne(k)* | (1.-falp32d)+alp3(k)*xin2p(k)*xne(k)*falp32d | +rb4(k)*xn2d(k)*xno(k)+rb5(k)*xn2d(k)*xne(k) | +rb7(k)*xn2d(k)+dnosrb(k)*xnno(k) xln4s(k)=rk4(k)*xio2p(k)+rb1(k)*xno2(k)+rb3(k)*xnno(k) | +rb8(k)*xnoh(k)+rb13(k)*xnno2(k) ! ! no production and loss: pno(k)=rb1(k)*xno2(k)*xn4s(k)+rb2(k)*xno2(k)*xn2d(k) | +rb8(k)*xnoh(k)*xn4s(k)+rb11(k)*xno(k)*xnno2(k) | +dno2t(k)*xnno2(k) xlno(k)=rk5(k)*xio2p(k)+rb3(k)*xn4s(k)+rb6(k)*xn2d(k) | +dnosrb(k)+dinolya(k)+rb9(k)*xno3(k)+rb10(k)*xnho2(k) | +dnoeuv(k) ! ! no2 production and loss: pno2(k)=rb9(k)*xno3(k)*xnno(k)+rb10(k)*xnho2(k)*xnno(k) xlno2(k)=rb11(k)*xno(k)+dno2t(k)+rb12(k)*xno3(k) | +rb13(k)*xn4s(k) ! ! nox=no+no2 partioning: ratio4=no2/no, ratio5=no/nox ! ratio4(k)=(rb9(k)*xno3(k)+rb10(k)*xnho2(k))/ | (rb11(k)*xno(k)+dno2t(k)+rb12(k)*xno3(k)+ | rb13(k)*xn4s(k)) ratio5(k)=1./(1.+ratio4(k)) ! ! nox production and loss: pnox(k)=rb1(k)*xno2(k)*xn4s(k)+rb2(k)*xno2(k)*xn2d(k)+ | rb8(k)*xnoh(k)*xn4s(k) xlnoxt(k)=(rk5(k)*xio2p(k)+rb3(k)*xn4s(k)+rb6(k)*xn2d(k)+ | dnosrb(k)+dinolya(k)+dnoeuv(k))*xnno(k)+ | (rb12(k)*xno3(k)+rb13(k)*xn4s(k))*xnno2(k) xlnox(k)=(rk5(k)*xio2p(k)+rb3(k)*xn4s(k)+rb6(k)*xn2d(k)+ | dnosrb(k)+dinolya(k)+dnoeuv(k))*ratio5(k)+ | (rb12(k)*xno3(k)+rb13(k)*xn4s(k))*ratio5(k)* | ratio4(k) ! ! Production and loss of o3: ! po3(k)=(rkm21(k)*xno2(k)*xno(k)+rkm22(k)*xno2(k)*xno2(k)+ | rkm23(k)*xno2(k)*xnn2(k))*xno(k) xlo3(k)=(do3t(k)+(rkm7a(k)+rkm7b(k))*xno1d(k)+ | rkm24(k)*xno(k)+rkm29(k)*xnoh(k)+rkm34(k)*xnho2(k)+ | rkm37(k)*xnh(k)+del1(k)*xncl(k)+rb9(k)*xnno(k)+ | rb12(k)*xnno2(k))*xno3(k) ! ! Production and loss of o: ! po(k)=2.*do2t(k)*xno2(k)+do3t(k)*xno3(k)+dh2ot(k)* | xnh2o(k)*eps+rkm30(k)*xnoh(k)*xnoh(k)+rkm40(k)* | xnh(k)*xnho2(k)+rb1(k)*xn4s(k)*xno2(k) | +rb2(k)*xn2d(k)*xno2(k)+rb3(k)*xn4s(k)*xnno(k) | +rb6(k)*xn2d(k)*xnno(k)+dnosrb(k)*xnno(k) | +dno2t(k)*xnno2(k)+rb13(k)*xn4s(k)*xnno2(k) | +rk1(k)*xiop(k)*xno2(k)+rk4(k)*xio2p(k)*xn4s(k) | +rk7(k)*xinp(k)*xno2(k)+rk10(k)*xiop(k)*xn2d(k) | +rk12(k)*xiop(k)*xnh(k)+rk15(k)*xiop(k)*xnh2o(k) | +rk16(k)*xiop2p(k)*xnn2(k)+rk23(k)*xiop2p(k)*xno2(k) | +rk24(k)*xiop2d(k)*xnn2(k)+rk27(k)*xiop2d(k)*xno2(k) | +alp1(k)*xinop(k)*xne(k)+2.*alp2(k)*xio2p(k)*xne(k) | +rkm7b(k)*xno3(k)*xno1d(k) xlo(k)=rkm3(k)*xnh2o(k)*xno1d(k)+rkm4(k)*xnh2(k)*xno1d(k) | +(rkm5a(k)+rkm5b(k))*xnch4(k)*xno1d(k) | +rkm7a(k)*xno3(k)*xno1d(k)+2.*rkm20(k)*xm*xno(k)*xno(k) | +rkm21(k)*xno2(k)*xno(k)*xno(k) | +rkm22(k)*xno2(k)*xno2(k)*xno(k) | +rkm23(k)*xno2(k)*xnn2(k)*xno(k) | +rkm24(k)*xno3(k)*xno(k)+rkm25(k)*xnoh(k)*xno(k) | +rkm26(k)*xnho2(k)*xno(k)+rkm27(k)*xnh2o2(k)*xno(k) | +rkm28(k)*xnh2(k)*xno(k)+gam4(k)*ch3(k)*xno(k) | +gam11(k)*ch2o(k)*xno(k)+gam13(k)*xnch4(k)*xno(k) | +gam14(k)*xnco(k)*xm*xno(k)+gam15(k)*xnco2(k)*xno1d(k) | +del2(k)*xnclo(k)*xno(k)+rb11(k)*xnno2(k)*xno(k) | +rk3(k)*xin2p(k)*xno(k)+rk8(k)*xinp(k)*xno(k) | +rk11(k)*xihp(k)*xno(k) ! ! ox=o+o3 partitioning ratio6=o3/o, ratio7=o/ox ! ratio6(k)=rkm21(k)*xno2(k)*xm/(do3t(k)+(rkm7a(k)+rkm7b(k)) | *xno1d(k)+rkm24(k)*xno(k)+rkm29(k)*xnoh(k)+rkm34(k) | *xnho2(k)+rkm37(k)*xnh(k)+del1(k)*xncl(k)+rb9(k)* | xnno(k)+rb12(k)*xnno2(k)) ratio7(k)=1./(1.+ratio6(k)) ! ! ox production and loss: a4(k)=eps*dh2olya(k)*xnh2o(k)+rkm30(k)*xnoh(k)*xnoh(k) | +rkm40(k)*xnh(k)*xnho2(k)+rb3(k)*xn4s(k)*xnno(k) | +rb6(k)*xn2d(k)*xnno(k)+dnosrb(k)*xnno(k) | +dno2t(k)*xnno2(k)+rb13(k)*xn4s(k)*xnno2(k) | +rkm7b(k)*xno3(k)*xno1d(k)+rk4(k)*xio2p(k)*xn4s(k) | +rk10(k)*xiop(k)*xn2d(k)+rk12(k)*xiop(k)*xnh(k) | +rk15(k)*xiop(k)*xnh2o(k)+rk16(k)*xiop2p(k)*xnn2(k) | +rk24(k)*xiop2d(k)*xnn2(k)+alp1(k)*xinop(k)*xne(k) | +2.*alp2(k)*xio2p(k)*xne(k) a5(k)=2.*(do2src(k)+do2srb(k)+do2lya(k)+do2hz(k)+do2euv(k)) | +rb1(k)*xn4s(k)+rb2(k)*xn2d(k)+rk1(k)*xiop(k) | +rk7(k)*xinp(k)+rk23(k)*xiop2p(k)+rk27(k)*xiop2d(k) | +dio2p(k) b1(k)=(rkm3(k)*xnh2o(k)+rkm4(k)*xnh2(k)+(rkm5a(k)+rkm5b(k)) | *xnch4(k)+gam15(k)*xnco2(k))*xno1d(k) b2(k)=2.*rkm24(k)*ratio6(k)*ratio7(k)**2 b3(k)=(rkm25(k)*xnoh(k)+rkm26(k)*xnho2(k)+rkm27(k)*xnh2o2(k) | +rkm28(k)*xnh2(k)+gam4(k)*ch3(k)+gam11(k)*ch2o(k) | +gam13(k)*xnch4(k)+gam14(k)*xnco(k)*xm+del2(k)*xnclo(k) | +rb11(k)*xnno2(k)+rk3(k)*xin2p(k)+rk8(k)*xinp(k) | +rk11(k)*xihp(k)+diop(k)+((rkm7a(k)+rkm7b(k))*xno1d(k) | +rkm29(k)*xnoh(k)+rkm34(k)*xnho2(k)+rkm37(k)*xnh(k) | +del1(k)*xncl(k)+rb9(k)*xnno(k)+rb12(k)*xnno2(k) | +rkm7a(k)*xno1d(k))*ratio6(k))*ratio7(k) c1(k)=2.*rkm24(k)*ratio6(k)*ratio7(k)**2 c2(k)=(rkm25(k)*xnoh(k)+rkm26(k)*xnho2(k)+del2(k)*xnclo(k) | +rb11(k)*xnno2(k)+(2.*rkm7a(k)*xno1d(k)+rkm7b(k)*xno1d(k) | +rkm29(k)*xnoh(k)+2.*rkm34(k)*xnho2(k)+rkm37(k)*xnh(k) | +del1(k)*xncl(k)+rb9(k)*xnno(k)+rb12(k)*xnno2(k)) | *ratio6(k))*ratio7(k) c3(k)=rkm31(k)*xnoh(k)*xnho2(k)+rkm35(k)*xnho2(k)*xnho2(k) | +rkm38(k)*xnh(k)*xnho2(k)+gam6(k)*ch3o2(k)*xnho2(k) | +gam7(k)*ch3o2(k)*ch3o2(k)+gam15(k)*xnco(k)*xno1d(k) | +rk5(k)*xio2p(k)*xnno(k) d1(k)=rkm36(k)*xnh(k)*xm+gam3(k)*ch3(k)*xm+gam9(k)*ch3o(k) | +gam12(k)*hco(k)+rb1(k)*xn4s(k)+rb2(k)*xn2d(k) | +rk1(k)*xiop4s(k)+(rk6(k)+rk7(k))*xinp(k) | +rk9(k)*xin2p(k)+rk23(k)*xiop2p(k)+rk27(k)*xiop2d(k) | +do2t(k)+dio2p(k) xnox(k)=xno(k)+xno3(k) ! ! fsxx are module data: fs11(k)=-d1(k) fs12(k)=c2(k)+c1(k)*xnox(k) fs21(k)=a5(k) fs22(k)=-b3(k)-b2(k)*xnox(k) fs1(k)=c3(k) fs2(k)=a4(k)-b1(k) ! write(6,"('chem: k=',i3,' fs11,12,21,22=',4e12.4)") ! | k,fs11(k),fs12(k),fs21(k),fs22(k) ! xnhox(k)=xnh(k)+xnoh(k)+xnho2(k) xnnoz(k)=xnno(k)+xnno2(k) rno2no(k)=ratio4(k) rnonoz(k)=ratio5(k) ho2hr(k) =ratio2(k) ohhr(k) =ratio1(k) hhoxr(k) =ratio3(k) o3or(k) =ratio6(k) ooxr(k) =ratio7(k) ! wl1d = 12700. wl1s = 7620. ene1d=12400.*erg/wl1d ene1s=12400.*erg/wl1s xno21dc(k)=adl*xno21d(k)*ene1d/rho(k) xno21sc(k)=asg*xno21s(k)*ene1s/rho(k) xmetah(k)=0. ! ! Neutral sodium production: ! (see old compo.F for alternative values) as=1.0e-2 zps= 89. h1s=5. #ifdef EXPO pna(k)=as*expo(-((zpht(k)-zps)/h1s)**2,0) #else pna(k)=as*exp(-((zpht(k)-zps)/h1s)**2) #endif pnas(k)=pna(k) enddo ! k=1,nlev (was do 1) ! ! Sodium, continued: sumnap = 0. do k=1,nlev sumnap=sumnap+(pna(k)+pna(k-1))*0.5*(zpht(k)-zpht(k-1))*1.e+5 enddo ! k=1,nlev call addfhist(xnhox,'HOX','HOX', | 'HOX (H+OH+HO2)','CM-3',nlev,1,1) call addfhist(xnox,'OX','OX', | 'OX (O+O3)','CM-3',nlev,1,1) ! call fprint8('chem: xnhox, xnox, xnnoz',nlev, ! | xnhox , xnox , xnnoz , phoxic , dum,dum,dum,dum, ! | 'xnhox','xnox','xnnoz','phoxic',' ',' ',' ',' ') ! ! call fprint8('chem: metastable n(2d)',nlev, ! | qn2d , qn4s , dn2 , pn2d , xln2d , xn2d , dum,dum, ! | 'qn2d','qn4s','dn2','pn2d','xln2d','xn2d',' ',' ') ! ! call fprint8('chem: metastable o2(1sg)',nlev, ! | po21sg , xlo21sg , xno21s , dum,dum,dum,dum,dum, ! | 'po21sg','xlo21sg','xno21s',' ',' ',' ',' ',' ') ! ! call fprint8('chem: metastable o2(1dl)',nlev, ! | po21dl , xlo21dl , xno21d ,dum,dum,dum,dum,dum, ! | 'po21dl','xlo21dl','xno21d',' ',' ',' ',' ',' ') ! ! call fprint8('chem: metastable o(1d)',nlev, ! | po1d , xlo1d , xno1d ,dum,dum,dum,dum,dum, ! | 'po1d','xlo1d','xno1d',' ',' ',' ',' ',' ') ! ! call fprint8('chem: methane oxidation (ch3,ch3o2)',nlev, ! | pch3 , xlch3 , ch3 , pch3o2 , xlch3o2 , ch3o2 ,dum,dum, ! | 'pch3','xlch3','ch3','pch3o2','xlch3o2','ch3o2',' ',' ') ! ! call fprint8('chem: methane oxidation (ch3ooh,ch3o)',nlev, ! | pch3oo , xlch3oo , ch3ooh , pch3o , xlch3o , ch3o ,dum,dum, ! | 'pch3oo','xlch3oo','ch3ooh','pch3o','xlch3o','ch3o',' ',' ') ! ! call fprint8('chem: methane oxidation (ch2o,hco)',nlev, ! | pch2o , xlch2o , ch2o , phco , xlhco , hco ,dum,dum, ! | 'pch2o','xlch2o','ch2o','phco','xlhco','hco',' ',' ') ! ! call fprint8('chem: h2o2,ch4,co2',nlev, ! | ph2o2 , xlh2o2 , xnh2o2 , pch4 , xlch4 , pco2 , xlco2,dum, ! | 'ph2o2','xlh2o2','xnh2o2','pch4','xlch4','pco2','xlco2',' ') ! ! call fprint8('chem: production and loss: co, h2o, h2',nlev, ! | pco , xlco , ph2o , xlh2o , ph2 , xlh2 , dum,dum, ! | 'pco','xlco','ph2o','xlh2o','ph2','xlh2',' ',' ') ! ! call fprint8('chem: production and loss: ho2,oh,h',nlev, ! | pho2 , xlho2 , poh , xloh , ph , xlh , dum,dum, ! | 'pho2','xlho2','poh','xloh','ph','xlh',' ',' ') ! ! call fprint8('chem: prod, loss, and partitioning of hox:',nlev, ! | phox , xlhoxt , ratio1 , ratio2 , ratio3 , ! | xlhox , dum,dum, ! | 'phox','xlhoxt','rat1:oh/h','rat2:ho2/h','rat3:h/hox', ! | 'xlhox',' ',' ') ! ! call fprint8('chem: production and loss: n4s,no,no2',nlev, ! | pn4s , xln4s , pno , xlno , pno2 , xlno2 , dum,dum, ! | 'pn4s','xln4s','pno','xlno','pno2','xlno2',' ',' ') ! ! call fprint8('chem: nox partitioning, production and loss',nlev, ! | ratio4 , ratio5 , pnox , xlnoxt , xlnox , dum,dum,dum, ! | 'ratio4','ratio5','pnox','xlnoxt','xlnox',' ',' ',' ') ! ! call fprint8('chem: o3,o production and loss',nlev, ! | po3 , xlo3 , po , xlo , dum,dum,dum,dum, ! | 'po3','xlo3','po','xlo',' ',' ',' ',' ') ! ! call fprint8('chem: ox production and loss',nlev, ! | ratio6 , ratio7 , a4 , a5 , b1 , b2 , b3 , dum, ! | 'ratio6','ratio7','a4','a5','b1','b2','b3',' ') ! ! call fprint8('chem: ox production and loss, cont',nlev, ! | c1 , c2 , c3 , d1 , xnox , dum , dum , dum, ! | 'c1','c2','c3','d1','xnox',' ',' ',' ') ! ! Photochemical equilibrium at lower boundary xm=xnn2(1)+xno2(1)+xno(1)+xno3(1) bb2=b2(1)+2.*rkm20(1)*xm*ratio7(1)**2 xnox(1)=(-b3(1)+sqrt(b3(1)**2-4.*bb2*b1(1)))/(2.*bb2) ! write(6,"('chem: xnox(1)=',e12.4)") xnox(1) ! ! Save xno and xno3 lbc for use in sub face: xno_lb = xno(1) xno3_lb = xno3(1) ! write(6,"('chem: xno_lb,xno3_lb=',2e12.4)") xno_lb,xno3_lb ! xno(1)=xnox(1)*ooxr(1) xno3(1)=o3or(1)*xno(1) ! write(6,"('chem: xno(1),xno3(1)=',2e12.4)") xno(1),xno3(1) ! ! ps2b is module data, and will be fb(2) ps2b=xm*xo3lb*48./(xno2(1)*32.+xnn2(1)*28.) aa=(2.*rkm31(1)*ratio1(1)*ratio2(1)+2.*(rkm38(1)+rkm40(1))* | ratio2(1)+2.*rkm35(1)*ratio2(1)**2+2.*rkm30(1)* | ratio1(1)**2+2.*rkm41(1)*xm)*ratio3(1)**2 bb=(gam6(1)*ch3o2(1)*ratio2(1)+(gam1(1)*xnch4(1)+gam8(1)* | ch3ooh(1)+gam10(1)*ch2o(1))*ratio1(1)+rk12(1)*xiop(1))* | ratio3(1) cc=-phox(1) xn4slb=pn4s(1)/xln4s(1) ! ! Set hox top boundary: call hoxubc(xnh(nlev),tn(nlev),gz(nlev),zpht(nlev)) ! ! Sub face calls composition routines: call face(ut,dt,dzp,zpb) do k=1,nlev xnn2(k) = max(xnn2(k) ,1.e-20) xnat(k) = max(xnat(k) ,1.e-20) xnh2o(k) = max(xnh2o(k) ,1.e-20) xnhox(k) = max(xnhox(k) ,1.e-20) xnh2o2(k) = max(xnh2o2(k),1.e-60) xno1d(k) = max(xno1d(k) ,1.e-20) xnh2(k) = max(xnh2(k) ,1.e-20) xnoh(k) = max(xnoh(k) ,1.e-20) xnho2(k) = max(xnho2(k) ,1.e-20) xno3(k) = max(xno3(k) ,1.e-20) xnh(k) = max(xnh(k) ,1.e-20) xnch4(k) = max(xnch4(k) ,1.e-20) xnco(k) = max(xnco(k) ,1.e-20) xnco2(k) = max(xnco2(k) ,1.e-20) xnno2(k) = max(xnno2(k) ,1.e-20) xnno(k) = max(xnno(k) ,1.e-20) xn2d(k) = max(xn2d(k) ,1.e-20) xn4s(k) = max(xn4s(k) ,1.e-20) xnnoz(k) = max(xnnoz(k) ,1.e-20) xnhtot(k) = 2.*xnh2o(k)+2.*xnh2o2(k)+2.*xnh2(k)+xnoh(k) | +xnho2(k)+xnh(k)+4.*xnch4(k) enddo ! k=1,nlev ! write(6,"(/,'chem: istep=',i5,' of ',i5,': xnhtot=',/,(6e12.4))") ! | istep,nstep,xnhtot ! write(6,"('chem: istep=',i5,': xnh2o=',/,(6e12.4))") istep,xnh2o ! write(6,"('chem: istep=',i5,': xnh2o2=',/,(6e12.4))")istep,xnh2o2 ! write(6,"('chem: istep=',i5,': xnh2=',/,(6e12.4))") istep,xnh2 ! write(6,"('chem: istep=',i5,': xnoh=',/,(6e12.4))") istep,xnoh ! write(6,"('chem: istep=',i5,': xnho2=',/,(6e12.4))") istep,xnho2 ! write(6,"('chem: istep=',i5,': xnh=',/,(6e12.4))") istep,xnh ! write(6,"('chem: istep=',i5,': xnch4=',/,(6e12.4))") istep,xnch4 ! call print_fhist('XNHTOT',f_hist,nf_hist,'end sub chem (comp.F)') ! write(char80,"('chem: updated number densities: iday=',i3, ! | ' ut=',f5.2)") iday,ut ! call fprint8(char80,nlev, ! | xno2 , xno , xnn2 ,dum,dum,dum,dum,dum, ! | 'xno2','xno','xnn2',' ' ,' ' , ' ',' ' ,' ') ! call fprint8('chem: updated number densities',nlev, ! | xno2 , xnox , xno , xno3 , xn2d , xn4s , xnnoz ,dum, ! | 'xno2' ,'xnox' ,'xno','xno3','xn2d', 'xn4s','xnnoz',' ') ! call fprint8('chem: updated number densities',nlev, ! | xnno , xnno2 , xnarg , xnhe , xnch4 ,xnh2 ,xnco2 ,dum, ! | 'xnno','xnno2','xnarg','xnhe','xnch4','xnh2','xnco2',' ') ! call fprint8('chem: updated number densities',nlev, ! | xnco , xnhox , xnh , xnho2 , xnoh , xnh2o , xnat ,dum, ! | 'xnco','xnhox','xnh','xnho2','xnoh','xnh2o','xnat',' ') ! call fprint8('chem: updated number densities',nlev, ! | xnn2 , xnh2o2 , dum , dum , dum , dum , dum ,dum, ! | 'xnn2','xnh2o2',' ',' ',' ',' ',' ',' ') call sodium end subroutine chemistry !----------------------------------------------------------------------- subroutine face(ut,step,dzp,zpb) use flds_atmos,only: tn,xno2,xno,xnarg,xnhe,xnnoz,xn4s,xnch4, | xnh2,xnco2,xnco,xnhox,xnh2o,xnn2,xnno2,xnno,xno3,xnox,xn2d, | xnh,xnho2,xnoh use flds_sodium,only: xnat use flds_fields,only: ooxr,rnonoz,ho2hr,ohhr,hhoxr,o3or,rno2no use cons,only: rmo2,rmo,rmn2,rmar,rmnat,rmhe,rmno,rmn4s,rmch4, | rmh2,rmco,rmco2,rmh2o,rmhox,p0,boltz ! use input,only: istep,nstep implicit none ! ! Args: real,intent(in) :: ut,step,dzp,zpb ! ! Local: integer :: k integer,save :: ncalls=0 real,dimension(nlev) :: rhoinv,frj,ftn,fw character(len=80) :: char80 real,parameter :: | rmn4sinv = 1./rmn4s, | rmnoinv = 1./rmno, | rmarinv = 1./rmar, | rmheinv = 1./rmhe, | rmch4inv = 1./rmch4, | rmh2inv = 1./rmh2, | rmco2inv = 1./rmco2, | rmcoinv = 1./rmco, | rmhoxinv = 1./rmhox, | rmh2oinv = 1./rmh2o, | rmnatinv = 1./rmnat ! ! Comp routines outputs (mmr): real :: xo2lb,xolb,xn2lb real :: dum,s,expds,factor ! ! pxxxb is output by sub minor real :: pn4sb,pnozb,parb,pheb,pch4b,ph2b,pco2b,pcob, | phoxb,ph2ob,pnatb ! ncalls = ncalls+1 ! ! Save number densities from field modules: ! (whole-array operations) ! xn_o2 = xno2 xn_o = xno+xno3 xn_n2 = xnn2 xn_ar = xnarg xn_he = xnhe xn_noz = xnnoz xn_n4s = xn4s xn_ch4 = xnch4 xn_h2 = xnh2 xn_co2 = xnco2 xn_co = xnco xn_hox = xnhox xn_h2o = xnh2o xn_nat = xnat ! ! Save o2,o,n2 number density lbc before update for sub photo: xo2lb = xn_o2(1) xolb = xn_o (1) xn2lb = xn_n2(1) ! ! In the old code, psi species (f-array) are calculated only ! in the first time step (see ISTAR in compo.F) -- why? This ! means the f-array psi species are always the same -- after ! the first step, they are not updated to use new values from ! the current step. ! if (ncalls==1) then ! write(6,"('begin face: istep=',i4,' xnh2o=',/,(6e12.4))") ! | istep,xnh2o ! ! Convert from number density to mass mixing ratio: ! (was DO 3 in old compo.F) do k=1,nlev rhoinv(k) = 1./(xno2(k)*rmo2 + xno(k)*rmo + xnn2(k)*rmn2) ps_o2(k) = xno2(k) *rmo2 *rhoinv(k) ps_ox(k) = (xno(k)+xno3(k))*rmo*rhoinv(k) ps_n2(k) = (1.-ps_o2(k)-ps_ox(k)) ps_ar(k) = xnarg(k) *rmar *rhoinv(k) ps_he(k) = xnhe(k) *rmhe *rhoinv(k) ps_noz(k) = (xnno(k)+xnno2(k))*rmno*rhoinv(k) ps_n4s(k) = xn4s(k) *rmn4s*rhoinv(k) ps_ch4(k) = xnch4(k) *rmch4*rhoinv(k) ps_h2(k) = xnh2(k) *rmh2 *rhoinv(k) ps_co2(k) = xnco2(k) *rmco2*rhoinv(k) ps_co(k) = xnco(k) *rmco *rhoinv(k) ps_hox(k) = xnhox(k) *rmhox*rhoinv(k) ps_h2o(k) = xnh2o(k) *rmh2o*rhoinv(k) ps_nat(k) = xnat(k) *rmnat*rhoinv(k) enddo ! k=1,nlev ! write(6,"('face: istep=',i4,' rhoinv=',/,(6e12.4))") ! | istep,rhoinv ! write(6,"('face: istep=',i4,' ps_h2o=',/,(6e12.4))") ! | istep,ps_h2o ps_ox(1) = (xno_lb+xno3_lb)*rmo*rhoinv(1) ! ! Convert to half-levels (was DO 4): do k=1,nlev-1 ps_o2 (k) = 0.5*(ps_o2 (k)+ps_o2 (k+1)) ps_ox (k) = 0.5*(ps_ox (k)+ps_ox (k+1)) ps_n2 (k) = 0.5*(ps_n2 (k)+ps_n2 (k+1)) ps_ar (k) = 0.5*(ps_ar (k)+ps_ar (k+1)) ps_he (k) = 0.5*(ps_he (k)+ps_he (k+1)) ps_noz(k) = 0.5*(ps_noz(k)+ps_noz(k+1)) ps_n4s(k) = 0.5*(ps_n4s(k)+ps_n4s(k+1)) ps_ch4(k) = 0.5*(ps_ch4(k)+ps_ch4(k+1)) ps_h2 (k) = 0.5*(ps_h2 (k)+ps_h2 (k+1)) ps_co2(k) = 0.5*(ps_co2(k)+ps_co2(k+1)) ps_co (k) = 0.5*(ps_co (k)+ps_co (k+1)) ps_hox(k) = 0.5*(ps_hox(k)+ps_hox(k+1)) ps_h2o(k) = 0.5*(ps_h2o(k)+ps_h2o(k+1)) ps_nat(k) = 0.5*(ps_nat(k)+ps_nat(k+1)) enddo ! k=1,nlev-1 ! ! Upper boundaries: ps_o2(nlev) = 1.5*(xno2(nlev) *rmo2*rhoinv(nlev))- | 0.5*(xno2(nlev-1) *rmo2*rhoinv(nlev-1)) ps_ox(nlev) = | 1.5*((xno(nlev) +xno3(nlev)) *rmo *rhoinv(nlev))- | 0.5*((xno(nlev-1)+xno3(nlev-1))*rmo *rhoinv(nlev-1)) ps_n2(nlev) = (1.-ps_o2(nlev)-ps_ox(nlev)) ps_ar(nlev) = 1.5*(xnarg(nlev) *rmar*rhoinv(nlev))- | 0.5*(xnarg(nlev-1)*rmar*rhoinv(nlev-1)) ps_he(nlev) = 1.5*(xnhe(nlev) *rmhe*rhoinv(nlev))- | 0.5*(xnhe(nlev-1) *rmhe*rhoinv(nlev-1)) ps_noz(nlev) = | 1.5*((xnno(nlev) +xnno2(nlev)) *rmno*rhoinv(nlev))- | 0.5*((xnno(nlev-1)+xnno2(nlev-1))*rmno*rhoinv(nlev-1)) ps_n4s(nlev) = 1.5*(xn4s(nlev) *rmn4s*rhoinv(nlev))- | 0.5*(xn4s(nlev-1) *rmn4s*rhoinv(nlev-1)) ps_ch4(nlev) = 1.5*(xnch4(nlev) *rmch4*rhoinv(nlev))- | 0.5*(xnch4(nlev-1)*rmch4*rhoinv(nlev-1)) ps_h2(nlev) = 1.5*(xnh2(nlev) *rmh2 *rhoinv(nlev))- | 0.5*(xnh2(nlev-1) *rmh2 *rhoinv(nlev-1)) ps_co2(nlev) = 1.5*(xnco2(nlev) *rmco2*rhoinv(nlev))- | 0.5*(xnco2(nlev-1)*rmco2*rhoinv(nlev-1)) ps_co(nlev) = 1.5*(xnco(nlev) *rmco *rhoinv(nlev))- | 0.5*(xnco(nlev-1) *rmco *rhoinv(nlev-1)) ps_hox(nlev) = 1.5*(xnhox(nlev) *rmhox*rhoinv(nlev))- | 0.5*(xnhox(nlev-1)*rmhox*rhoinv(nlev-1)) ps_h2o(nlev) = 1.5*(xnh2o(nlev) *rmh2o*rhoinv(nlev))- | 0.5*(xnh2o(nlev-1)*rmh2o*rhoinv(nlev-1)) ps_nat(nlev) = 1.5*(xnat(nlev) *rmnat*rhoinv(nlev))- | 0.5*(xnat(nlev-1) *rmnat*rhoinv(nlev-1)) ! write(6,"('face: istep=',i4,' ps_h2o=',/,(6e12.4))") ! | istep,ps_h2o ! ! write(char80,"('face: ut=',f5.2,' psi species')") ut ! call fprint8(char80,nlev, ! | ps_o2 , ps_ox, ps_n2 , dum , dum , dum , dum ,dum, ! | 'ps_o2','ps_ox','ps_n2', ' ' ,' ' ,' ' ,' ' ,' ') ! ! call fprint8(char80,nlev, ! | ps_o2 , ps_ox, ps_noz , ps_n4s , ps_ch4 , ps_h2 , ps_co2 ,dum, ! | 'ps_o2','ps_ox','ps_noz','ps_n4s','ps_ch4','ps_h2','ps_co2',' ') ! ! call fprint8('face: psi species',nlev, ! | ps_o2 , ps_ox, ps_n2 , ps_ar , ps_he , ps_noz , ps_n4s ,dum, ! | 'ps_o2','ps_ox','ps_n2','ps_ar','ps_he','ps_noz','ps_n4s',' ') ! call fprint8('face: psi species, cont',nlev, ! | ps_ch4 , ps_h2 , ps_co2 , ps_co , ps_hox , ps_h2o , ps_nat,dum, ! | 'ps_ch4','ps_h2','ps_co2','ps_co','ps_hox','ps_h2o','ps_nat',' ') endif ! ncalls==1 ! From sub comp: expdz = exp(-.5*dzp) ! was C(86) (dzp was C(3)) expdzinv = 1./expdz ! was C(87) s = zpb+.5*dzp exps(1) = exp(-s) expds = exp(-dzp) do k=2,nlev exps(k) = exps(k-1)*expds enddo ! k=2,nlev do k=1,nlev-1 ftn(k) = 0.5*(tn(k)+tn(k+1)) ! was f(i,ntk) enddo ! k=1,nlev-1 ftn(nlev) = tn(1) ! bottom boundary in top slot?? (change later) do k=1,nlev fw(k) = 0. frj(k) = 0. enddo ! k=1,nlev ! ! New JPL-85 rates (formerly sub rates in compo.F): do k=1,nlev rk12_comp(k) = ooxr(k)*ooxr(k)*9.59e-34*exp(480./tn(k)) rk13_comp(k) = 0. rk14_comp(k) = 0. rk15_comp(k) = 0. enddo ! k=1,nlev ! ! Edy is in module flds_fields, and was calculated by sub init_edy ! in cons.F. difk(:) = edy(:) ! ! Call composition routines (mmr): call comp(ut,step,ftn,fw,ps_o2,ps_ox) call photo(xo2lb,xolb,xn2lb) call compn2d(xn2d) ! ! Minor species are returned by subs compxxx: ps_xxx(nlev) is inout, ! and pxxxb boundary scalar is output. ! call compn4s(ps_n4s ,ftn,fw,pn4sb) call compnoz(ps_noz ,ftn,fw,pnozb) call compar (ps_ar ,ftn,fw,parb) call comphe (ps_he ,ftn,fw,pheb) call compch4(ps_ch4 ,ftn,fw,pch4b) call comph2 (ps_h2 ,ftn,fw,ph2b) call compco2(ps_co2 ,ftn,fw,pco2b) call compco (ps_co ,ftn,fw,pcob) call comphox(ps_hox ,ftn,fw,phoxb) ! write(6,"('before comph2o: istep=',i4,' ps_h2o=',/,(6e12.4))") ! | istep,ps_h2o call comph2o(ps_h2o ,ftn,fw,ph2ob) ! write(6,"('after comph2o: istep=',i4,' ps_h2o=',/,(6e12.4))") ! | istep,ps_h2o call compnat(ps_nat ,ftn,fw,pnatb) ! ! Convert updated species back to number density: do k=1,nlev-1 ! was DO 7 xno2 (k+1) = .5*(ps_o2 (k)+ps_o2 (k+1)) xnox (k+1) = .5*(ps_ox (k)+ps_ox (k+1)) xn4s (k+1) = .5*(ps_n4s(k)+ps_n4s(k+1)) xnnoz(k+1) = .5*(ps_noz(k)+ps_noz(k+1)) xnarg(k+1) = .5*(ps_ar (k)+ps_ar (k+1)) xnhe (k+1) = .5*(ps_he (k)+ps_he (k+1)) xnch4(k+1) = .5*(ps_ch4(k)+ps_ch4(k+1)) xnh2 (k+1) = .5*(ps_h2 (k)+ps_h2 (k+1)) xnco2(k+1) = .5*(ps_co2(k)+ps_co2(k+1)) xnco (k+1) = .5*(ps_co (k)+ps_co (k+1)) xnhox(k+1) = .5*(ps_hox(k)+ps_hox(k+1)) xnh2o(k+1) = .5*(ps_h2o(k)+ps_h2o(k+1)) xnat(k+1) = .5*(ps_nat(k)+ps_nat(k+1)) enddo ! k=1,nlev ! ! Lower boundary: xno2(1) = ps1b xnox(1) = ps2b xn4s(1) = .5*(ps_n4s(1)+pn4sb) xnnoz(1) = .5*(ps_noz(1)+pnozb) xnarg(1) = .5*(ps_ar (1)+parb) xnhe (1) = .5*(ps_he (1)+pheb) xnch4(1) = .5*(ps_ch4(1)+pch4b) xnh2 (1) = .5*(ps_h2 (1)+ph2b) xnco2(1) = .5*(ps_co2(1)+pco2b) xnco (1) = .5*(ps_co (1)+pcob) xnhox(1) = .5*(ps_hox(1)+phoxb) xnh2o(1) = .5*(ps_h2o(1)+ph2ob) xnat(1) = .5*(ps_nat(1)+pnatb) ! ! Convert to number density: do k=1,nlev ! was DO 8 factor=p0*exps(k)*expdzinv/(boltz*tn(k)*(xno2(k)*rmo2inv+ | xnox(k)*rmoinv+(1.-xno2(k)-xnox(k))*rmn2inv)) xnn2(k) = (1.-xno2(k)-xnox(k))*factor*rmn2inv xno2(k) = xno2 (k)*factor*rmo2inv xnox(k) = xnox (k)*factor*rmoinv xno(k) = ooxr(k)*xnox(k) xno3(k) = o3or(k)*xno(k) xn4s(k) = xn4s (k)*factor*rmn4sinv xnnoz(k) = xnnoz(k)*factor*rmnoinv xnno(k) = xnnoz(k)*rnonoz(k) xnno2(k) = xnno(k)*rno2no(k) xnarg(k) = xnarg(k)*factor*rmarinv xnhe(k) = xnhe (k)*factor*rmheinv xnch4(k) = xnch4(k)*factor*rmch4inv xnh2 (k) = xnh2 (k)*factor*rmh2inv xnco2(k) = xnco2(k)*factor*rmco2inv xnco (k) = xnco (k)*factor*rmcoinv xnhox(k) = xnhox(k)*factor*rmhoxinv xnh(k) = xnhox(k)*hhoxr(k) xnho2(k) = xnh(k)*ho2hr(k) xnoh(k) = xnh(k)*ohhr(k) xnh2o(k) = xnh2o(k)*factor*rmh2oinv xnat(k) = xnat(k)*factor*rmnatinv enddo ! k=1,nlev ! ! if (istep==nstep) then ! write(char80,"('face: istep=',i4,' number densities:')") istep ! call fprint8(char80,nlev, ! | xno2 , xnox , xno , xno3 , xn2d , xn4s , xnnoz ,dum, ! | 'xno2' ,'xnox' ,'xno','xno3','xn2d', 'xn4s','xnnoz',' ') ! call fprint8(char80,nlev, ! | xnno , xnno2 , xnarg , xnhe , xnch4 ,xnh2 ,xnco2 ,dum, ! | 'xnno','xnno2','xnarg','xnhe','xnch4','xnh2','xnco2',' ') ! call fprint8(char80,nlev, ! | xnco , xnhox , xnh , xnho2 , xnoh , xnh2o , xnat ,dum, ! | 'xnco','xnhox','xnh','xnho2','xnoh','xnh2o','xnat',' ') ! call fprint8(char80,nlev, ! | xnn2 , dum , dum , dum , dum , dum , dum ,dum, ! | 'xnn2',' ',' ',' ',' ',' ',' ',' ') ! endif end subroutine face !----------------------------------------------------------------------- subroutine comp(ut,step,ftn,fw,o2out,oxout) use cons,only: rmo2,rmo,rmn2,p0,boltz implicit none ! ! Args: real,intent(in) :: ut,step real,dimension(nlev),intent(in) :: ftn,fw real,dimension(nlev),intent(out) :: o2out,oxout ! ! Local: integer :: k,kk,km,kp,ktmp,l,ll,m real,dimension(nlev) :: embar real :: ps0(2),ep(2,2),ak(2,2,2), | delta(2,2),gama(nlev,2,2),wkm1(2,2),wkm2(2,2), | wkv1(2),pk(2,2),qk(2,2),rk(2,2),fk(2) real :: embar0,enmbar,rmass(2),fo2ox(nlev,2), | wk1,wk2,wk3,st,fs11st,fs21st,fs12st,fs22st, | z(nlev,2),fo2oxout(nlev,2) ! ! write(6,"('Enter comp: ut=',f5.2)") ut t0(:) = 0. do k=1,nlev ! was do 1 in old sub comp embar(k) = 1./(ps_o2(k)*rmo2inv + ps_ox(k)*rmoinv + | (1.-ps_o2(k)-ps_ox(k))*rmn2inv) enddo ! k=1,nlev ! write(6,"('comp: embar=',/,(6e12.4))") embar ps1b = 0.24 fb(1) = 2.*ps1b ! was C(48)=2.*C(46), where C(46)=ps1b fb(2) = ps2b ! was C(49) b(1,1) = -1. ! B(2,2) was equivalenced to C(40-43) b(1,2) = 0. b(2,1) = 0. b(2,2) = -1. ps0(1) = b(1,1)*ps_o2(1)+b(1,2)*ps_ox(1)+fb(1) ps0(2) = b(2,1)*ps_o2(1)+b(2,2)*ps_ox(1)+fb(2) embar0 = 1./(ps0(1)*rmo2inv + ps0(2)*rmoinv + | (1.-ps0(1)-ps0(2))*rmn2inv) wk3 = (embar(1)-embar0)/(dzp*(embar0+embar(1))) rmass(1) = rmo2 ; rmass(2) = rmo ! write(6,"('ps0=',2e12.4,' embar0=',e12.4)") ps0,embar0 ! ! Calculate ep(2,2) and ak(2,2,2) at level 1/2 (set z(0) to zero) km = 1 kp = 2 do l=1,2 ! was DO 2 ep(l,kp) = 1.-(2./(embar0+embar(1)))*(rmass(l)+ | (embar(1)-embar0)/dzp) z(1,l) = 0. enddo ! l=1,2 ! write(6,"('comp: kp=',i2,' ep(1:2,kp)=',2e12.4)") kp,ep(1:2,kp) ! ! delta(2,2): delta(1,1) = 1. delta(2,1) = 0. delta(1,2) = 0. delta(2,2) = 1. ! ! phi(2,3): phi(1,1)=0. phi(2,1)=0.673 phi(1,2)=1.35 phi(2,2)=0. phi(1,3)=1.11 phi(2,3)=0.769 ! fo2ox(:,1) = ps_o2(:) fo2ox(:,2) = ps_ox(:) do l=1,2 ! was DO 3 ll = 3-l do m=1,2 ak(l,m,kp) = -delta(l,m)*(phi(ll,3)+(phi(ll,l)-phi(ll,3))* | .5*(ps0(l)+fo2ox(1,l)))-(1.-delta(l,m))*(phi(l,m)-phi(l,3))* | .5*(ps0(l)+fo2ox(1,l)) enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: kp=',i2,' ak(:,:,kp)=',4e12.4)") kp,ak(:,:,kp) ! ! wk1 = mbar/m3*(t00/(t0+t))*0.25/(tau*det) at level 1/2 wk1 = 0.5*(embar0+embar(1))*rmn2inv*(t00/(t0(1)+ftn(nlev)))**0.25/ | (tau*(ak(1,1,kp)*ak(2,2,kp)-ak(1,2,kp)*ak(2,1,kp))) ! ! Complete calculation of ak(1/2) do l=1,2 ! was DO 5 do m=1,2 ak(l,m,kp) = ak(l,m,kp)*wk1 gama(1,l,m) = 0. enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: embar=',/,(6e12.4))") embar ! write(6,"('ps0=',2e12.4,' embar0=',e12.4)") ps0,embar0 ! write(6,"('comp: kp=',i2,' ep(1:2,kp)=',2e12.4)") kp,ep(1:2,kp) ! write(6,"('comp: kp=',i2,' ak(:,:,kp)=',4e12.4)") kp,ak(:,:,kp) ! ! Main k-loop (was DO 6): do k=1,nlev-1 ! ! Cycle km and kp levels: ktmp = km km = kp kp = ktmp ! ! ep(k+1/2): do l=1,2 ! was DO 7 ep(l,kp) = 1.-(2./(embar(k)+embar(k+1)))*(rmass(l)+ | (embar(k+1)-embar(k))/dzp) enddo ! l=1,2 ! write(6,"('comp: k=',i3,' kp=',i3,' ep(1:2,kp)=', ! | 2e12.4)") k,l,ep(1:2,kp) ! ! ak(k+1/2): do l=1,2 ! was DO 8 ll = 3-l do m=1,2 ak(l,m,kp) = -delta(l,m)*(phi(ll,3)+(phi(ll,l)-phi(ll,3))* | .5*(fo2ox(k,l)+fo2ox(k+1,l)))-(1.-delta(l,m))* | (phi(l,m)-phi(l,3))*.5*(fo2ox(k,l)+fo2ox(k+1,l)) enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: k=',i3,' kp=',i3,' ak(:,:,kp)=',4e12.4)") ! | k,kp,ak(:,:,kp) ! ! wk1=mbar/m3*(t00/(t0+t))**0.25/(tau*det(alfa)) wk1=0.5*(embar(k)+embar(k+1))*rmn2inv*(t00/(t0(k+1)+.5* | (ftn(k)+ftn(k+1))))**0.25/(tau*(ak(1,1,kp)*ak(2,2,kp)- | ak(1,2,kp)*ak(2,1,kp))) ! ! Eddy diffusion terms in wk1,wk2: wk2 = wk3 wk3 = (embar(k+1)-embar(k))/(dzp*(embar(k)+embar(k+1))) ! write(6,"('comp: k=',i3,' wk1,2,3=',3e12.4)") k,wk1,wk2,wk3 ! ! Calculate augmented source terms: enmbar = p0*exps(k)*embar(k)/(boltz*ftn(k)) ! write(6,"('comp: k=',i3,' exps,embar=',2e12.4,' boltz=',e12.4, ! | ' tn=',e12.4,' enmbar=',e12.4)") ! | k,exps(k),embar(k),boltz,ftn(k),enmbar ! ! Initial fs11st and fs12st are zero because rk13,14,15 are zero ! (see end of sub face): fs11st=enmbar**2*(-.5/.3*(rk13_comp(k)+rk13_comp(k+1))* | (ps_ox(k)*rmoinv)**2-1./3.*(rk14_comp(k)+rk14_comp(k+1))* | ps_o2(k)*rmo2inv*ps_ox(k)*rmoinv+.5*(rk15_comp(k)+ | rk15_comp(k+1))*(-.5*ps_ox(k)/(rmo*embar(k))+1./3.* | (ps_ox(k)*rmoinv)**2+2./3.*ps_o2(k)*rmo2inv*ps_ox(k)*rmoinv)) fs21st = .5*(fs21(k)+fs21(k+1))+fs11st fs11st = .5*(fs11(k)+fs11(k+1))+fs11st fs12st = | -1./3.*(rk13_comp(k)+rk13_comp(k+1))* ps_o2(k)*rmo2inv* | ps_ox(k)*rmoinv-1./6.*(rk14_comp(k)+rk14_comp(k+1))* | (ps_o2(k)*rmo2inv)**2+.5*(rk15_comp(k)+rk15_comp(k+1))* | ps_o2(k)*rmo2inv*(-.5/embar(k)+1./3.*ps_o2(k)*rmo2inv+ | 2./3.*ps_ox(k)*rmoinv) fs22st=.5*(fs22(k)+fs22(k+1))+enmbar**2*(-(rk12_comp(k)+ | rk12_comp(k+1))*ps_ox(k)/(rmo*embar(k))+fs12st) fs12st=.5*(fs12(k)+fs12(k+1))+enmbar**2*(.5*(rk12_comp(k)+ | rk12_comp(k+1))*ps_ox(k)/(rmo*embar(k))+fs12st) ! ! write(6,"('comp: k=',i3,' fs11st,fs12st,fs21st,fs22st=', ! | 4e12.4)") k,fs11st,fs12st,fs21st,fs22st ! ! Total source (complicated calculation of st is commented out in old code) st = 0. ! ! Set up source matrix in wkm2(2,2): wkm2(1,1) = fs11st-st wkm2(1,2) = 2.*fs12st wkm2(2,1) = .5*fs21st wkm2(2,2) = fs22st-st ! ! Finish calculating ak(k+1/2) and generate pk, qk, rk: stepinv = 1./step ! module data (was C(7)) do l=1,2 ! was DO 10 do m=1,2 ak(l,m,kp) = ak(l,m,kp)*wk1 pk(l,m)=(ak(l,m,km)*(1./dzp+ep(m,km)/2.)-exps(k) | *(expdzinv*difk(k)*dfact*(1./dzp-wk2)+0.25* | (fw(k)+fw(k+1)))*delta(l,m))/dzp rk(l,m)=(ak(l,m,kp)*(1./dzp-ep(m,kp)/2.)-exps(k) | *(expdz*difk(k+1)*dfact*(1./dzp+wk3)-0.25*(fw(k)+ | fw(k+1)))*delta(l,m))/dzp qk(l,m)=-(ak(l,m,km)*(1./dzp-ep(m,km)/2.)+ak(l,m,kp)* | (1./dzp+ep(m,kp)/2.))/dzp+exps(k)*(((expdz*difk(k+1)* | (1./dzp-wk3)+expdzinv*difk(k)*(1./dzp+wk2))*dfact/dzp+ | stepinv)*delta(l,m)-wkm2(l,m)) enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: k=',i3,' pk=',4e12.4)") k,pk ! write(6,"('comp: k=',i3,' rk=',4e12.4)") k,rk ! write(6,"('comp: k=',i3,' qk=',4e12.4)") k,qk ! ! Add explicit source terms and complete calculation of rhs in fk ! (old code comments that shapiro smoother was eliminated) fk(:) = 0. fk(1)=fk(1)- (fs1(k)+fs1(k+1))*rmo/enmbar fk(2)=fk(2)-.5*(fs2(k)+fs2(k+1))*rmo/enmbar do l=1,2 ! was DO 12 fk(l) = exps(k)*(fo2ox(k,l)*stepinv-fk(l)) enddo ! l=1,2 ! write(6,"('comp: k=',i3,' fk=',2e12.4)") k,fk ! ! Lower boundary: if (k==1) then do l=1,2 ! was DO 16 do m=1,2 do kk=1,2 qk(l,m) = qk(l,m)+pk(l,kk)*b(kk,m) enddo ! kk=1,2 enddo ! m=1,2 enddo ! l=1,2 do l=1,2 ! was DO 33 do m=1,2 fk(l) = fk(l)-pk(l,m)*fb(m) pk(l,m) = 0. enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp lower boundary: qk=',4e12.4)") qk ! write(6,"('comp lower boundary: fk=',4e12.4)") fk ! ! Upper boundary: elseif (k==nlev-1) then do l=1,2 ! was DO 17 do m=1,2 qk(l,m)=qk(l,m)+(1.+.5*ep(m,kp)*dzp)/(1.-.5*ep(m,kp)*dzp)* | rk(l,m) rk(l,m) = 0. enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp upper boundary: qk=',4e12.4)") qk endif ! ! qk=alfak=qk-pk*gama(k-1) do l=1,2 ! was DO 18 do m=1,2 do kk=1,2 qk(l,m)=qk(l,m)-pk(l,kk)*gama(k,kk,m) enddo ! kk=1,2 enddo ! m=1,2 enddo ! l=1,2 ! ! wk1 = det(alfa), wkm1 = alfai wk1 = qk(1,1)*qk(2,2)-qk(1,2)*qk(2,1) do l=1,2 ! was DO 20 ll = 3-l do m=1,2 wkm1(l,m)=(delta(l,m)*qk(ll,ll)-(1.-delta(l,m))*qk(l,m))/wk1 enddo ! m=1,2 enddo ! l=1,2 ! ! gama(k+1)=alfai*rk, wkv1=fk-pk*z(k) do l=1,2 ! was DO 21 wkv1(l) = fk(l) do m=1,2 gama(k+1,l,m) = 0. wkv1(l) = wkv1(l)-pk(l,m)*z(k,m) do kk=1,2 gama(k+1,l,m)=gama(k+1,l,m)+wkm1(l,kk)*rk(kk,m) enddo ! kk=1,2 enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: k=',i3,' gama(k+1,:,:)=',4e12.4)") k, ! | gama(k+1,:,:) ! ! z(k+1) = wkm1*wkv1 do l=1,2 ! was DO 231 z(k+1,l) = 0. do m=1,2 z(k+1,l) = z(k+1,l)+wkm1(l,m)*wkv1(m) enddo ! m=1,2 enddo ! l=1,2 ! write(6,"('comp: k=',i3,' z(k+1,:)=',2e12.4)") k,z(k+1,:) enddo ! end main k-loop: k=1,nlev-1 ! ! Update o2,o: ! Set upper boundary zero: fo2oxout(nlev,:) = 0. ! ! Downward sweep (was DO 25): do kk=1,nlev-1 ! 1->96 (nlev==97) k = nlev-kk ! 96->1 do l=1,2 fo2oxout(k,l) = z(k+1,l) do m=1,2 fo2oxout(k,l) = fo2oxout(k,l)-gama(k+1,l,m)*fo2oxout(k+1,m) enddo ! m=1,2 enddo ! l=1,2 enddo ! kk=1,nlev-1 ! ! Insert value of ps(nlev) using boundary condition: do l=1,2 ! was DO 27 fo2oxout(nlev,l) = (1.+.5*ep(l,kp)*dzp)/(1.-.5*ep(l,kp)*dzp)* | fo2oxout(nlev-1,l) enddo ! ! Insure non-negative psi: do l=1,2 do k=1,nlev if (fo2oxout(k,l) < 0.) fo2oxout(k,l) = 0. enddo enddo ! ! Transfer to output arrays: o2out(:) = fo2oxout(:,1) oxout(:) = fo2oxout(:,2) ! write(6,"(/,'End comp: o2out=',/,(6e12.4))") o2out ! write(6,"(/,'End comp: oxout=',/,(6e12.4))") oxout end subroutine comp !----------------------------------------------------------------------- subroutine photo(xo2lb,xolb,xn2lb) use cons,only: xnoxlb,xch4lb,xn4slb,xh2lb,xco2lb,xcolb,xh2olb, | xhelb,xarglb implicit none ! ! Args: real :: xo2lb,xolb,xn2lb ! cm-3 lbc saved before comp update ! ! Local: real :: fm ! fm = xo2lb+xolb+xn2lb phno = xnoxlb*fm phn4s = xn4slb*fm phch4 = xch4lb*fm phh2 = xh2lb *fm phco2 = xco2lb*fm phco = xcolb *fm phh2o = xh2olb*fm phhe = xhelb *fm pharg = xarglb*fm ! write(6,"('photo: fm=',e12.4)") fm ! write(6,"('photo: xnoxlb=',e12.4,' phno =',e12.4)") xnoxlb,phno ! write(6,"('photo: xn4slb=',e12.4,' phn4s=',e12.4)") xn4slb,phn4s ! write(6,"('photo: xch4lb=',e12.4,' phch4=',e12.4)") xch4lb,phch4 ! write(6,"('photo: xh2lb =',e12.4,' phh2 =',e12.4)") xh2lb ,phh2 ! write(6,"('photo: xco2lb=',e12.4,' phco2=',e12.4)") xco2lb,phco2 ! write(6,"('photo: xcolb =',e12.4,' phco= ',e12.4)") xcolb ,phco ! write(6,"('photo: xh2olb=',e12.4,' phh2o=',e12.4)") xh2olb,phh2o ! write(6,"('photo: xhelb =',e12.4,' phhe= ',e12.4)") xhelb ,phhe ! write(6,"('photo: xarglb=',e12.4,' pharg=',e12.4)") xarglb,pharg end subroutine photo !----------------------------------------------------------------------- subroutine minor(finout,fname,rmass,src1,src2,lbc,ubc,phix, | alfax,ftn,fw,ibnd,ibndb,pxb,iprint) use cons,only: rmo2,rmo,rmn2,boltz,p0,grav,avo implicit none ! ! Args: real,intent(inout) :: finout(nlev) real,intent(out) :: pxb real,intent(in) :: rmass,phix(3),alfax real,intent(inout) :: lbc(3),ubc real,dimension(nlev),intent(in) :: src1,src2,ftn,fw integer,intent(in) :: ibnd,ibndb,iprint character(len=*),intent(in) :: fname ! ! Local: integer :: k real,dimension(nlev) :: | hadvec, ! horizontal advection (was S13) | src1i,src2i, ! k+1/2 from src2 (was S14,S15) | mbari, ! mbar at k+1/2 (was S12) | mbar, ! mbar at k (was S11) | dmdz, ! dm/dz(k) (was S10) | dmdzmbar, ! (dm/dz(k))/mbar (was S9) | pso, ! o (was S9) | pso2, ! o2 (was S8) | do, ! do/dz (was S7) | do2, ! do2/dz (was S6) | tntot, ! t0+tn (was S5) (t0 should be zero, so tntot==tn) | alpha11,alpha12,alpha21,alpha22, ! S1,S2,S3,S4 | thermdif, ! thermal diffusion term (was S12) | ex, ! was S12 | ax, ! was S11 | pcoef,qcoef,rcoef,fcoef ! coefficients for trsolv (S1,S2,S3,S4) real :: | pso2lb,psolb, ! o2,o lbc at k-1/2 (was T4,T5) | mbarlb, ! mean molecular weight at k-1/2 (bottom) (was T6) | alfa12,alfa21,alfax1,alfax2 real :: dum,small real,parameter :: dzpinv = 1./dzp character(len=80) :: char80 ! if (iprint > 0) then write(6,"('Enter Minor for field ',a,': rmass=', | f8.2,' phix=',3e12.4,' alfax=',e12.4,' ibnd,ibndb=',2i3)") | trim(fname),rmass,phix,alfax,ibnd,ibndb write(6,"('Minor for field ',a,' src1=',/,(6e12.4))") | trim(fname),src1 write(6,"('Minor for field ',a,' src2=',/,(6e12.4))") | trim(fname),src2 endif ! ! 8/24/05 btf: changed k=1,nlev to k=1,nlev-1. hadvec(:) = 0. do k=1,nlev-1 src1i(k) = .5*(src1(k)+src1(k+1)) ! S14 src2i(k) = .5*(src2(k)+src2(k+1)) ! S15 enddo ! k=1,nlev if (iprint > 0) then write(6,"('Minor for field ',a,': src1i=',/,(6e12.4))") | trim(fname),src1i write(6,"('Minor for field ',a,': src2i=',/,(6e12.4))") | trim(fname),src2i endif pso2lb = b(1,1)*ps_o2(1)+b(1,2)*ps_ox(1)+fb(1) psolb = b(2,1)*ps_o2(1)+b(2,2)*ps_ox(1)+fb(2) mbarlb = | 1./(pso2lb*rmo2inv+psolb*rmoinv+(1.-pso2lb-psolb)*rmn2inv) if (iprint > 0) then write(6,"('Minor: fb=',2e12.4)") fb write(6,"('Minor: b=',4e12.4)") b write(6,"('Minor: ps_o2(1)=',e12.4, | ' ps_ox(1)=',e12.4)") ps_o2(1),ps_ox(1) write(6,"('Minor for field ',a,': pso2lb,psolb,mbarlb=', | 3e12.4)") trim(fname),pso2lb,psolb,mbarlb endif do k=1,nlev mbari(k) = 1./(ps_o2(k)*rmo2inv+ps_ox(k)*rmoinv+ | (1.-ps_o2(k)-ps_ox(k))*rmn2inv) enddo ! k=1,nlev if (iprint > 0) | write(6,"('Minor for field ',a,': mbari=',/,(6e12.4))") | trim(fname),mbari mbar(1) = .5*(mbarlb+mbari(1)) ! S11 dmdz(1) = (mbari(1)-mbarlb)*dzpinv ! S10 pso (1) = .5*(psolb +ps_ox(1)) ! S9 pso2(1) = .5*(pso2lb+ps_o2(1)) ! S8 do (1) = (ps_ox(1)-psolb )*dzpinv ! S7 do2 (1) = (ps_o2(1)-pso2lb)*dzpinv ! S6 do k=2,nlev ! was DO 5 mbar(k) = .5*(mbari(k)+mbari(k-1)) ! S11 dmdz(k) = (mbari(k)-mbari(k-1))*dzpinv ! S10 pso (k) = .5*(ps_ox(k)+ps_ox(k-1)) ! S9 pso2(k) = .5*(ps_o2(k)+ps_o2(k-1)) ! S8 do (k) = (ps_ox(k)-ps_ox(k-1))*dzpinv ! S7 do2 (k) = (ps_o2(k)-ps_o2(k-1))*dzpinv ! S6 enddo ! k=2,nlev if (iprint > 0) then write(char80,"('Minor: field ',a)") trim(fname) call fprint8(char80,nlev, | mbar , dmdz , pso , pso2 , do , do2 , ftn,dum, | 'mbar','dmdz','pso','pso2','do','do2','tn',' ') endif tntot(1) = t0(1)+ftn(nlev) ! lower boundary in top slot? tntot(nlev) = t0(nlev)+ftn(nlev-1) do k=2,nlev-1 ! was DO 7 tntot(k) = t0(k)+.5*(ftn(k)+ftn(k-1)) enddo ! k=2,nlev if (iprint > 0) | write(6,"('Minor for field ',a,': tntot=',/,(6e12.4))") | trim(fname),tntot alfa12 = phi(1,2)-phi(1,3) alfa21 = phi(2,1)-phi(2,3) alfax1 = phix(1)-phix(3) alfax2 = phix(2)-phix(3) if (iprint > 0) | write(6,"('Minor for field ',a,': alfa12,21,x1,x2=',4e12.4)") | trim(fname),alfa12,alfa21,alfax1,alfax2 ! 8/24/05 btf: changed k=1,nlev to k=1,nlev-1. do k=1,nlev-1 ! was DO 8 src2i(k) = src2i(k)*rmass*boltz*(.5*(t0(k)+t0(k+1))+ftn(k))/ | (p0*exps(k)*mbari(k)) enddo ! k=1,nlev if (iprint > 0) | write(6,"('Minor for field ',a,': src2i=',/,(6e12.4))") | trim(fname),src2i if (ibndb==1) then ! ibndb is normally 0 lbc(1) = 1. ! T1 lbc(2) = dmdz(1)/mbar(1) ! T2 lbc(3) = lbc(3)*grav*rmass/(difk(1)*p0*exps(1)*expdzinv*avo) ! T3 endif if (iprint > 0) | write(6,"('Minor for field ',a,': lbc=',3e12.4)") | trim(fname),lbc ! do k=1,nlev ! was DO 9 alpha11(k) = -(phi(1,3)+alfa12*pso(k)) ! S1 alpha12(k) = alfa12*pso2(k) ! S2 alpha21(k) = alfa21*pso(k) ! S3 alpha22(k) = -(phi(2,3)+alfa21*pso2(k)) ! S4 ! ! ex==S12, dmdz==S10, ax==S11 ex(k) = ((alfax1*alpha22(k)-alfax2*alpha21(k))*(do2(k)- | (1.-(rmo2+dmdz(k))/mbar(k))*pso2(k))+(alfax2*alpha11(k)- | alfax1*alpha12(k))*(do(k)-(1.-(rmo+dmdz(k))/mbar(k))*pso(k)))/ | (alpha11(k)*alpha22(k)-alpha12(k)*alpha21(k))+ | 1.-(rmass+dmdz(k))/mbar(k) dmdzmbar(k) = dmdz(k)/mbar(k) ! S10 ax(k) = -mbar(k)/(tau*rmn2)*(t00/tntot(k))**0.25/(phix(3)+ | alfax1*pso2(k)+alfax2*pso(k)) enddo ! k=1,nlev if (iprint > 0) then write(char80,"('Minor alphas: field ',a)") trim(fname) call fprint8(trim(char80),nlev, | alpha11 , alpha12 , alpha21 , alpha22 , ex , dmdzmbar, ax ,dum, | 'alpha11','alpha12','alpha21','alpha22','ex','dmdzmbar','ax', | ' ') endif ! ! thermdif = ex-alfax*d/ds(ln(t(tot)) (thermal diffusion term) (S12) do k=2,nlev-1 thermdif(k) = ex(k)-alfax*(tntot(k+1)-tntot(k-1))/ | (2.*dzp*tntot(k)) enddo ! k=2,nlev-1 thermdif(1) = ex(1)-alfax*(tntot(2)-tntot(1))/ | (dzp*tntot(1)) thermdif(nlev) = ex(nlev)-alfax*(tntot(nlev)-tntot(nlev-1))/ | (dzp*tntot(nlev)) if (iprint > 0) | write(6,"('Minor for field ',a,': thermdif=',/,(6e12.4))") | trim(fname),thermdif ! ! (Here, old code sets s9=s10 (s9 is now dmdzmbar), and s10=finout) ! ! pcoef==S1, qcoef==S2, rcoef==S3, fcoef==S4 ! C(86)==expdz, C(87)==expdzinv, S9==dmdzmbar do k=1,nlev-1 ! was DO 13 pcoef(k) = ax(k)/dzp*(1./dzp+.5*thermdif(k))-exps(k)* | (expdzinv*difk(k)*dfact*(1./dzp-.5*dmdzmbar(k))+ | 0.25*(fw(k)+fw(k+1)))/dzp rcoef(k) = ax(k+1)/dzp*(1./dzp-.5*thermdif(k+1))-exps(k)* | (expdz*difk(k+1)*dfact*(1./dzp+.5*dmdzmbar(k+1))- | 0.25*(fw(k)+fw(k+1)))/dzp qcoef(k) = -(ax(k)/dzp*(1./dzp-.5*thermdif(k))+ax(k+1)/dzp* | (1./dzp+.5*thermdif(k+1)))+exps(k)*((expdzinv*difk(k)* | (1./dzp+.5*dmdzmbar(k))+expdz*difk(k+1)*(1./dzp-.5* | dmdzmbar(k+1)))*dfact/dzp-src1i(k)+stepinv) fcoef(k) = exps(k)*(finout(k)*stepinv-hadvec(k)+src2i(k)) enddo ! k=1,nlev ! if (iprint > 0) write(6,"('Minor field ',a,' pcoef=',/,(6e12.4))") ! | trim(fname),pcoef ! if (iprint > 0) write(6,"('Minor field ',a,' qcoef=',/,(6e12.4))") ! | trim(fname),qcoef ! if (iprint > 0) write(6,"('Minor field ',a,' rcoef=',/,(6e12.4))") ! | trim(fname),rcoef ! if (iprint > 0) write(6,"('Minor field ',a,' fcoef=',/,(6e12.4))") ! | trim(fname),fcoef ! ! Boundaries: modify ex(nlev) if ibnd==1: if (ibnd==1) then ! ibnd is normally zero thermdif(nlev) = thermdif(nlev)-(exps(nlev-1)*expdz*fw(nlev)/ | ax(nlev)) endif ! ! Lower boundary qcoef(1) = qcoef(1)+pcoef(1)*(lbc(1)+.5*lbc(2)*dzp)/ | (lbc(1)-.5*lbc(2)*dzp) if (iprint > 0) write(6,"('Minor field ',a,' pcoef(1)=',e12.4, | ' lbc(1:2)=',2e12.4,' qcoef(1)=',e12.4)") trim(fname), | pcoef(1),lbc(1:2),qcoef(1) fcoef(1) = fcoef(1)-pcoef(1)*lbc(3)*dzp/(lbc(1)-.5*lbc(2)*dzp) pcoef(1) = 0. ! ! Upper boundary (upward flux given) qcoef(nlev-1) = qcoef(nlev-1)+rcoef(nlev-1)*(1.+.5*thermdif(nlev)* | dzp)/(1.-.5*thermdif(nlev)*dzp) fcoef(nlev-1) = fcoef(nlev-1)-rcoef(nlev-1)*grav*rmass/ | (p0*ax(nlev)*avo)*ubc*dzp/(1.-.5*thermdif(nlev)*dzp) rcoef(nlev-1) = 0. if (iprint > 0) then call fprint8(trim(char80),nlev-1, | thermdif , pcoef , qcoef , rcoef , fcoef ,dum,dum,dum, | 'thermdif','pcoef','qcoef','rcoef','fcoef',' ',' ',' ') endif ! ! Tri-diagonal solver (note nlon==i1==i2==1): ! subroutine trsolv(a,b,c,f,x,nlev,k1,k2,nlon,i1,i2,iprint) ! call trsolv(pcoef,qcoef,rcoef,fcoef,finout,nlev,1,nlev-1,1,1,1,0) ! ! Upper boundary: finout(nlev) = (grav*rmass/(p0*ax(nlev)*avo)*ubc*dzp+ | (1.+.5*thermdif(nlev)*dzp)*finout(nlev-1))/ | (1.-.5*thermdif(nlev)*dzp) if (iprint > 0) | write(6,"('Minor for field ',a,': trsolv output (including ', | 'upper boundary) x=',/,(6e12.4))") trim(fname),finout ! ! pbx = psx(-1/2) pxb = (finout(1)*(lbc(1)+.5*lbc(2)*dzp)+lbc(3)*dzp)/ | (lbc(1)-.5*lbc(2)*dzp) ! ! Insure non-negative density: small = 0. do k=1,nlev if (finout(k) < small) then if (iprint > 0) write(6,"('>>> Minor for field ',a,': ', | 'finout(k=',i4,')=',e12.4,' is too small: setting to ', | e12.4)") k,finout(k) finout(k) = small endif enddo if (iprint > 0) | write(6,"('End Minor for field ',a,': pxb=',e12.4, | ' trsolv output x=',/,(6e12.4))") trim(fname),pxb,finout end subroutine minor !----------------------------------------------------------------------- subroutine compn2d(xn2dout) ! ! Calculate n(n2d) assuming photochemical equilibrium. ! use flds_prodloss,only: pn2d,xln2d implicit none ! ! Args: real,dimension(nlev),intent(out) :: xn2dout ! ! Local: integer :: k ! do k=1,nlev xn2dout(k) = pn2d(k)/xln2d(k) enddo ! k=1,nlev ! write(6,"(/,'compn2d: xn2dout=',/,(6e12.4))") xn2dout end subroutine compn2d !----------------------------------------------------------------------- subroutine compn4s(fn4s,tn,w,pn4sb) ! use flds_prodloss,only: pn4s,xln4s use cons,only: rmn4s,rmo,rmo2,rmn2 implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fn4s real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pn4sb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real :: phin4s(3) = (/0.651,0.731,0.741/) ! ! Upper boundary: diffusive equilibrium ! Lower boundary: photochemical equilibrium ! phn4s = pn4s(1)/xln4s(1) ! production/loss ubc = 0. ! was T4 lbc(1) = 0. ! was T1 lbc(2) = 1. ! was T2 lbc(3) = -phn4s*rmn4s/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) ! write(6,"('compn4s: phn4s=',e12.4,' xn_o2,xno_lb,xn_n2(1)=', ! | 3e12.4,' lbc(3)=',e12.4)") phn4s,xn_o2(1),xno_lb,xn_n2(1), ! | lbc(3) source1(:) = -xln4s(:) source2(:) = pn4s(:) ibnd = 0 ibndb = 0 alfa = 0. call minor(fn4s,'N4S',rmn4s,source1,source2,lbc,ubc, | phin4s,alfa,tn,w,ibnd,ibndb,pn4sb,0) ! write(6,"(/,'compn4s after call minor: pn4sb=',e12.4,' fn4s=',/, ! | (6e12.4))") pn4sb,fn4s end subroutine compn4s !----------------------------------------------------------------------- subroutine compnoz(fnoz,tn,w,pnozb) use cons,only: rmno,rmo,rmo2,rmn2 use flds_prodloss,only: xlnox,pnox implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fnoz real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pnozb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real :: phinoz(3) = (/0.814,0.866,0.926/) real,parameter :: fluxno = 0. ! ! ibndno=0 -> flux at lower boundary, ! ibndno=1 -> photochemical equilibrium integer :: ibndno = 1 ! ubc = 0. do k=1,nlev source1(k) = -xlnox(k) source2(k) = pnox(k) enddo ! k=1,nlev ! if (ibndno==0) then lbc(1:2) = 0. lbc(3) = fluxno ibndb = 1 else ubc = 0. lbc(1) = 0. lbc(2) = 1. lbc(3) = -phno*rmno/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) ibndb = 0 endif ibnd =0 alfa = 0. call minor(fnoz,'NOZ',rmno,source1,source2,lbc,ubc, | phinoz,alfa,tn,w,ibnd,ibndb,pnozb,0) ! write(6,"(/,'compnoz after call minor: pnozb=',e12.4,' fnoz=',/, ! | (6e12.4))") pnozb,fnoz end subroutine compnoz !----------------------------------------------------------------------- subroutine compar(far,tn,w,parb) use cons,only: rmar,rmo,rmo2,rmn2 use flds_prodloss,only: xlnox,pnox implicit none ! ! Args: real,dimension(nlev),intent(inout) :: far real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: parb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real :: phiar(3) = (/1.042,1.509,1.176/) ! ubc = 0. lbc(1) = 0. lbc(2) = 1. lbc(3) = -pharg*rmar/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = 0.17 ! ! Sources are zero: source1 = 0. source2 = 0. ibnd = 0 ibndb = 0 call minor(far,'AR',rmar,source1,source2,lbc,ubc, | phiar,alfa,tn,w,ibnd,ibndb,parb,0) ! write(6,"(/,'compar after call minor: parb=',e12.4,' far=',/, ! | (6e12.4))") parb,far end subroutine compar !----------------------------------------------------------------------- subroutine comphe(fhe,tn,w,pheb) use cons,only: rmhe,rmo,rmo2,rmn2 implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fhe real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pheb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phihe(3) = (/.270,.404,.322/) ! ubc = 0. lbc(1) = 0. lbc(2) = 1. lbc(3) = -phhe*rmhe/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = -0.38 source1(:) = 0. source2(:) = 0. ibnd = 1 ibndb = 0 call minor(fhe,'HE',rmhe,source1,source2,lbc,ubc, | phihe,alfa,tn,w,ibnd,ibndb,pheb,0) ! write(6,"(/,'comphe after call minor: pheb=',e12.4,' fhe=',/, ! | (6e12.4))") pheb,fhe end subroutine comphe !----------------------------------------------------------------------- subroutine compch4(fch4,tn,w,pch4b) use cons,only: rmch4,rmo,rmo2,rmn2 use flds_prodloss,only: xlch4 implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fch4 real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pch4b ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phich4(3) = (/0.921,0.846,1.077/) ! ! Upper boundary: diffusive equilibrium ubc = 0. ! ! Sources: source1(:) = -xlch4(:) source2(:) = 0. ! ! Density specified at lower boundary: lbc(1) = 0. lbc(2) = 1. lbc(3) = -phch4*rmch4/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) ! alfa = 0. ibnd = 0 ibndb = 0 call minor(fch4,'CH4',rmch4,source1,source2,lbc,ubc, | phich4,alfa,tn,w,ibnd,ibndb,pch4b,0) ! write(6,"(/,'compch4 after call minor: pch4b=',e12.4,' fch4=',/, ! | (6e12.4))") pch4b,fch4 end subroutine compch4 !----------------------------------------------------------------------- subroutine comph2(fh2,tn,w,ph2b) use cons,only: rmh2,rmo,rmo2,rmn2 use flds_prodloss,only: xlh2,ph2 implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fh2 real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: ph2b ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phih2(3) = (/0.226,0.321,0.282/) ! ubc = 0. source1(:) = -xlh2(:) ! xlh2 is loss rate (module flds_prodloss) source2(:) = ph2(:) ! ph2 is production (module flds_prodloss) lbc(1) = 0. lbc(2) = 1. lbc(3) = -phh2*rmh2/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = 0. ibnd = 0 ibndb = 0 call minor(fh2,'H2',rmh2,source1,source2,lbc,ubc, | phih2,alfa,tn,w,ibnd,ibndb,ph2b,0) ! write(6,"(/,'comph2 after call minor: ph2b=',e12.4,' fh2=',/, ! | (6e12.4))") ph2b,fh2 end subroutine comph2 !----------------------------------------------------------------------- subroutine compco2(fco2,tn,w,pco2b) use cons,only: rmco2,rmo,rmo2,rmn2 use flds_prodloss,only: xlco2,pco2 implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fco2 real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pco2b ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phico2(3) = (/1.199,3.91,1.201/) ! ubc = 0. source1(:) = -xlco2(:) ! xlco2 is loss rate (module flds_prodloss) source2(:) = pco2(:) ! pco2 is production (module flds_prodloss) lbc(1) = 0. lbc(2) = 1. lbc(3) = -phco2*rmco2/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = 0. ibnd = 0 ibndb = 0 call minor(fco2,'CO2',rmco2,source1,source2,lbc,ubc, | phico2,alfa,tn,w,ibnd,ibndb,pco2b,0) ! write(6,"(/,'compco2 after call minor: pco2b=',e12.4,' fco2=',/, ! | (6e12.4))") pco2b,fco2 end subroutine compco2 !----------------------------------------------------------------------- subroutine compco(fco,tn,w,pcob) use cons,only: rmco,rmo,rmo2,rmn2 use flds_prodloss,only: xlco,pco implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fco real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pcob ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phico(3) = (/0.833,1.427,0.852/) integer :: ibndco = 1 ! ubc = 0. source1(:) = -xlco(:) ! xlco is loss rate (module flds_prodloss) source2(:) = pco(:) ! pco is production (module flds_prodloss) if (ibndco == 0) then lbc(1) = 0. lbc(2) = 0. lbc(3) = 0. ! was fluxco ibndb = 1 else lbc(1) = 0. lbc(2) = 1. lbc(3) = -phco*rmco/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) ibndb = 0 endif alfa = 0. ibnd = 0 call minor(fco,'CO',rmco,source1,source2,lbc,ubc, | phico,alfa,tn,w,ibnd,ibndb,pcob,0) ! write(6,"(/,'compco after call minor: pcob=',e12.4,' fco=',/, ! | (6e12.4))") pcob,fco end subroutine compco !----------------------------------------------------------------------- subroutine hoxubc(xnh_top,tn_top,gz_top,zpht_top) use cons,only: boltz implicit none ! ! Args: real,intent(in) :: xnh_top,tn_top,gz_top,zpht_top ! ! Local: real :: alfa,hc1,xnc1,pi,b,hc,alamc,u,fac ! hc1=boltz*1000./(1.66e-24*gz_top) xnc1=2.75e+4 pi=3.1415926 b=0.72 hc=boltz*tn_top/(1.66e-24*gz_top) alamc=(6371.e+5+zpht_top*1.e+5)/hc u=sqrt(2.*boltz*tn_top/1.66e-24) fac=b*u*(1.+alamc)*exp(-alamc)/(2.*sqrt(pi)) ! ! phij,phip,phie, and fluxhox are in module data: ! ! write(6,"('hoxubc: xnh_top=',e12.4,' tn_top=',e12.4, ! | ' gz_top=',e12.4)") xnh_top,tn_top,gz_top ! write(6,"(' fac=',e12.4,' hc=',e12.4,' xnc1=',e12.4, ! | ' hc1=',e12.4)") fac,hc,xnc1,hc1 phij=xnh_top*fac phip=0.36*phij phie=2.8e+8*xnh_top*hc*phij/(5.e+7*xnc1*hc1) fluxhox=phij+phip+phie ! write(6,"('hoxubc: phij,phip,phie=',3e12.4,' fluxhox=',e12.4)") ! | phij,phip,phie,fluxhox end subroutine hoxubc !----------------------------------------------------------------------- subroutine comphox(fhox,ftn,w,phoxb) use cons,only: rmhox,rmo,rmo2,rmn2,boltz use flds_prodloss,only: xlhox,phox use flds_modelz,only: gz,zpht use flds_atmos,only: tn,xnh implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fhox real,dimension(nlev),intent(in) :: ftn,w real,intent(out) :: phoxb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | xhox,alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phihox(3) = (/0.146,0.243,0.162/) ! source1(:) = -xlhox(:) source2(:) = phox(:) lbc(1) = 0. lbc(2) = 1. xhox=phox(1)/xlhox(1) lbc(3) = -xhox/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = -0.38 ibnd = 1 ibndb = 0 call minor(fhox,'HOX',rmhox,source1,source2,lbc,fluxhox, | phihox,alfa,ftn,w,ibnd,ibndb,phoxb,0) ! write(6,"(/,'comphox after call minor: phoxb=',e12.4,' fhox=',/, ! | (6e12.4))") phoxb,fhox end subroutine comphox !----------------------------------------------------------------------- subroutine comph2o(fh2o,tn,w,ph2ob) use cons,only: rmh2o,rmo,rmo2,rmn2 use flds_prodloss,only: xlh2o,ph2o implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fh2o real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: ph2ob ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phih2o(3) = (/0.817,0.922,0.920/) ! ubc = 0. source1(:) = -xlh2o(:) source2(:) = ph2o(:) lbc(1) = 0. lbc(2) = 1. lbc(3) = -phh2o*rmh2o/(xn_o2(1)*rmo2+xno_lb*rmo+xn_n2(1)*rmn2) alfa = 0. ibnd = 0 ibndb = 0 ! write(6,"(/,'comph2o before minor: lbc(3)=',e12.4,' fh2o=',/, ! | (6e12.4))") lbc(3),fh2o call minor(fh2o,'H2O',rmh2o,source1,source2,lbc,ubc, | phih2o,alfa,tn,w,ibnd,ibndb,ph2ob,0) ! write(6,"(/,'comph2o after call minor: ph2ob=',e12.4,' fh2o=',/, ! | (6e12.4))") ph2ob,fh2o end subroutine comph2o !----------------------------------------------------------------------- subroutine compnat(fnat,tn,w,pnatb) use cons,only: rmnat,rmo,rmo2,rmn2 use flds_prodloss,only: xlna,pna implicit none ! ! Args: real,dimension(nlev),intent(inout) :: fnat real,dimension(nlev),intent(in) :: tn,w real,intent(out) :: pnatb ! ! Local: integer :: k,ibnd,ibndb real :: | lbc(3), ! was t1,t2,t3 lower boundary | ubc, ! was t4 lower boundary | alfa real,dimension(nlev) :: | source1, ! was s1 | source2 ! was s2 real,parameter :: phinat(3) = (/1.042,1.509,1.176/) ! ubc = 0. source1(:) = 0. source2(:) = pna(:) lbc(1) = 0. lbc(2) = 1. lbc(3) = -1.65e+3 alfa = 0. ibnd = 0 ibndb = 1 ! flux boundary condition call minor(fnat,'NAT',rmnat,source1,source2,lbc,ubc, | phinat,alfa,tn,w,ibnd,ibndb,pnatb,0) ! write(6,"(/,'compnat after call minor: pnatb=',e12.4,' fnat=',/, ! | (6e12.4))") pnatb,fnat end subroutine compnat !----------------------------------------------------------------------- subroutine sodium use flds_ratecoef,only: rkn1,rkn2,rkn3a,rkn3b,rkn4,rkn5,rkn6,rkn7, | rkn8,rkn9,rkn11,rkn13,rkn15,rkn16,rkn17,rkn19,rkn20,rkn21,rkn22, | rkn24a,rkn24b,rkn25,rkn27a,rkn27b,rkn28,rkn29,rkn30,xjn10,xjn12, | xjn14,xjn26 use flds_atmos,only: xno,xno3,xnn2,xno2,xnh,xnco2,xnh2o,xnh2, | xnh2o use flds_ionatm,only: xne,xinop,xio2p use flds_sodium,only: xnas,xnao,xnao3,xnaco3,xnaoh,xnao2,xnahco3, | xnan2p,xnaop,xnasp,xnaco2p,xnah2op,xnat ! ! Local: integer :: k real,dimension(nlev) :: an,bn,cn,dn,en,fn,gn,hn,denom,denon,nn real,dimension(nlev) :: r1,r2,r3,r4,r5,r6,r7,r8,r9,rt,aa,bb,cc,dd real :: dum ! ! Sodium partitioning: do k=1,nlev an(k) = rkn5(k)*xno2(k)*xnn2(k)/(rkn19(k)*xno(k)+xjn26(k)) bn(k)=rkn6(k)*xnco2(k)*xnn2(k)/(rkn20(k)*xno(k)+rkn21(k)*xnh(k)) cn(k)=(rkn4(k)*xnh2o(k)+rkn24a(k)*xnh2(k)+rkn21(k)*xnh(k)* | bn(k))/(rkn11(k)*xnh(k)+xjn12(k)+rkn13(k)* | xnn2(k)*xnco2(k)) denom(k)=rkn8(k)*xno(k)+rkn9(k)*xnh(k)+xjn10(k) dn(k)=(rkn3a(k)*xno3(k)+rkn20(k)*xno(k)*bn(k))/denom(k) en(k)=rkn7(k)*xno2(k)*xnn2(k)/denom(k) fn(k)=rkn13(k)*xnco2(k)*xnn2(k)*cn(k)/(rkn22(k)*xnh(k)) denon(k)=rkn2(k)*xno(k)+(rkn3a(k)+rkn3b(k))*xno3(k)+ | rkn4(k)*xnh2o(k)+ | rkn5(k)*xno2(k)*xnn2(k)+rkn6(k)*xnco2(k)*xnn2(k)+ | (rkn24a(k)+rkn24b(k))*xnh2(k)+rkn25(k)*xnh(k) gn(k)=(rkn1(k)*xno3(k)+rkn8(k)*xno(k)*en(k))/denon(k) hn(k)=(xjn26(k)*an(k)+rkn8(k)*xno(k)*dn(k))/denon(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(k)+rkn15(k)*xio2p(k)+rkn16(k)*xinop(k) bb(k)=rkn29(k)*xno(k)*rkn28(k)*xno(k)/(rkn29(k)*xno(k)+ | rkn30(k)*xne(k)) cc(k)=rkn17(k)*xnn2(k)*xnn2(k) dd(k)=rkn27a(k)*xnco2(k)+rkn27b(k)*xnh2o(k)+rkn28(k)*xno(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)*xno(k)/(rkn29(k)*xno(k)+rkn30(k)*xne(k))+ | (rkn27a(k)*xnco2(k)+rkn27b(k)*xnh2o(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)) xnas(k)=rt(k)*xnat(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) xnan2p(k)=r9(k)*r7(k)*xnas(k) xnaop(k)=rkn28(k)*xno(k)/(rkn29(k)*xno(k)+rkn30(k)*xne(k)) | *xnan2p(k) xnasp(k)=(aa(k)*xnas(k)+rkn29(k)*xno(k)*xnaop(k))/cc(k) xnaco2p(k)=rkn27a(k)*xnco2(k)/(rkn30(k)*xne(k))*xnan2p(k) xnah2op(k)=rkn27b(k)*xnh2o(k)/(rkn30(k)*xne(k))*xnan2p(k) enddo ! k=1,nlev ! call fprint8('sodium partitioning',nlev, ! | an , bn , cn , denom , dn , en , fn , dum, ! | 'an','bn','cn','denom','dn','en','fn',' ') ! call fprint8('sodium partitioning',nlev, ! | denon , gn , hn , nn , r1 , r2 , r3 , dum, ! | 'denon','gn','hn','nn','r1','r2','r3',' ') ! call fprint8('sodium partitioning',nlev, ! | r4 , r5 , r6 , aa , bb , cc , dd , dum, ! | 'r4','r5','r6','aa','bb','cc','dd',' ') ! call fprint8('sodium partitioning',nlev, ! | r7 , r8 , r9 , rt , dum , dum , dum , dum, ! | 'r7','r8','r9','rt',' ',' ',' ',' ') ! call fprint8('sodium fields',nlev, ! | xnas , xnao , xnao3 , xnaco3 , xnaoh , xnao2 , xnahco3 , dum, ! | 'xnas','xnao','xnao3','xnaco3','xnaoh','xnao2','xnahco3',' ') ! call fprint8('sodium fields',nlev, ! | xnan2p , xnaop , xnasp , xnaco2p , xnah2op , dum , dum , dum, ! | 'xnan2p','xnaop','xnasp','xnaco2p','xnah2op',' ',' ',' ') end subroutine sodium !----------------------------------------------------------------------- end module composition