!
      subroutine newton(tn,o2,o1,no,barm,cp,cool_implicit,cool_explicit,
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculate implicit and explicit cooling terms.
!
      use params_module,only: dz
      use cons_module,only: rmass_co2,p0,expz,expzmid_inv,boltz,
     |  avto,avo,rmassinv_o2,rmassinv_o1,rmassinv_n2,rmassinv_no,
     |  gask
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  o2,      ! molecular oxygen (mmr)
     |  o1,      ! atomic oxygen (mmr)
     |  no,      ! nitric oxide (mmr)
     |  cp,      ! specific heat
     |  barm     ! mean molecular weight
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  cool_implicit, cool_explicit ! output
!
! Local:
      real,parameter :: co2u=350.e-6
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  xmco2,
     |  nco2,aco2,bco2,  ! n(co2), a(co2), b(co2) (s10,s11,s12)
     |  ano,             ! a(no) (s11)
     |  co2_cool,        ! co2 contribution to cooling
     |  no_cool          ! no contribution to cooling
      integer :: k,i
      integer :: nk,nkm1 ! for addfsech
!
! Integral(m(co2)*dz/mbar) (s10).
! Low boundary: 
      do i=lon0,lon1
        xmco2(1,i) = .25*rmass_co2*(1./barm(1,i)+1./barm(2,i))*dz
!
! Integrate upwards:
        do k=lev0+1,lev1-1
          xmco2(k,i) = xmco2(k-1,i)+rmass_co2/barm(k,i)*dz
        enddo ! k=lev0+1,lev1-1
!
! n(co2) (s10):
        do k=lev0,lev1-1
          nco2(k,i) = co2u*p0*expz(1)*expzmid_inv/(boltz*tn(k,i))*
     |      exp(-xmco2(k,i))
!
! a(co2) (s11):
          if (tn(k,i) >= 200.) then
            aco2(k,i) = 2.5e-15*(1.+0.03*(tn(k,i)-200.))
          else
            aco2(k,i) = 2.5e-15
          endif
!
! b(co2) (s12):
          if (tn(k,i) < 260.) then
             bco2(k,i) = 1.56e-12
          elseif (tn(k,i) >= 260..and.tn(k,i) <= 300.) then
             bco2(k,i) = (2.6-0.004*tn(k,i))*1.0e-12
          else
             bco2(k,i) = 1.4e-12
          endif
!
! Co2 cooling (s12):
          co2_cool(k,i) = 2.65e-13*nco2(k,i)*exp(-960./tn(k,i))*
     |      avo*((o2(k,i)*rmassinv_o2+(1.-o2(k,i)-o1(k,i))*rmassinv_n2)*
     |      aco2(k,i)+o1(k,i)*rmassinv_o1*bco2(k,i))
!
! a(no) (s11):
          ano(k,i) = 
     |      .5*(barm(k,i)+barm(k+1,i))*avo*p0*expz(k)/(gask*tn(k,i))*
     |      (4.2e-11*o1(k,i)*rmassinv_o1+
     |       2.4e-14*o2(k,i)*rmassinv_o2)
!
! No cooling (s11):
          no_cool(k,i) = 4.956e-12*(avo*no(k,i)*rmassinv_no)*
     |      (ano(k,i)/(ano(k,i)+13.3))*exp(-2700./tn(k,i))
!
! Total cooling terms.
          cool_implicit(k,i) = (2700.*no_cool(k,i)+960.*co2_cool(k,i))/
     |      (.5*(cp(k,i)+cp(k+1,i))*tn(k,i)**2)
          cool_explicit(k,i) = ((1.-2700./tn(k,i))*no_cool(k,i)+
     |      (1.-960./tn(k,i))*co2_cool(k,i))*expz(k)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      nk = lev1-lev0+1
      nkm1 = nk-1
!     call addfsech('XMCO2',' ',' ',xmco2(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!     call addfsech('CO2_COOL',' ',' ',co2_cool(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!     call addfsech('NO_COOL' ,' ',' ',no_cool(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!     call addfsech('COOL_IMP' ,' ',' ',cool_implicit(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
      call addfsech('COOL_EXP' ,' ',' ',cool_explicit(:,lon0:lon1),
     |  lon0,lon1,nk,nkm1,lat)
      end subroutine newton
!-----------------------------------------------------------------------
      subroutine newto3p(tn,o1,cp,cool_implicit,cool_explicit,
     |  lev0,lev1,lon0,lon1,lat)
!
! Add O(3P) cooling to implicit and explicit cooling terms.
! This is called from dynamics after newton.
!
      use params_module,only: dz
      use cons_module,only: avo,rmassinv_o1,expz
      implicit none
!
! 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)
     |  cp       ! specific heat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) ::
     |  cool_implicit,
     |  cool_explicit  
!
! Local:
      integer :: k,i
      integer :: nk,nkm1 ! for addfsech
!
! Change O(3P) cooling factor from 0.5 to 1.0 of the Bates expression.
      real,parameter :: 
     |  an(3) = (/0.835E-18,0.6,0.2/),
     |  bn(3) = (/228.,228.,325./)
      real :: xfac(lev0:lev1)
      real,dimension(lev0:lev1,lon0:lon1) :: fac1,fac2
!
      if (dz==0.5) then
        xfac =
     |    (/0.1000E-01, 0.1000E-01, 0.1000E-01, 0.5000E-01, 0.1000E+00,
     |      0.2000E+00, 0.4000E+00, 0.5500E+00, 0.7000E+00, 0.7500E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00/)
      elseif (dz==0.25) then
        xfac =
     |    (/0.1000E-01, 0.1000E-01, 0.1000E-01, 0.1000E-01, 0.1000E-01,
     |      0.3000E-01, 0.5000E-01, 0.7500E-01, 0.1000E+00, 0.1500E+00,
     |      0.2000E+00, 0.3000E+00, 0.4000E+00, 0.4750E+00, 0.5500E+00,
     |      0.6250E+00, 0.7000E+00, 0.7250E+00, 0.7500E+00, 0.7750E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00, 0.8000E+00,
     |      0.8000E+00, 0.8000E+00/)
      else
        write(6,"('>>> WARNING newto3p: unsupported dz=',f10.3)") dz
        call shutdown('newto3p')
      endif
!
      do i=lon0,lon1
        do k=lev0,lev1-1
!       
! fac1 = AN(1)*XFAC*N0*PSI2/M2*EXP(-BN(1)/T) (s11)
          fac1(k,i) = an(1)*.5*(xfac(k)+xfac(k+1))*avo*rmassinv_o1
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          fac1(k,i) = fac1(k,i)*o1(k,i)*exp(-bn(1)/tn(k,i)) ! s11
!
! fac2 = 1.+AN(2)*EXP(-BN(2)/T)+AN(3)*EXP(-BN(3)/T) (s12)

          fac2(k,i) = 1.+an(2)*exp(-bn(2)/tn(k,i))+
     |                   an(3)*exp(-bn(3)/tn(k,i))
          fac1(k,i) = fac1(k,i)/fac2(k,i)
!
! fac2 = DF/DT=S11/S12/T**2*(BN(1)+(BN(1)-BN(3))*AN(3)*EXP(-BN(3)/T)) (s12)
          fac2(k,i) = fac1(k,i)/fac2(k,i)/tn(k,i)**2*
     |      (bn(1)+(bn(1)-bn(3))*an(3)*exp(-bn(3)/tn(k,i)))
!
! Add implicit contributions to cool_implicit:
          cool_implicit(k,i) = cool_implicit(k,i)+fac2(k,i)/
     |      (.5*(cp(k,i)+cp(k+1,i)))
!
! Add explicit contributions to cool_explicit:
          fac1(k,i) = fac1(k,i)-tn(k,i)*fac2(k,i)
          cool_explicit(k,i) = cool_explicit(k,i)+fac1(k,i)*expz(k)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      nk = lev1-lev0+1
      nkm1 = nk-1
!     call addfsech('COOL_IMP' ,' ',' ',cool_implicit(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!     call addfsech('COOL_EXP' ,' ',' ',cool_explicit(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!
      end subroutine newto3p
