! module ox_module ! ! Calculate matrix coefficients in FS for partitioning of OX into O and O3. ! This must be called before the composition module (comp_major.F). ! use addfld_module,only: addfld implicit none contains !----------------------------------------------------------------------- subroutine comp_ox(tn,o2,ox,op,o1,o1d,o3,no,n4s,n2d,o2p,n2p,nplus, | nop,ne,xiop2p,xiop2d,barm,xnmbar,lev0,lev1,lon0,lon1,lat) use params_module,only: spval use cons_module,only: | rmassinv_o1 ,rmassinv_h , | rmassinv_no ,rmassinv_n4s,rmassinv_n2d, | rmassinv_o3 ,rmassinv_o2 ,rmassinv_n2 , | rmass_o1 ,rmass_o2 ,rmass_o3 use qrj_module,only: | pdo3d, ! photodissociation of o3(d) | pdo3p, ! photodissociation of o3(p) ! | pdo2, ! photodissociation of o2 (not yet available in tiegcm) | qo2p, ! o2+ ionization | qop ! o+ ionization ! | pdnosrb ! photodissociation of no by srb (not yet available in tiegcm) use chemrates_module,only: | rkm7a,rkm7b,rkm21,rkm24,rkm37,beta9,rk12,rk16,rk24,ra1,ra2,rk5, | beta3,beta6,rk10,rk4,rk7,rk23,rk27,beta1,beta2,rk1,rk8,rk11,rk3, | rkm20,rk6,rk9,rkm22,rkm23, | fs ! (fs(lon0:lon1,nlevp1,2,0:2,lat0:lat1) (defined by this routine) use init_module,only: istep ! ! xoxlb(nlonp4,nlat) bottom boundary of OX (bndry module data in bndry.F), ! is defined by this routine, then used by sub bndry_comp in bndry.F. ! use bndry_module,only: xoxlb ! (nlonp4,nlat) implicit none ! ! Input args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | tn, ! neutral temperature | o2, ! molecular oxygen | ox, ! ox from previous step | op, ! o+ ion | o1, ! atomic oxygen | o1d, ! o(1d) | o3, ! ozone | no, ! no | n4s, ! n(4s) | n2d, ! n(2d) | o2p, ! o2+ from elden | n2p, ! n2+ from elden | nplus, ! n+ from oplus | nop, ! no+ from elden | ne, ! ne from elden | xiop2p, ! O+(2p) from oplus module | xiop2d, ! O+(2d) from oplus module | barm, ! mbar | xnmbar ! p0e(-z), etc. ! ! Local: integer :: i,k,i0,i1,k0,k1 real,dimension(lev0:lev1,lon0-2:lon1+2) :: | h, ! h (h and h+ will be added to main code later) | hplus ! h+ real,dimension(lev0:lev1,lon0:lon1) :: | o3_o, ! o3/o ratio | o_ox, ! o/ox ratio | barmm, ! barm at midpoints | n2, ! n2 (mmr) | rmtru_ox, ! "true" mass of ox ! Photodissociation: | pdo3m, ! pdo3d+pdo3p at midpoints | pdo2m, ! pdo2 at midpoints ! Ionization: | qo2pm, ! qo2p at midpoints | qopm, ! qop at midpoints ! Production and loss rates: | ox_prod1, ! ox production | ox_prod2, ! ox production | ox_loss1, ! ox loss | ox_loss2, ! ox loss | o2_prod1, ! o2 production | o2_prod2, ! o2 production | o2_prod3, ! o2 production | o2_loss ! o2 loss ! ! Diagnostics for secondary history (see addfld calls): real,dimension(lev0:lev1,lon0:lon1) :: | fs11, ! fs(i,k,1,1,lat) | fs12, ! fs(i,k,1,2,lat) | fs21, ! fs(i,k,2,1,lat) | fs22, ! fs(i,k,2,2,lat) | fs10, ! fs(i,k,1,0,lat) | fs20 ! fs(i,k,2,0,lat) ! ! Local hydrogen is zero for now (s.a., compart.F): h = 0. hplus = 0. ! ! Init local diagnostics: fs11 = spval fs12 = spval fs21 = spval fs22 = spval fs10 = spval fs20 = spval ! i0=lon0 ; i1=lon1 ; k0=lev0 ; k1=lev1 ! ! Check some inputs: call addfld('NO_ox','NO_ox',' ',no(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('N4S_ox','N4S_ox',' ',n4s(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('N2D_ox','N2D_ox',' ',n2d(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('N2P_ox','N2P_ox',' ',n2p(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('XIOP2P','XIOP2P',' ',xiop2p(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('XIOP2D','XIOP2D',' ',xiop2d(:,i0:i1), | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('BETA1','BETA1',' ',beta1(k0:k1-1,:,lat), | 'lev',k0,k1-1,'lon',i0,i1,lat) call addfld('RK1','RK1',' ',rk1(k0:k1-1,:,lat), | 'lev',k0,k1-1,'lon',i0,i1,lat) call addfld('BARM_ox','BARM_ox',' ',barm(k0:k1-1,:), | 'lev',k0,k1-1,'lon',i0,i1,lat) ! ! Convert several fields from full-levels (interfaces) to half-levels ! (midpoints). Most of these are from the qrj module. ! do i=lon0,lon1 do k=lev0,lev1-1 pdo3m(k,i) = 0.5*(pdo3d (k,i,lat)+pdo3d(k+1,i,lat))+ | 0.5*(pdo3p (k,i,lat)+pdo3p(k+1,i,lat)) ! pdo2m not yet available: ! pdo2m(k,i) = 0.5*(pdo2 (k,i,lat)+pdo2 (k+1,i,lat)) pdo2m(k,i) = 0. qo2pm(k,i) = 0.5*(qo2p (k,i,lat)+qo2p (k+1,i,lat)) qopm (k,i) = 0.5*(qop (k,i,lat)+qop (k+1,i,lat)) barmm(k,i) = 0.5*(barm(k,i)+barm(k+1,i)) enddo ! k=lev0,lev1-1 pdo3m(lev1,i) = spval pdo2m(lev1,i) = spval qo2pm(lev1,i) = spval qopm (lev1,i) = spval barmm(lev1,i) = spval enddo ! i=lon0,lon1 call addfld('PDO3M','PDO3M',' ',pdo3m,'lev',k0,k1,'lon',i0,i1,lat) call addfld('QO2PM','QO2PM',' ',qo2pm,'lev',k0,k1,'lon',i0,i1,lat) call addfld('QOPM' ,'QOPM' ,' ',qopm ,'lev',k0,k1,'lon',i0,i1,lat) call addfld('BARMM','BARMM',' ',barmm,'lev',k0,k1,'lon',i0,i1,lat) do i=lon0,lon1 do k=lev0,lev1-1 ! ! Ratio O3/O: ! o3_o(k,i) = xnmbar(k,i)*rkm21(k,i,lat)*o2(k,i)*rmassinv_o2* | xnmbar(k,i)/barmm(k,i) / (xnmbar(k,i)* | (rkm24(k,i,lat)*o1 (k,i)*rmassinv_o1+ | (rkm7a+rkm7b) *o1d(k,i)*rmassinv_o1+ | rkm37(k,i,lat) *h (k,i)*rmassinv_h+ | beta9(k,i,lat) *no (k,i)*rmassinv_no)+ | pdo3m(k,i)) ! ! Ratio O/OX: o_ox(k,i) = 1./(1.+o3_o(k,i)) ! ! OX production terms: ! (pdnosrb not yet available in tiegcm) ! n2(k,i) = 1.-o2(k,i)-ox(k,i) ox_prod1(k,i) = xnmbar(k,i)**2* | (beta3(k,i,lat)*n4s(k,i)*rmassinv_n4s*no (k,i)*rmassinv_no + | beta6 *n2d(k,i)*rmassinv_n2d*no (k,i)*rmassinv_no)+ ! | 0.5*(pdnosrb(k,i,lat)+pdnosrb(k+1,i,lat))* ! | no(k,i)*rmassinv_no*xnmbar(k,i)+ | xnmbar(k,i)*(rk4*o2p(k,i)*n4s(k,i)*rmassinv_n4s | +rk10*op(k,i)*n2d(k,i)*rmassinv_n2d+ | rk12*op(k,i)*h(k,i)*rmassinv_h+ | (rk16*xiop2p(k,i)+rk24*xiop2d(k,i))*n2(k,i)*rmassinv_n2)+ | (ra1(k,i,lat)*nop(k,i)+2.*ra2(k,i,lat)*o2p(k,i))* | sqrt(ne(k,i)*ne(k+1,i)) ! ox_prod2(k,i) = xnmbar(k,i)* | (beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+ | beta2*n2d(k,i)*rmassinv_n2d)+ | rk1(k,i,lat)*op(k,i)+rk7*nplus(k,i)+ | rk23*xiop2p(k,i)+rk27*xiop2d(k,i)+ | 2.*pdo2m(k,i) ! ! OX loss terms: ! ox_loss1(k,i) = | (2.*rkm20(k,i,lat)*xnmbar(k,i)/barmm(k,i)+ | 2.*rkm24(k,i,lat)*o3_o(k,i))*o_ox(k,i)**2 ! ox_loss2(k,i) = (xnmbar(k,i)* | (rkm37(k,i,lat)*h (k,i)*rmassinv_h *o3_o(k,i)+ | beta9 (k,i,lat)*no (k,i)*rmassinv_no *o3_o(k,i))+ | rk3(k,i,lat)*n2p(k,i) +rk8 *nplus(k,i)+ | rk11(k,i,lat)*hplus(k,i))*o_ox(k,i) ! ! O2 production: ! o2_prod1(k,i) = | (rkm20(k,i,lat)*xnmbar(k,i)/barmm(k,i)+ | 2.*rkm24(k,i,lat)*o3_o(k,i))*o_ox(k,i)**2 ! o2_prod2(k,i) = | ((rkm37(k,i,lat)*h (k,i)*rmassinv_h *o3_o(k,i)+ | beta9(k,i,lat)*no (k,i)*rmassinv_no *o3_o(k,i))* | xnmbar(k,i)+pdo3m(k,i)*o3_o(k,i))*o_ox(k,i) ! o2_prod3(k,i) = | rk5*no(k,i)*rmassinv_no*o2p(k,i)*xnmbar(k,i) ! ! O2 loss: ! o2_loss(k,i) = xnmbar(k,i)* | (beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+ | beta2 *n2d(k,i)*rmassinv_n2d)+ | rk1(k,i,lat)*op(k,i)+(rk6+rk7)*nplus(k,i)+ | rk9*n2p(k,i)+pdo2m(k,i)+rk23*xiop2p(k,i)+rk27*xiop2d(k,i) ! ! "True" OX mass: ! rmtru_ox(k,i) = o_ox(k,i)*(rmass_o1+o3_o(k,i)*rmass_o3) ! write(6,"('comp_ox: lat=',i3,' k=',i3,' i=',i3, ! | ' o3_o=',e12.4,' o_ox=',e12.4,' ox_prod1=',e12.4, ! | ' ox_prod2=',e12.4,' ox_loss1=',e12.4,' ox_loss2=',e12.4, ! | /,' o2_prod1,2,3=',3e12.4,' o2_loss=',e12.4,' rmtru=', ! | e12.4)") lat,k,i,o3_o(k,i),o_ox(k,i),ox_prod1(k,i), ! | ox_prod2(k,i),ox_loss1(k,i),ox_loss2(k,i),o2_prod1(k,i), ! | o2_prod2(k,i),o2_prod3(k,i),o2_loss(k,i),rmtru_ox(k,i) ! ! Matrix coefficients (fs(lon0:lon1,nlevp1,2,0:2,lat0:lat1)) ! fs(i,k,1,1,lat) = -o2_loss(k,i)-xnmbar(k,i)**2*o_ox(k,i)* | ox(k,i)/rmtru_ox(k,i)*(1./3.*rkm21(k,i,lat)*o_ox(k,i)* | ox(k,i)/rmtru_ox(k,i)+ 2./3.*rkm22(k,i,lat)* | o2(k,i)*rmassinv_o2 + 1./2.*rkm23(k,i,lat)* | n2(k,i)*rmassinv_n2) fs(i,k,1,2,lat) = (o2_prod2(k,i)+xnmbar(k,i)*o2_prod1(k,i)* | ox(k,i)/rmtru_ox(k,i)-xnmbar(k,i)**2*o_ox(k,i)* | o2(k,i)*rmassinv_o2*(2./3.*rkm21(k,i,lat)*o_ox(k,i)* | ox(k,i)/rmtru_ox(k,i)+1./3.*rkm22(k,i,lat)*o2(k,i)* | rmassinv_o2+1./2.*rkm23(k,i,lat)*n2(k,i)*rmassinv_n2))* | rmass_o2/rmtru_ox(k,i) fs(i,k,2,1,lat) = ox_prod2(k,i)*rmtru_ox(k,i)*rmassinv_o2 fs(i,k,2,2,lat) = -ox_loss2(k,i)-ox_loss1(k,i)*ox(k,i)/ | rmtru_ox(k,i)*xnmbar(k,i) fs(i,k,1,0,lat) = (o2_prod3(k,i)-qo2pm(k,i))*rmass_o2/ | xnmbar(k,i) fs(i,k,2,0,lat) = (ox_prod1(k,i)-qopm(k,i))*rmtru_ox(k,i)/ | xnmbar(k,i) ! ! Save fs to secondary history: fs11(k,i) = fs(i,k,1,1,lat) fs12(k,i) = fs(i,k,1,2,lat) fs21(k,i) = fs(i,k,2,1,lat) fs22(k,i) = fs(i,k,2,2,lat) fs10(k,i) = fs(i,k,1,0,lat) fs20(k,i) = fs(i,k,2,0,lat) enddo ! k=lev0,lev1-1 ! ! Upper boundaries (nlevp1) are not calculated: o3_o(lev1,i) = spval ox_prod1(lev1,i) = spval ox_prod2(lev1,i) = spval ox_loss1(lev1,i) = spval ox_loss2(lev1,i) = spval o2_prod1(lev1,i) = spval o2_prod2(lev1,i) = spval o2_prod3(lev1,i) = spval o2_loss (lev1,i) = spval rmtru_ox(lev1,i) = spval ! ! Solve quadratic for lower boundary of OX: xoxlb(i,lat) = 2.e-5 enddo ! i=lon0,lon1 call addfld('O3_O','O3_O',' ',o3_o,'lev',k0,k1,'lon',i0,i1,lat) call addfld('OX_PROD1','OX_PROD1',' ',ox_prod1, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('OX_PROD2','OX_PROD2',' ',ox_prod2, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('OX_LOSS1','OX_LOSS1',' ',ox_loss1, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('OX_LOSS2','OX_LOSS2',' ',ox_loss2, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('O2_PROD1','O2_PROD1',' ',o2_prod1, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('O2_PROD2','O2_PROD2',' ',o2_prod2, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('O2_PROD3','O2_PROD3',' ',o2_prod3, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('O2_LOSS' ,'O2_LOSS' ,' ',o2_loss, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('RMTRU_OX' ,'RMTRU_OX' ,' ',rmtru_ox, | 'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS11_ox',' ',' ',fs11,'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS12_ox',' ',' ',fs12,'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS21_ox',' ',' ',fs21,'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS22_ox',' ',' ',fs22,'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS10_ox',' ',' ',fs10,'lev',k0,k1,'lon',i0,i1,lat) call addfld('FS20_ox',' ',' ',fs20,'lev',k0,k1,'lon',i0,i1,lat) end subroutine comp_ox end module ox_module