! module o2o_module use params_module,only: nlonp4 implicit none real :: xolb(nlonp4) contains !----------------------------------------------------------------------- subroutine comp_o2o(tn,o2,o1,barm,op,no,n4s,n2d,o2p,ne,o1d, | o3,n2p,nplus,nop,xiop2p,xiop2d,lev0,lev1,lon0,lon1,lat) ! ! Calculates fs array which gives sources and sinks for comp (o2,o). ! Chemrates: real,allocatable :: fs(:,:,:,:,:) ! (i,k,2,0:2,j) ! use cons_module,only: rmassinv_o2,rmassinv_o1,rmassinv_n2,p0, | expz,boltz,rmass_no,rmass_n4s,rmass_n2d,rmass_o2,rmass_o1, | rmassinv_o3 use qrj_module,only: rj,qo2p,qop,pdo3d use chemrates_module,only: rk4,rk5,rk6,rk7,rk8,rk9,rk10,beta2, | beta6,ra1,ra2,beta3,beta8,rk1,beta1,rkm12,rk3,fs,rk16,rk24, | rk23,rk27,rkm7a,rkm7b,rkm24,rkm21,beta9 use addfld_module,only: addfld implicit none ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat ! ! All input arrays are at current time step (itp): real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | tn, ! neutral temperature (deg K) | o2, ! O2 (mmr) | o1, ! O1 (mmr) | barm, ! mean molecular weight | op, ! O+ | no, ! nitric oxide | n4s, ! N(4s) | n2d, ! N(2d) | o2p, ! O2+ | ne, ! electron density | o1d, ! mmr (from comp_o1d) | o3, ! mmr (from comp_o1d) | n2p, ! N(2p) | nplus, ! N+ | nop, ! NO+ | xiop2p,! from oplus | xiop2d ! from oplus ! ! Local: integer :: k,i,nk real :: acof,bcof,ccof real,dimension(lev0:lev1,lon0:lon1) :: | rji, ! RJ at interfaces (qrj module) | qo2pi, ! O2+ ionization at interfaces (qrj module) | qopi, ! O+ ionization at interfaces | xnmbar, ! p0*e(-z)*barm | pox1,pox2, ! OX production terms | lox1,lox2,lox3, ! OX loss terms | po21,po22,po23, ! O2 production terms | lo21,lo22, ! O2 loss terms | n2, ! N2 mmr | o1d_cm3, ! o1d (cm3) | o3_cm3 ! o3 (cm3) ! nk = lev1-lev0+1 ! ! Loop over subdomain at current latitude: do i=lon0,lon1 do k=lev0,lev1-1 ! ! Qrj ionization rates at interfaces: rji (k,i) = 0.5*(rj(k,i,lat)+rj(k+1,i,lat)) qo2pi(k,i) = 0.5*(qo2p(k,i,lat)+qo2p(k+1,i,lat)) qopi (k,i) = 0.5*(qop(k,i,lat)+qop(k+1,i,lat)) ! ! n*mbar (k+1/2): xnmbar(k,i) = expz(k)*p0/(boltz*tn(k,i)*(o2(k,i)*rmassinv_o2+ | o1(k,i)*rmassinv_o1+(1.-o2(k,i)-o1(k,i))*rmassinv_n2)) o3_cm3(k,i) = o3(k,i)*xnmbar(k,i)*rmassinv_o3 o1d_cm3(k,i) = o1d(k,i)*xnmbar(k,i)*rmassinv_o1 ! ! OX production: n2(k,i) = 1.-o2(k,i)-o1(k,i) pox1(k,i) = xnmbar(k,i)**2* | (beta3(k,i,lat)*n4s(k,i)/rmass_n4s*no(k,i)/rmass_no+ | beta6*n2d(k,i)/rmass_n2d*no(k,i)/rmass_no)+ | 0.5*(beta8(k,i,lat)+beta8(k+1,i,lat))* | no(k,i)/rmass_no*xnmbar(k,i)+xnmbar(k,i)* | (rk4*o2p(k,i)*n4s(k,i)/rmass_n4s+ | rk10*op(k,i)*n2d(k,i)/rmass_n2d)+ | (ra1(k,i,lat)*nop(k,i)+2.*ra2(k,i,lat)*o2p(k,i))* | sqrt(ne(k,i)*ne(k+1,i))+ | pdo3d(k,i,lat)*o3_cm3(k,i)+ | (rk16*xiop2p(k,i)*n2(k,i)*rmassinv_n2+rk24*xiop2d(k,i)* | n2(k,i)*rmassinv_n2)*xnmbar(k,i) pox2(k,i) = xnmbar(k,i)*(beta1(k,i,lat)*n4s(k,i)/rmass_n4s+ | beta2*n2d(k,i)/rmass_n2d)+rk1(k,i,lat)*op(k,i)+rk7* | nplus(k,i)+2.*rji(k,i)+ | rk23*xiop2p(k,i)+rk27*xiop2d(k,i) ! ! OX loss: lox1(k,i) = 2.*rkm12(k,i,lat)*xnmbar(k,i)/ | (.5*(barm(k,i)+barm(k+1,i))) lox2(k,i) = rk3(k,i,lat)*n2p(k,i)+rk8*nplus(k,i)+ | rkm21(k,i,lat)*o2(k,i)*rmassinv_o2*xnmbar(k,i)**2/ | (.5*(barm(k,i)+barm(k+1,i)))+rkm24(k,i,lat)*o3_cm3(k,i) lox3(k,i) = qopi(k,i)+rkm7a*o3_cm3(k,i)*o1d_cm3(k,i) ! ! O2 production: po21(k,i) = rkm12(k,i,lat)*xnmbar(k,i)/ | (.5*(barm(k,i)+barm(k+1,i))) po22(k,i) = 2.*rkm24(k,i,lat)*o3_cm3(k,i) po23(k,i) = rk5*no(k,i)/rmass_no*o2p(k,i)*xnmbar(k,i)+ | (2.*rkm7a+rkm7b)*o1d_cm3(k,i)*o3_cm3(k,i)+ | pdo3d(k,i,lat)*o3_cm3(k,i)+beta9(k,i,lat)*no(k,i)/rmass_no* | xnmbar(k,i)*o3_cm3(k,i) ! ! O2 loss: lo21(k,i) = xnmbar(k,i)*(beta1(k,i,lat)*n4s(k,i)/rmass_n4s+ | beta2*n2d(k,i)/rmass_n2d)+rk1(k,i,lat)*op(k,i)+(rk6+rk7)* | nplus(k,i)+rk9*n2p(k,i)+rji(k,i)+ | rkm21(k,i,lat)*o1(k,i)*rmass_o1*xnmbar(k,i)*xnmbar(k,i)/ | (.5*(barm(k,i)+barm(k+1,i)))+rk23*xiop2p(k,i)+rk27* | xiop2d(k,i) ! lo22(k,i) = qo2pi(k,i) ! ! Matrix coefficients for O-O2-N2 solution: fs(i,k,1,1,lat) = -lo21(k,i) fs(i,k,1,2,lat) = xnmbar(k,i)*po21(k,i)*o1(k,i)*rmassinv_o1* | rmass_o2*rmassinv_o1 fs(i,k,2,1,lat) = pox2(k,i)*rmass_o1*rmassinv_o2 fs(i,k,2,2,lat) = -lox2(k,i)-lox1(k,i)*o1(k,i)*rmassinv_o1* | xnmbar(k,i) fs(i,k,1,0,lat) = (po23(k,i)-lo22(k,i))*rmass_o2/xnmbar(k,i) fs(i,k,2,0,lat) = (pox1(k,i)-lox3(k,i))*rmass_o1/xnmbar(k,i) enddo ! k=lev0,lev1-1 ! ! This defines the lower boundary condition for bndry.F but only calculated ! at the k-1 level. acof = lox1(lev0,i) bcof = lox2(lev0,i) ccof = lox3(lev0,i)-pox1(lev0,i)-pox2(lev0,i)*o2(lev0,i)* | rmassinv_o2*xnmbar(lev0,i) xolb(i) = (-bcof-sqrt(bcof**2+4.*acof*ccof))/(2.*acof) ! ! We may have to change this to mass mixing ratio for bndry.F xolb(i) = xolb(i)*rmassinv_o1/xnmbar(lev0,i) enddo ! i=lon0,lon1 ! write(6,"('comp_o2o: xolb=',/,(6e12.4))") xolb(lon0:lon1) ! write(6,"('comp_o2o: fs min,max=',2e12.4)") ! | minval(fs),maxval(fs) ! do k=lev0,lev1-1 ! write(6,"('comp_o2o: lat=',i3,' k=',i3)") lat,k ! write(6,"('fs(:,k,1,1,lat)=',/,(6e12.4))") fs(:,k,1,1,lat) ! write(6,"('fs(:,k,1,2,lat)=',/,(6e12.4))") fs(:,k,1,2,lat) ! write(6,"('fs(:,k,2,1,lat)=',/,(6e12.4))") fs(:,k,2,1,lat) ! write(6,"('fs(:,k,2,2,lat)=',/,(6e12.4))") fs(:,k,2,2,lat) ! write(6,"('fs(:,k,1,0,lat)=',/,(6e12.4))") fs(:,k,1,0,lat) ! write(6,"('fs(:,k,2,0,lat)=',/,(6e12.4))") fs(:,k,2,0,lat) ! enddo ! k=lev0,lev1-1 call addfld('XNMBAR',' ',' ',xnmbar, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('POX1',' ',' ',pox1, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('POX2',' ',' ',pox2, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('LOX1',' ',' ',lox1, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('LOX2',' ',' ',lox2, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('LOX3',' ',' ',lox3, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('PO21',' ',' ',po21, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('PO22',' ',' ',po22, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('PO23',' ',' ',po23, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('LO21',' ',' ',lo21, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) call addfld('LO22',' ',' ',lo22, | 'lev',lev0,lev1,'lon',lon0,lon1,lat) end subroutine comp_o2o end module o2o_module