!
      module radcool_module
!
! 1/27/06 btf: sub radcool replaces newton and newto3p for zp-10 lbc.
!
      contains
!-----------------------------------------------------------------------
      subroutine radcool(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,spval
      use cons_module,only: rmass_co2,p0,expz,expzmid_inv,boltz,
     |  avto,avo,rmassinv_o2,rmassinv_o1,rmassinv_n2,rmassinv_no,
     |  gask
      use init_module,only: zpmid
      use addfld_module,only: addfld
      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:
      integer :: k,i
      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)
     |  cool_co2,        ! co2 contribution to cooling
     |  cool_no,         ! no contribution to cooling
     |  cool_o3p,        ! o3p cooling
     |  o3p_implicit,    ! implicit o3p cooling (output by o3pcool)
     |  o3p_explicit,    ! explicit o3p cooling (output by o3pcool)
     |  cool_total       ! total cooling
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  vo2,vo,vn2
      real,dimension(lon0:lon1) :: cppf
      real,parameter :: ur=8.3144e+7
!
! Cooling in degrees per day for secondary history diagnostics:
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  cooltot_deg,   ! total cooling (deg/day)
     |  coolco2_deg,   ! 15 micron co2
     |  coolno_deg,    ! 5.3 micron no
     |  coolo3_deg,    ! 9.6 micron o3
     |  coolo3p_deg,   ! 63 micron o3+
     |  cool_total_deg ! 
!
! 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):
          cool_co2(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))
!
! 1/30/06 btf: Set cool_co2 to 1.e-20 below zp -8 (Roble, personal comm)
! 2/17/06 btf: changed cutoff from -8 to -7 (as per Roble modsrc)
!         if (zpmid(k) < -8.) cool_co2(k,i) = 1.e-20
          if (zpmid(k) < -7.) cool_co2(k,i) = 1.e-20

!         if (zpmid(k) < -8.) then
!           if (i==1) write(6,"('radcool: k=',i3,' cool_co2=',e12.4,
!    |        ' -- resetting to 1.e-20.')") k,cool_co2(k,i)
!           cool_co2(k,i) = 1.e-20
!         else
!           if (i==1) write(6,"('radcool: k=',i3,' cool_co2=',e12.4)")
!    |        k,cool_co2(k,i)
!         endif
!
! 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):
          cool_no(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.*cool_no(k,i)+960.*cool_co2(k,i))/
     |      (.5*(cp(k,i)+cp(k+1,i))*tn(k,i)**2)
          cool_explicit(k,i) = ((1.-2700./tn(k,i))*cool_no(k,i)+
     |      (1.-960./tn(k,i))*cool_co2(k,i))*expz(k)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Add o(3p) cooling to total implicit and explicit terms.
! Sub o3pcool also returns o3+ implicit and explicit terms:
!
      call o3pcool(tn,o1,cp,cool_implicit,cool_explicit,
     |  o3p_implicit,o3p_explicit,lev0,lev1,lon0,lon1,lat)
!
! Combine implicit and explicit terms for total and o3+ cooling:
      do i=lon0,lon1
        do k=lev0,lev1-1
          cool_total(k,i) = cool_explicit(k,i) / expz(k) +
     |                      cool_implicit(k,i)*cp(k,i)*tn(k,i)
          cool_o3p(k,i) = o3p_explicit(k,i) / expz(k) +
     |                    o3p_implicit(k,i)*cp(k,i)*tn(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      do i=lon0,lon1
        do k=lev0,lev1
          vo2 (k,i) = o2 (k,i)*barm(k,i)*rmassinv_o2
          vo  (k,i) = o1 (k,i)*barm(k,i)*rmassinv_o1
          vn2 (k,i) = (1.-o2(k,i)-o1(k,i))*barm(k,i)*rmassinv_n2
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
! Save cooling terms on secondary history (deg/day):
! cppf: conversion factor from ergs to deg/day
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          cppf(i) = 86400./(0.5*ur*(7.*rmassinv_o2*vo2(k,i)+
     |                              7.*rmassinv_n2*vn2(k,i)+
     |                              5.*rmassinv_o1*vo (k,i)))
          cooltot_deg(k,i) = cool_total(k,i)*cppf(i)
          coolco2_deg(k,i) = cool_co2(k,i)*cppf(i)
          coolno_deg (k,i) = cool_no (k,i)*cppf(i)
          coolo3p_deg(k,i) = cool_o3p(k,i)*cppf(i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
      cool_implicit(lev1,:) = spval
      cool_explicit(lev1,:) = spval
      cooltot_deg(lev1,:) = spval
      coolco2_deg(lev1,:) = spval
      coolno_deg (lev1,:) = spval
      coolo3p_deg(lev1,:) = spval
!
      call addfld('COOL_IMP' ,'Implicit Cooling',' ',
     |  cool_implicit(:,lon0:lon1),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('COOL_EXP' ,'Explicit Cooling',' ',
     |  cool_explicit(:,lon0:lon1),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('COOL_DEG' ,'Total Cooling','deg/day',
     |  cooltot_deg(:,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('CO2_DEG' ,'15 micron co2 cooling','deg/day',
     |  coolco2_deg(:,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('NO_DEG' ,'5.3 micron no cooling','deg/day',
     |  coolno_deg(:,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('O3P_DEG' ,'63 micron o3+ cooling','deg/day',
     |  coolo3p_deg(:,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)

      end subroutine radcool
!-----------------------------------------------------------------------
#include <defs.h>
      subroutine o3pcool(tn,o1,cp,cool_implicit,cool_explicit,
     |  o3p_implicit,o3p_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,nlev,spval
      use cons_module,only: avo,rmassinv_o1,expz
      use addfld_module,only: addfld
      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  
!
! Output args:
      real,dimension(lev0:lev1,lon0:lon1),intent(out) ::
     |  o3p_implicit,   ! o3+ implicit cooling
     |  o3p_explicit    ! o3+ explicit cooling
!
! Local:
      integer :: k,i
      real,dimension(lev0:lev1,lon0:lon1) :: cool_o3p
!
! 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 (NLEV==28)
! lbc=-7, dz=0.5, nlev=28, nlevp1=29:
        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/)
#elif (NLEV==56)
! lbc=-7, dz=0.25, nlev=56, nlevp1=57:
        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/)
#elif (NLEV==34)
! lbc=-10, dz=0.5, nlev=34, nlevp1=35:
        xfac =
     |    (/0.1000E-01, 0.1000E-01, 0.1000E-01, 0.1000E-01, 0.1000E-01, 
     |      0.1000E-01, 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/)
#elif (NLEV==68)
! 12/19/05 btf: adding double res with new lbc
! lbc=-10, dz=0.25, nlev=68, nlevp1=69:
        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, 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,"('>>> newto3p: unknown NLEV=',i3)") nlev
        call shutdown('newto3p NLEV')
#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)))
          o3p_implicit(k,i) = fac2(k,i)/(0.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)
          o3p_explicit(k,i) = fac1(k,i)*expz(k)
!
! Calculate o3p cooling:
          cool_o3p(k,i) = o3p_explicit(k,i) / expz(k) +
     |                    o3p_implicit(k,i)*cp(k,i)*tn(k,i)

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      end subroutine o3pcool
      end module radcool_module
