!
      module so_module
      use params_module,only: nlevp1,nlonp4,nlat,dz
      implicit none
! -------------------------------------------------------------
! Boundary conditions, production and loss for SO are defined
! by comp_so, and referenced by minor_so. Comp_so is called
! from a latitude loop in dynamics. After comp_so, dynamics calls
! minor_so, which passes this module data to sub minor. Sub 
! minor contains 3d mpi calls and its own latitude loops.
!
!     real,dimension(nlonp4,nlat)   :: so_ubc ! upper boundary
!     real,dimension(nlonp4,3,nlat) :: so_lbc ! lower boundary
!     real,dimension(nlevp1,nlonp4,nlat) ::
!    |  so_prod,     ! production of so
!    |  so_loss      ! loss of so
!
! Boundary conditions and production and loss terms are allocated
! subdomains by sub alloc_so (called from allocdata.F).
! - last modified by SWB:  04/29/13
! 04/29/13 swb: comp_so.F and minor_so.F structure created
! 04/29/13 swb: P&L reactions and rates to be added
! -------------------------------------------------------------
!
      real,allocatable,dimension(:,:)   :: so_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: so_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  so_prod,     ! production of so (k,i,j)
     |  so_loss      ! loss of so       (k,i,j)
!  Mutual diffusion coefficients (normalized) for O-SO, CO-SO, CO2-SO
!  Ask Ben: no m. diffusion below 115 knm where SOx reactions are important
      real :: phi_so(3) =  (/0., 0., 0./)
      real,parameter :: alfa_so = 0. ! thermal diffusion coefficient
!  Mass mixing ratio lower boundary condition for SO (at ~69 km)
!  massmix = mix*mso/mbar = mix*48./44.
!  (1)  mix1 = 3.0E-09 ! Zhang et al (2011) Figure 8 eyeball value
!  (2)  mix2 = 2.1E-09 ! KINETICS reference simulation
!     real :: pssob =  3.27E-09   !  From Zhang et al (2011) Figure 8 
      real :: pssob =  2.29E-09   !  KINETICS reference simulation
!
      contains
!-----------------------------------------------------------------------
      subroutine alloc_so(lon0,lon1,lat0,lat1)
!
! Allocate subdomains (without ghost cells) to module data for boundary 
! conditions and production and loss terms. This is called once per run 
! from sub allocdata (allocdata.F).
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate subdomains to boundary conditions:
      allocate(so_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_so: error allocating',
     |  ' so_ubc: stat=',i3)") istat
      allocate(so_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_so: error allocating',
     |  ' so_lbc: stat=',i3)") istat
!
! Allocate subdomains to production and loss:
      allocate(so_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_so: error allocating',
     |  ' so_prod: stat=',i3)") istat
      allocate(so_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_so: error allocating',
     |  ' so_loss: stat=',i3)") istat

!     write(6,"('alloc_so: allocated module data')")
!
      end subroutine alloc_so
!
!-----------------------------------------------------------------------
      subroutine comp_so(tn,o1,co,co2,n2,o2,so,so2,barm,
     |  lev0,lev1,lon0,lon1,lat)
!
!-----------------------------------------------------------------------
! Advance SO (minor): 
! Arguments for comp_so (14).
! ----------------------------------------------------------------------
! Advance so by one time step. This is called from driver at
! each subdomain latitude.
! Last modified: swb (05/29/13)
!              : add new rates from chemrates_sox
!              : finalized sosrc, soloss terms
!              : next interp trace species from KINETICS to VTGCM grid
! ----------------------------------------------------------------------
!
      use qrj_module,only: qjso, qjso2a, mk_xnmbari
      use fields_module,only: tlbc
      use cons_module,only: rmassinv_o1,rmassinv_co,rmassinv_n2,p0,
     |  rmass_so, rmass_so2, rmassinv_so, rmassinv_so2,
     |  expz,expzmid_inv,boltz,rmassinv_co2,rmassinv_o2,rmass_o1,
     |  rmass_co,rmass_o2,rmass_so,rmass_so2,gask,grav
      use chemrates_module,only:s114,s133,s166,s240,S244,s249,s258
      use addfld_module,only: addfld
!
!     write(6,"('enter comp_so: lat=',i2)") lat
! ----------------------------------------------------------------------
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  o1,      ! atomic oxygen (mmr)
     |  co,      ! carbon monoxide (mmr)
     |  co2,     ! carbon dioxide (mmr)
     |  n2,      ! molecular nitrogen (mmr)
     |  o2,      ! molecular oxygen (mmr)
     |  so,      ! so (mmr)
     |  so2,     ! so2 (mmr)
     |  barm     ! mean molecular weight
!
! ----------------------------------------------------------------------
! Local:
      integer :: k,i,k0,k1,i0,i1,nk  ! for addfld (ask for check)
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  tni,            ! TN at interfaces
     |  qjsoi,          ! Jso at interfaces (qrj module): 1/sec units from #/cm3/sec
     |  qjso2ai,        ! Jso2a->SO+O at interfaces (qrj module): 1/sec units from #/cm3/sec
     |  xnmbari,        ! xnmbar at interfaces (qrj module)
     |  xnmbar,         ! xnmbar at midlevels (qrj module)
     |  xtot,           ! Total number density (#/cm3)
     |  xo,             ! O1 number density (#/cm3)
     |  xco,            ! CO number density (#/cm3)
     |  xco2,           ! CO2 number density (#/cm3)
     |  xn2,            ! N2 number density (#/cm3)
     |  xo2,            ! O2 number density (#/cm3)
     |  xso,            ! SO number density (#/cm3)
     |  xso2,           ! SO2 number density (#/cm3)
     |  soprod,         ! SO production (#/cm3/s)
     |  soloss1,        ! SO loss1 (#/cm3/s)
     |  soloss2,        ! SO loss2 (#/cm3/s)
     |  xcl,            ! CL (#/cm3) 
     |  xcl2,           ! CL2 (#/cm3) 
     |  xclo,           ! CLO (#/cm3) 
     |  xclco,          ! CLCO (#/cm3) 
     |  xclco3,         ! CLCO3 (#/cm3) 
     |  xs,             ! S (#/cm3) 
     |  xso3,           ! SO3 (#/cm3) 
     |  xh,             ! H (#/cm3) 
     |  xh2,            ! H2 (#/cm3)
     |  xoh,            ! OH (#/cm3) 
     |  xho2,           ! HO2 (#.cm3) 
     |  xo3,            ! O3 (#/cm3) 
     |  xhcl            ! hcl (#/cm3)
!
! From recent 1-D KINETICS calculation for nightside
! 29-MAY-2013 (private communication, C. Parkinson)
      real,dimension(lev0:lev1) ::
     |  cl,            ! CL (vmr) 
     |  cl2,           ! CL2 (vmr) 
     |  clo,           ! CLO (vmr) 
     |  clco,          ! CLCO (vmr) 
     |  clco3,         ! CLCO3 (vmr) 
     |  so3,           ! SO3 (vmr) 
     |  s,             ! S (vmr) 
     |  h,             ! H (vmr) 
     |  h2,            ! H2 (vmr)
     |  oh,            ! OH (vmr) 
     |  ho2,           ! HO2 (vmr) 
     |  o3,            ! O3 (vmr) 
     |  hcl            ! hcl (vmr)
!
! ----------------------------------------------------------------------
! For addfld calls: Check what is needed!
      k0=lev0
      k1=lev1
      i0=lon0
      i1=lon1
      nk = lev1-1
! *******************************************************************
! Low Altitude Chemistry from Chris Parkinson Venus KINETICS code
! Nightside profiles in VMR: 22 (values) + 47 (zeros) = 69 levels
!  (30-MAY-2013 Update to VMR for 13 species)
! Alt rage = ~68 - ~112km (zp = -16 to -5.5)
!********************************************************************
      cl(:) = (/
     | 1.286E-08, 1.798E-08, 2.365E-08, 2.923E-08, 3.368E-08, 3.566E-08,
     | 3.478E-08, 3.343E-08, 3.557E-08, 3.900E-08, 4.136E-08, 4.611E-08,
     | 5.683E-08, 7.439E-08, 1.022E-07, 1.266E-07, 1.527E-07, 1.732E-07,
     | 1.904E-07, 2.049E-07, 2.174E-07, 2.283E-07, 0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      cl2(:) = (/
     | 1.325E-08, 1.119E-08, 8.692E-09, 6.115E-09, 3.924E-09, 2.448E-09,
     | 1.395E-09, 7.241E-10, 3.837E-10, 2.051E-10, 9.209E-11, 4.165E-11,
     | 2.230E-11, 1.295E-11, 6.294E-12, 2.714E-12, 1.283E-12, 6.303E-13,
     | 2.832E-13, 1.217E-13, 4.915E-14, 2.037E-14, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      clo(:) = (/
     | 5.191E-11, 6.758E-11, 1.026E-10, 1.467E-10, 2.420E-10, 8.057E-10,
     | 2.309E-09, 3.010E-09, 1.523E-09, 1.527E-10, 2.993E-11, 1.790E-11,
     | 1.154E-11, 7.656E-12, 4.651E-12, 2.817E-12, 1.581E-12, 8.268E-13,
     | 3.914E-13, 1.777E-13, 7.843E-14, 3.567E-14, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      clco(:) = (/
     | 8.730E-12, 1.213E-11, 1.784E-11, 2.989E-11, 6.085E-11, 1.677E-10,
     | 4.096E-10, 6.221E-10, 5.948E-10, 3.412E-10, 8.294E-11, 3.258E-11,
     | 1.507E-11, 6.807E-12, 2.572E-12, 1.007E-12, 4.311E-13, 2.090E-13,
     | 1.026E-13, 5.269E-14, 2.755E-14, 1.484E-14, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      clco3(:) = (/
     | 2.198E-11, 2.045E-11, 2.747E-11, 4.997E-11, 1.196E-10, 4.374E-10,
     | 1.550E-09, 3.034E-09, 2.623E-09, 8.529E-10, 3.115E-11, 4.126E-12,
     | 7.018E-13, 1.026E-13, 1.073E-14, 1.223E-15, 1.608E-16, 3.078E-17,
     | 5.836E-18, 1.200E-18, 2.546E-19, 5.816E-20, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      h(:) = (/
     | 1.061E-13, 3.074E-13, 6.720E-13, 1.371E-12, 2.602E-12, 4.251E-12,
     | 6.333E-12, 1.004E-11, 3.831E-11, 5.822E-10, 3.793E-09, 1.445E-08,
     | 3.849E-08, 8.110E-08, 1.611E-07, 2.433E-07, 3.491E-07, 4.466E-07,
     | 5.364E-07, 6.128E-07, 6.738E-07, 7.157E-07, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      h2(:) = (/
     | 1.287E-11, 2.901E-11, 5.708E-11, 1.397E-10, 3.857E-10, 8.965E-10,
     | 1.672E-09, 2.667E-09, 3.899E-09, 5.399E-09, 7.296E-09, 9.767E-09,
     | 1.272E-08, 1.580E-08, 1.940E-08, 2.221E-08, 2.486E-08, 2.685E-08,
     | 2.875E-08, 3.081E-08, 3.329E-08, 3.634E-08, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      oh(:) = (/
     | 4.890E-15, 5.446E-15, 8.956E-15, 1.976E-14, 4.176E-14, 6.652E-14,
     | 1.079E-13, 2.043E-13, 8.698E-13, 5.159E-12, 1.572E-11, 2.906E-11,
     | 3.619E-11, 3.185E-11, 2.142E-11, 1.391E-11, 8.243E-12, 4.586E-12,
     | 2.389E-12, 1.263E-12, 7.513E-13, 6.064E-13, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      ho2(:) = (/
     | 2.135E-14, 3.644E-14, 6.299E-14, 1.189E-13, 2.215E-13, 4.003E-13,
     | 7.293E-13, 1.309E-12, 2.256E-12, 4.252E-12, 8.215E-12, 1.287E-11,
     | 1.414E-11, 1.091E-11, 6.178E-12, 3.371E-12, 1.664E-12, 8.168E-13,
     | 3.859E-13, 1.825E-13, 8.552E-14, 4.049E-14, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      o3(:) = (/
     | 2.250E-11, 7.002E-12, 5.795E-12, 9.205E-12, 1.532E-11, 1.734E-11,
     | 2.097E-11, 4.133E-11, 2.288E-10, 1.575E-09, 4.075E-09, 6.283E-09,
     | 8.300E-09, 1.060E-08, 1.304E-08, 1.330E-08, 1.304E-08, 1.061E-08,
     | 7.141E-09, 4.196E-09, 2.189E-09, 1.070E-09, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      hcl(:) = (/
     | 3.788E-07, 3.778E-07, 3.771E-07, 3.766E-07, 3.764E-07, 3.763E-07,
     | 3.765E-07, 3.767E-07, 3.772E-07, 3.775E-07, 3.766E-07, 3.721E-07,
     | 3.614E-07, 3.439E-07, 3.162E-07, 2.919E-07, 2.659E-07, 2.455E-07,
     | 2.284E-07, 2.142E-07, 2.020E-07, 1.915E-07, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      s(:) = (/
     | 3.000E-13, 3.088E-13, 1.910E-13, 1.366E-13, 1.309E-13, 9.668E-14,
     | 6.398E-14, 7.870E-14, 4.326E-13, 1.345E-12, 2.612E-12, 4.512E-12,
     | 8.457E-12, 1.642E-11, 3.708E-11, 6.248E-11, 1.053E-10, 1.800E-10,
     | 3.178E-10, 5.583E-10, 9.727E-10, 1.655E-09, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
      so3(:) = (/
     | 6.769E-09, 3.535E-10, 2.492E-11, 6.631E-12, 4.240E-12, 3.696E-12,
     | 4.278E-12, 4.037E-12, 2.244E-12, 5.432E-13, 1.014E-12, 1.392E-11,
     | 2.556E-10, 2.682E-09, 2.163E-08, 2.303E-08, 1.431E-08, 8.000E-09,
     | 4.234E-09, 2.331E-09, 1.437E-09, 1.137E-09, 0.0,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |  0.,0.,0.,0.,0.,0./)
!
! ----------------------------------------------------------------------
! Species Number densities at interfaces (#/cm3)
! Given Qrj dissociation rates (#/cm3/sec) at midlevels from qrj
!
! ----------------------------------------------------------------------
! Lower boundary:
! so_lbc(:,1)=A, so_lbc(:,2)=B, so_lbc(:,3)=C define lower boundary 
!   condition where A*DPSX/DZ + B*PSX + C = 0.
!
      do i=lon0,lon1
!
! Value at bottom given by specified mass missing ratio (start)
        so_lbc(i,1,lat) = 0.
        so_lbc(i,2,lat) = 1.
        so_lbc(i,3,lat) = -pssob
!
! Upper boundary: Zero diffusive flux at top
!       so_ubc(i,lat) = 0.
      enddo ! i=lon0,lon1
! --------------------------------------------------------------------
! Calculate tn at interfaces:
      do i=lon0,lon1
      tni(lev0,i) = tlbc(i,lat)
      	do k=lev0+1,lev1-1
         	tni(k,i) = .5*(tn(k-1,i)+tn(k,i))
      	enddo
      tni(lev1,i) = tn(lev1-1,i) ! nlevp1 <- nlev
      enddo
! --------------------------------------------------------------------
!
! Calculate p0*e(-z)*barm/kT (1/(N*MBAR))at interfaces
!     subroutine mk_xnmbari(tni,barm,xnmbari,lev0,lev1,lon0,lon1,lat)
!
      call mk_xnmbari(tni,barm(:,lon0:lon1),
     |  xnmbari,
     |  lev0,lev1,lon0,lon1,lat)
! ----------------------------------------------------------------------
! Species Number densities at interfaces (#/cm3)
! Given Qrj dissociation rates (#/cm3/sec) at midlevels from qrj
! Calculate Qrj dissociation rates (#/sec) at interfaces: Jso, Jso2
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          xo(k,i)   = xnmbari(k,i)*o1(k,i)*rmassinv_o1
          xco(k,i)  = xnmbari(k,i)*co(k,i)*rmassinv_co
          xco2(k,i) = xnmbari(k,i)*co2(k,i)*rmassinv_co2
          xo2(k,i) = xnmbari(k,i)*o2(k,i)*rmassinv_o2
	  if (xco2(k,i) < 1.0) xco2(k,i) = 1.0
          xn2(k,i)  = xnmbari(k,i)*n2(k,i)*rmassinv_n2
          xso(k,i)  = xnmbari(k,i)*so(k,i)*rmassinv_so
          xso2(k,i)  = xnmbari(k,i)*so2(k,i)*rmassinv_so2
          qjsoi(k,i) = 0.5*(qjso(k,i,lat)+qjso(k+1,i,lat))/xso(k,i)
          qjso2ai(k,i) = 0.5*(qjso2a(k,i,lat)+qjso2a(k+1,i,lat))
     |                   /xso2(k,i)
!
!   LOG10 DENSITIES FOR PLOTTING 
!         xco2l(k,i)=log10(xco2(k,i))
!         xcol(k,i)=log10(xco(k,i))
!         xn2l(k,i)=log10(xo2(k,i))
!         xol(k,i)=log10(xo(k,i))

!  Calcuklate Total Density
          xtot(k,i) = xo(k,i)+xco(k,i)+xco2(k,i)+xn2(k,i)
!
!  Converting vmr to number density
          xcl(k,i)   = xtot(k,i)*cl(k)
          xcl2(k,i)  = xtot(k,i)*cl2(k)
          xclo(k,i)  = xtot(k,i)*clo(k)
          xclco(k,i) = xtot(k,i)*clco(k)
          xclco3(k,i) = xtot(k,i)*clco3(k)
          xs(k,i) = xtot(k,i)*s(k)
          xso3(k,i) = xtot(k,i)*so3(k)
          xh(k,i)   = xtot(k,i)*h(k)
          xh2(k,i)  = xtot(k,i)*h2(k)
          xoh(k,i)   = xtot(k,i)*oh(k)
          xho2(k,i)  = xtot(k,i)*ho2(k)
          xo3(k,i) = xtot(k,i)*o3(k)
          xhcl(k,i) = xtot(k,i)*hcl(k)

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! ----------------------------------------------------------------------
! (1)  Major Sources and Losses for SO (before NOX Chemistry Added)
!      From Zhang et al (2012)  Figures 12a-b.
! ---------------------------------------------------------------------
      do i=lon0,lon1
        do k=lev0,lev1-1
! ---------------------------------------------------------------------
!     SOURCES (#/cm3.sec) AND SINKS(#/sec): Prognostics
!     SOURCES (#/cm3.sec) AND SINKS(#/cm3/sec): Diagnostics
! ---------------------------------------------------------------------
!
          so_prod(k,i,lat) = qjso2ai(k,i)*xso2(k,i) + 
     |       xs(k,i)*xo2(k,i)*s166                  

          write(6,"('comp_so: k,i,lat=',3i4,' qjso2ai(k,i)=',e12.4,
     |      ' xso2(k,i)=',e12.4,' xs(k,i)=',e12.4,' xo2(k,i)=',e12.4,
     |      ' s166=',e12.4,' so_prod(k,i,lat)=',e12.4)") k,i,lat,
     |      qjso2ai(k,i),xso2(k,i),xs(k,i),xo2(k,i),s166,
     |      so_prod(k,i,lat)
!
          soprod(k,i) = so_prod(k,i,lat)
! --------------------------------------
!
	  so_loss(k,i,lat) = - (qjsoi(k,i) + xo(k,i)*xco2(k,i)*s240 +
     |       xso(k,i)*xco2(k,i)*s249                                +
     |       xclo(k,i)*s133                                         +
     |       xclco3(k,i)*s244)
!
	  soloss1(k,i) = qjsoi(k,i)*xso(k,i)
	  soloss2(k,i) = xo(k,i)*xso(k,i)*xco2(k,i)*s240            +
     |                   xso(k,i)*xso(k,i)*xco2(k,i)*s249           +
     |                   xclo(k,i)*xso(k,i)*s133                    +
     |                   xclco3(k,i)*xso(k,i)*s244
!
         enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
! ----------------------------------------------------------------------
!   * O2 Production and Losses

      write(6,"('comp_so call addfld for SOPROD: lat=',i4,
     |  ' soprod mn,mx=',2e12.4)") lat,minval(soprod),maxval(soprod)

      call addfld('SOPROD',' ',' ',soprod,'lev',k0,k1,'lon',i0,i1,
     |  lat)
      call addfld('SOLOSS1',' ',' ',soloss1,'lev',k0,k1,'lon',i0,i1,
     |  lat)
      call addfld('SOLOSS2',' ',' ',soloss2,'lev',k0,k1,'lon',i0,i1,
     |  lat)
!
      end subroutine comp_so
!-----------------------------------------------------------------------
      subroutine minor_so(tn,difk,o1,co,co2,n2,so,so_nm,so_out,
     |  so_nm_out,lev0,lev1,lon0,lon1,lat0,lat1)

!     Arguments (16): o1, co, co2, n2 needed by minor.F for barm, etc.
      use cons_module,only: rmass_so, rmass_so2
!
! -----------------------------------------------------------------------
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  difk,    ! eddy diffusion
     |  o1,      ! atomic oxygen (mmr)
     |  co,      ! carbon monoxide (mmr)
     |  co2,     ! carbon dioxide (mmr)
     |  n2,      ! molecular nitrogen (mmr)
     |  so,      ! so (mmr)
     |  so_nm    ! so at time n-1 
!
! -----------------------------------------------------------------------
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  so_out,    ! SO output 
     |  so_nm_out  ! SO output at time n-1
!
! -----------------------------------------------------------------------
! Local:
      integer :: lat
      integer :: i0,i1,nk ! for addfld
!
!     write(6,"('enter minor_so')")
      i0 = lon0
      i1 = lon1
      nk = lev1-1
! -----------------------------------------------------------------------
! Minor returns so_out and so_nm_out. Module data so_prod,
! so_loss, etc, were defined by comp_so.  
!
!     subroutine minor(tn,difk,o2,o1,fcomp,fcomp_tm1,fcomp_out,
!    |  fcomp_tm1_out,sloss,sprod,flbc,fubc,rmx,phix,alfax,
!    |  lev0,lev1,lon0,lon1,lat0,lat1,idebug)
!
      call minor(tn,difk,o1,co,co2,n2,so,so_nm,so_out,so_nm_out,
     |  so_loss,so_prod,so_lbc,so_ubc,rmass_so,phi_so,
     |  alfa_so,lev0,lev1,lon0,lon1,lat0,lat1,0)

!
!       call addfld('SO_OUT' ,' ',' ',so_out(:,i0:i1,lat),
!    |    i0,i1,k0,k1,lat)
!       call addfld('SO_TM1' ,' ',' ',so_nm_out(:,i0:i1,lat),
!    |    i0,i1,k0,k1,lat)
      end subroutine minor_so
      end module so_module
