!
      module radcool_module
!
! 1/27/06 btf: sub radcool replaces newton and newto3p for zp-10 lbc.
! 7/06 btf: Add recur and vicool for 15 micron co2 cooling.
!
      use params_module,only: dz,spval,zibot
      use init_module,only: zpmid,zpint
      implicit none
      real,parameter ::
     |  zpm9   = -9.,  ! zp -9.0 level
     |  zpm5   = -5.,  ! zp -5.0 level
     |  zpsrf  = -21.5 ! zp at surface
      integer ::
     |  kzpm9,        ! k index to zp -9
     |  kzpm5         ! k index to zp -5
      integer,parameter ::
     |  lmax=59,      ! max number of wavelengths (from timegcm)
     |  lam0=36,      ! wavelength index used for zp -10.75
     |  lam1=43       ! wavelength index used for zp -5
      contains
!-----------------------------------------------------------------------
      subroutine radcool(tn,o2,o1,o3,no,barm,xnmbar,cp,cool_implicit,
     |  cool_explicit,nzpm9m5,lev0,lev1,lon0,lon1,lat)
!
! Calculate implicit and explicit cooling terms.
!
      use cons_module,only: rmass_co2,p0,expz,expzmid_inv,boltz,
     |  avto,avo,rmassinv_o2,rmassinv_o1,rmassinv_n2,rmassinv_no,
     |  gask,grav,expzmid,rmassinv_o3,rmassinv_co2
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: nzpm9m5,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)
     |  o3,      ! ozone (mmr)
     |  no,      ! nitric oxide (mmr)
     |  cp,      ! specific heat
     |  barm,    ! mean molecular weight
     |  xnmbar   ! p0e(-z)*barm/(boltz*T)
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  cool_implicit, cool_explicit ! output
!
! Local:
      integer :: k,kk,i,lam
      real,parameter :: co2u=370.e-6
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  xmco2,
     |  aco2,bco2,       ! a(co2), b(co2) (s10,s11,s12)
     |  co2_cm3,         ! co2 (cm3)
     |  co2_mmr,         ! co2 (mmr)
     |  colco2,
     |  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,vco2,vo3
      real,dimension(lon0:lon1) :: cppf
      real,parameter :: ur=8.3144e+7
      real :: factor,dfactor
!
      real,dimension(nzpm9m5,lon0:lon1) ::
     |  zuco2,   ! zp -9 to zp -5 (recur input)
     |  alrecur, ! zp -9 to zp -5 (recur output)
     |  zo3  ,   ! zp -9 to zp -5
     |  zsn2 ,   ! zp -9 to zp -5
     |  zso2 ,   ! zp -9 to zp -5
     |  zo   ,   ! zp -9 to zp -5
     |  zco2 ,   ! zp -9 to zp -5
     |  am   ,   ! zp -9 to zp -5
     |  zden     ! zp -9 to zp -5
!
! vicool output:
      real ::
     |  zhco2(lmax,lon0:lon1),
     |  zho3 (lmax,lon0:lon1),
     |  zflux(lon0:lon1)
!
! 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 ! 
!
! 7/24/06 btf: This new version of radcool for tiegcm assumes
! double vertical res (dz==0.25)
      if (dz /= 0.25) then
        write(6,"(/,'>>> radcool works only with dz 0.25: dz=',f8.2)")
     |    dz
        call shutdown('radcool')
      endif
!
! 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
          dfactor = p0*expz(lev0)*expzmid_inv/(boltz*tn(k,i))
          co2_cm3(k,i) = co2u*dfactor*exp(-xmco2(k,i))
          co2_mmr(k,i) = co2_cm3(k,i)*rmass_co2/xnmbar(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 to space approximation:
          cool_co2(k,i) = 2.65e-13*co2_cm3(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
!
! This can be overwritten after vicool below:
          if (zpmid(k) < -7.) cool_co2(k,i) = 1.e-20
!
! 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.
! 7/06 btf: this moved to after vicool call, to use new cool_co2 below zp -5 (or -7)
!         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
        co2_mmr(lev1,i) = co2_mmr(lev1-1,i)
        co2_cm3(lev1,i) = co2_cm3(lev1-1,i)
      enddo ! i=lon0,lon1

      call addfld('CO2_CM3' ,' ',' ',
     |  co2_cm3(lev0:lev1,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('CO2_MMR' ,' ',' ',
     |  co2_mmr(lev0:lev1,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      call addfld('CO2' ,' ',' ',
     |  co2_mmr(lev0:lev1,:),'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Upper boundary:
      factor = avo*p0*expz(lev1-1)*expzmid/(rmass_co2**2*grav)
      do i=lon0,lon1
        colco2(lev1,i) = factor*.5*(co2_mmr(lev1-1,i)+co2_mmr(lev1,i))*
     |    barm(lev1,i)
      enddo
!
! Integrate down to level 1:
! (loop order must be reversed from normal, i.e., i loop on inside):
!
      do k=lev1-1,lev0,-1
        factor = avo*p0*expz(k)/(rmass_co2*grav)*dz
        do i=lon0,lon1
          colco2(k,i) = colco2(k+1,i)+factor*co2_mmr(k,i)
        enddo ! i=lon0,lon1
!       write(6,"('radcool: lat=',i3,' k=',i3,' colco2(k,:)=',/,
!    |    (6e12.4))") lat,k,colco2(k,:)
      enddo
!
      kzpm9 = int((zpm9-zibot)/dz)+1   ! col index at zp -9 from model lb
      kzpm5 = int((zpm5-zibot)/dz)+1   ! col index at zp -5 from model lb
!
      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
          vo3 (k,i) = o3(k,i)*barm(k,i)*rmassinv_o3
          vn2 (k,i) = (1.-o2(k,i)-o1(k,i))*barm(k,i)*rmassinv_n2
          vco2(k,i) = co2_mmr(k,i)*barm(k,i)*rmassinv_co2
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
! Prepare inputs to recur and vicool (zp -9 to zp -5):
      do i=lon0,lon1
        do k=kzpm9,kzpm5 ! 5-21 for dz=0.25 model
          kk = k-kzpm9+1 ! 1-17 for dz=0.25 model
          zso2 (kk,i) = vo2   (k,i)
          zo   (kk,i) = vo    (k,i)
          zsn2 (kk,i) = vn2   (k,i)
          zco2 (kk,i) = vco2  (k,i)
          zuco2(kk,i) = colco2(k,i)                ! input for recur
          zden (kk,i) = xnmbar(k,i)/barm(k,i)
          am   (kk,i) = barm  (k,i)

        enddo ! k=kzpm9,kzpm5
      enddo ! i=lon0,lon1
!
! recur: al is output, zuco2 is input. Both are dimensioned (nzpm9m5,lon0:lon1)
      call recur(alrecur,zuco2,nzpm9m5,lev0,lev1,lon0,lon1,lat)

!     do k=1,nzpm9m5
!       write(6,"('radcool: lat=',i3,' k=',i3,' alrecur min,max=',
!    |    2e12.4)") lat,k,minval(alrecur(k,:)),maxval(alrecur(k,:))
!     enddo

!                        36:43               44:59
! vicool returns zhco2(lam0:lam1,:) and (lam1+1:lmax,:)
! This corresponds to zp -10.75 to -5.0
      call vicool(
     |  tn(:,lon0:lon1),zsn2,zso2,zo,zco2,alrecur,am,zden, ! input to vicool
     |  zflux, zhco2, zho3,                             ! output from vicool
     |  nzpm9m5,lev0,lev1,lon0,lon1,lat)

!     write(6,"(/,'radcool: lat=',i3,' zhco2(36:43,lon0)=',/,(6e12.4))")
!    |  lat,zhco2(36:43,lon0)
!     write(6,"('radcool: lat=',i3,' zhco2(36:43,lon0)=',/,(6e12.4))")
!    |  lat,zhco2(44:59,lon0)

!     do lam=lam0,lmax
!       write(6,"('radcool: lat=',i3,' lam=',i3,' zuco2 min,max=',
!    |    2e12.4)") lat,lam,minval(zuco2(lam,:)),maxval(zuco2(lam,:))
!     enddo
!
! Overwrite cool_co2 from bottom up to zp -5, with vicool output:
! (above zp -5, cool_co2 remains cool-to-space approx, see above)
      do i=lon0,lon1
        do k=1,kzpm5       ! zp -10 to -5
          kk = lam0+3+k-1  ! 39 -> 59 (zp -10 -> -5)
          cool_co2(k,i) = -zhco2(kk,i)

!         write(6,"('radcool: k=',i3,' kk=',i3,' i=',i3,' lat=',i3,
!    |      ' zhco2(kk,i)=',e12.4,' zp(k)=',f6.2)") k,kk,i,lat,
!    |      zhco2(kk,i),float(k-1)/4.-10.

        enddo ! k=kzpm9,kzpm5
      enddo ! i=lon0,lon1
!
      do i=lon0,lon1
	do k=1,kzpm9   ! zp -10 to -9
	  cool_co2(k,i) = cool_co2(kzpm9,i)
	enddo ! k=1,kzpm9
      enddo ! i=lon0,lon1
!
!
      call addfld('COOL_CO2' ,' ',' ',
     |  cool_co2(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!
! Total cooling terms.
      do i=lon0,lon1
        do k=lev0,lev1-1
          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
!
! 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
!-----------------------------------------------------------------------
      subroutine recur(al,zuco2,nzpm9m5,lev0,lev1,lon0,lon1,lat)
      use cco2gr_module,only:
     |  cor150, cor360, cor540, cor720, uco2co, uco2ro, alo
      implicit none
!
! Args:
      integer,intent(in) :: nzpm9m5,lev0,lev1,lon0,lon1,lat
      real,intent(in) ::
     |  zuco2(nzpm9m5,lon0:lon1)    ! zp -9 to zp -5
      real,intent(out) ::
     |  al   (nzpm9m5,lon0:lon1)    ! zp -9 to zp -5
!
! Local:
      integer :: k,i,ii,nlons
      real :: co2int(4),uref(4)
      real,dimension(lon1-lon0+1) :: a, cor, uc
 
!     do k=1,nzpm9m5
!       write(6,"('recur: lat=',i3,' k=',i3,' zuco2=',/,(6e12.4))")
!    |    lat,k,zuco2(k,:)
!     enddo ! k=1,nzpm9m5
!
      nlons = lon1-lon0+1
      do k=1,6
        co2int(1) = cor150(k)
        co2int(2) = cor360(k)
        co2int(3) = cor540(k)
        co2int(4) = cor720(k)
!
        uref(1) = uco2co(k)*150./360.
        uref(2) = uco2co(k)
        uref(3) = uco2co(k)*540./360.
        uref(4) = uco2co(k)*720./360.
!              
        do i=lon0,lon1
          ii = i-lon0+1
          uc(ii) = zuco2(k,i)
        enddo ! i=lon0,lon1
 
!     SUBROUTINE A18LINV(X,Y,XN,YN,N,IMAX)
!     integer,intent(in) :: n,imax
!     real,intent(out)   :: Y(IMAX)
!     real,intent(in)    :: X(IMAX), XN(N), YN(N)
 
!       write(6,"('recur before a18linv: lat=',i3,' k=',i3,' uc=',/,
!    |    (6e12.4))") lat,k,uc
 
        call a18linv(uc,a,uco2ro,alo,51,nlons)   ! a(nlons) is output
!       write(6,"('recur: lat=',i3,' k=',i3,' a=',/,(6e12.4))")
!    |    lat,k,a
 
        call a18linv(uc,cor,uref,co2int,4,nlons) ! cor(nlons) is output
!       write(6,"('recur: lat=',i3,' k=',i3,' cor=',/,(6e12.4))")
!    |    lat,k,cor
 
        do i=lon0,lon1
          ii = i-lon0+1
          al(k,i) = exp(cor(ii)+a(ii))
        enddo ! i=lon0,lon1
      enddo ! k=1,6
      do k=7,nzpm9m5
        do i=lon0,lon1
          ii = i-lon0+1
          uc(ii) = zuco2(k,i)
        enddo
        call a18linv(uc,a,uco2ro,alo,51,nlons) ! a(nlons) is output
!       write(6,"('recur: lat=',i3,' k=',i3,' a=',/,(6e12.4))")
!    |    lat,k,a

        do i=lon0,lon1
          ii = i-lon0+1
          al(k,i) = exp(a(ii))
        enddo
      enddo ! do k=7,nzpm9m5

!     do k=1,nzpm9m5
!       write(6,"('recur: lat=',i3,' k=',i3,' zuco2(k,:)=',/,
!    |    (6e12.4))") lat,k,zuco2(k,:)
!       write(6,"('recur: lat=',i3,' lon0,1=',2i3,' k=',i3,' al(k,:)=',
!    |    /,(6e12.4))") lat,lon0,lon1,k,al(k,:)
!     enddo ! k=1,nzpm9m5

      end subroutine recur
!-----------------------------------------------------------------------
      SUBROUTINE A18LINV(X,Y,XN,YN,N,IMAX)
      implicit none
C     ****
C     ****     This procedure performs linear interpolation within the
C     ****     table defined by the N points (XN(NN),Y(NN)).
C     ****     Where:
C     ****
C     ****       NN = 1,N,1
C     ****
C     ****       XN(NN) < XN(NN+1) for NN = 1,N-1
C     ****
C     ****     Parameters:
C     ****
C     ****       X(IMAX) = array of IMAX x-values at which linear
C     ****                 interpolation is required
C     ****
C     ****       XN(N) = array of N abscissae at which function values
C     ****               are given
C     ****
C     ****       YN(N) = function values corresponding to abscissae,
C     ****               XN(N)
C     ****
C     ****     Output:
C     ****
C     ****      Y(IMAX)  The IMAX interpolated values are
C     ****                     returned in this array
C     ****
!
! Args:
      integer,intent(in) :: n,imax
      real,intent(out)   :: Y(IMAX)
      real,intent(in)    :: X(IMAX), XN(N), YN(N) 
!
! Local:
      integer :: KK(IMAX),i,nn
C     ****
C     ****     Where:
C     ****       Y(IMAX) is vector output
C     ****
C     ****       KK is work space
C     ****
C     ****
C     ****     Initialize array KK
C     ****
      DO I = 1,IMAX
	KK(I) = 0
      ENDDO
C     ****
C     ****     Locate interval in (XN,YN) in table containing X(I)
C     ****
      DO NN = 1,N-1
	DO I = 1,IMAX
          if ((xn(nn+1)-x(i))*(x(i)-xn(nn)) >= 0.) kk(i) = nn+1
	ENDDO
      ENDDO
C     ****
C     ****     Check for
C     ****
C     ****       X(I) < XN(1),  X(I) > X(N)
C     ****
C     ****       and use linear extrapolation if necessary
C     ****
      DO I = 1,IMAX
        if (xn(1)-x(i) >= 0.) kk(i) = 2
        if (x(i) >= xn(n)) kk(i) = n
      ENDDO
C     ****
C     ****     Perform interpolation prescribed above
C     ****
      DO I = 1,IMAX
	Y(I) = (YN(KK(I)-1)*(XN(KK(I))-X(I)) + YN(KK(I))*
     1               (X(I)-XN(KK(I)-1)))/(XN(KK(I))-XN(KK(I)-1))
      ENDDO
      END subroutine a18linv
!-----------------------------------------------------------------------
      subroutine vicool(ztv,zsn2,zso2,zo,zco2,al,am,den,
     |  zflux,zhco2,zho3,nzpm9m5,lev0,lev1,lon0,lon1,lat)
      use cco2gr_module,only: amat,bmat,ao3,ig
      use params_module,only: nlevp1
      implicit none
!
! Args:
      integer,intent(in) :: nzpm9m5,lev0,lev1,lon0,lon1,lat
!
! Inputs are for zp -9 to zp -5:
      real,dimension(nzpm9m5,lon0:lon1),intent(in) ::
     |  zsn2,zso2,zo,zco2,al,am,den
      real,dimension(lev0:lev1,lon0:lon1),intent(in) :: ztv
!
! Output:
      real,intent(out) ::
     |  zflux(lon0:lon1),
     |  zhco2(lmax,lon0:lon1),
     |  zho3 (lmax,lon0:lon1)
!
! Local:
      real,parameter :: 
     |  a10 = 1.5988,
     |  const  = 2.63187e11,
     |  constb = 9.08795e9
      integer :: i,k,j,js,ks,km,kt
      real,dimension(lev0:lev1,lon0:lon1) :: fu,fo3
      real,dimension(nlevp1,lon0:lon1) :: alam
      real,dimension(lon0:lon1) :: h1,h2,h3,tt,zo2,zn2,y,zz,
     |  aa1,aa2,d1,d2
      real :: rko
!
! Grid levels for height integration   (p.s.h. distance = 0.25*IG)
! FU - the LTE source functions for 15 um CO2 band at
!   p.s.h.=0-16.5,
! FO3 - the same for 9.6 um O3 band at p.s.h. = 0-11, 0.25
! AL - the quantum survival probability for
!   p.s.h. = 12.5-16.5, 0.25
!
      do i=lon0,lon1
        do k=lev0,lev1
          fu (k,i) = exp(-960.217/ztv(k,i))
          fo3(k,i) = exp(-1500.  /ztv(k,i))
        enddo
      enddo ! i=lon0,lon1
!
! Cooling rate in co2 bands (x=10.75-12.5, matrix approach)
! cco2gr.F: real :: amat(43,9),bmat(43,9)
!
      do k=lam0,lam1       ! lamdas 36,43
        ks = k-38          ! -2,5 (zp -10.75 to -9)
!       if (ks < 1) ks = 1 ! use -10 for 3 levs from -10.75 up to -10.25 
        do i=lon0,lon1
          h2(i) = 0.
        enddo ! i=lon0,lon1
        do j=1,9
!
! 7/24/06 btf: ig(9) has been changed in cco2gr.F to work for tiegcm here
! (assuming dz=0.25)
          js = ks+ig(j)

!         if (js > 0) then
          if (js > 0 .and. ks > 0) then
            do i=lon0,lon1
              h2(i) = h2(i)+(amat(k,j)+bmat(k,j)*fu(ks,i))*fu(js,i)
            enddo ! i=lon0,lon1

!           write(6,"('vicool: lat=',i3,' k=',i3,' ks=',i4,' zp(ks)=',
!    |        f6.2,' j=',i4,' js=',i4,' zp(js)=',f6.2,' fu(ks,lon0)=',
!    |        e12.4,' fu(js,lon0)=',e12.4,' amat,bmat(k,j)=',2e12.4,
!    |        ' h2(lon0)=',e12.4)")
!    |        lat,k,ks,zpk(ks,zpint(1),dz),j,js,zpk(js,zpint(1),dz),
!    |        fu(ks,lon0),fu(js,lon0),amat(k,j),bmat(k,j),h2(lon0)

          endif
        enddo ! j=1,9
!
! Define output for lam 36,43 (lam0,lam1)
        if (ks > 0) then
          do i=lon0,lon1 
            zhco2(k,i) = h2(i)
            zho3(k,i) = 0.
          enddo ! i=lon0,lon1 
        endif
      enddo ! k=lam0,lam1 ! (36,43)
!
! Determine alam:
      do k=1,nzpm9m5    ! zp -9 to +5 for arrays zsn2,zso2,zo, etc
        do i=lon0,lon1
!
! co2-o2 and co2-n2 v-t constants.
!
          kt = k+4  ! ztv is full column from lbc, so zp -9 is k==5 (assumes dz=0.25)
          tt(i) = ztv(kt,i)
          rko = setrko(tt(i)) ! function setrko is in this module
          y(i) = tt(i)**(-1./3.)
          zn2(i) = 5.5e-17*sqrt(tt(i))+6.7e-10*exp(-83.8*y(i))
          zo2(i) = 1.e-15*exp(23.37-230.9*y(i)+564.*y(i)*y(i))
          zz(i) = (zsn2(k,i)*zn2(i)+zso2(k,i)*zo2(i)+zo(k,i)*rko)*
     |      den(k,i)
          alam(k,i) = a10/(a10+zz(i))
        enddo ! i=lon0,lon1
      enddo ! k=1,nzpm9m5
!
      do i=lon0,lon1
        h1(i) = zhco2(lam1,i)/(zco2(1,i)*(1.-alam(1,i))*constb)
      enddo ! i=lon0,lon1
      do k=2,nzpm9m5 
        km = k-1
        do i=lon0,lon1
          aa1(i)= 1.-alam(km,i)*(1.-.25*al(k,i)-.75*al(km,i))
          aa2(i)= 1.-alam(k ,i)*(1.-.75*al(k,i)-.25*al(km,i))
          d1(i) = -.25*(   al(k,i)+3.*al(km,i))
          d2(i) =  .25*(3.*al(k,i)+   al(km,i))
          h2(i) = (aa1(i)*h1(i)-d1(i)*fu(kzpm9+k-2,i)-
     |                          d2(i)*fu(kzpm9+k-1,i))/aa2(i)
!
! Define output for lam 44,59
          zhco2(lam1+k-1,i) = h2(i)*zco2(k,i)*(1.-alam(k,i))/am(k,i)*
     |      const

!         write(6,"('vicool: i=',i4,' k=',i4,' lam1+k-1=',i4,
!    |      ' zhco2(lam1+k-1,i)=',e12.4,' zp(k)=',f6.2)") i,k,
!    |      lam1+k-1,zhco2(lam1+k-1,i),(kzpm9+k-1)/4.-10.

          zho3(lam1+k-1,i) = 0.
          h1(i) = h2(i)
        enddo ! i=lon0,lon1

      enddo ! k=2,nzpm9m5 
!
! zflux is in longitude only:
      do i=lon0,lon1
        zflux(i) = h2(i)+fu(kzpm5,i)
      enddo ! i=lon0,lon1

!     write(6,"(/,'vicool: lat=',i3,' zflux=',/,(6e12.4))") lat,zflux

!     do k=lam0,lmax
!       write(6,"('vicool: lat=',i3,' k=',i3,' lam0=',i3,' lmax=',i3,
!    |    ' zhco2(k,:)=',/,(6e12.4))") lat,k,lam0,lmax,zhco2(k,:)
!     enddo ! k=1,lmax

      end subroutine vicool
!-----------------------------------------------------------------------
      real function setrko(tn)
      implicit none
!
! Set temperature dependent O-CO2 rate (as per Vladimir Ogibalov, 11/8/02).
! V-T constant for co2(0110) quenching by O
! (measurements by Khvorostovskaya et al. 2002)
! This should be inlined.
!
      real,intent(in) :: tn
!
      if (tn < 260.) then
        setrko = 1.56e-12
      elseif (tn >= 260..and.tn <= 300.) then
        setrko = (2.6-0.004*tn)*1.0e-12
      else
        setrko = 1.4e-12
      endif
      end function setrko
!-----------------------------------------------------------------------
      integer function kzp(zp,zplb,dz)
      implicit none
!
! Given a log pressure zp, log pressure at lower boundary (k==1),
! and delta zp (dz), return k-index corresponding to zp.
!
      real,intent(in) :: zp,zplb,dz
      kzp = int((zp-zplb)/dz)+1
      end function kzp
!-----------------------------------------------------------------------
      real function zpk(k,zplb,dz)
      implicit none
!
! Given a column index k, log pressure at lower boundary (k==1),
! and delta zp (dz), return zp corresponding to index k
! (opposite of function kzp)
!
      integer,intent(in) :: k
      real,intent(in) :: zplb,dz
      zpk = float(k-1)*dz+zplb
      end function zpk
!-----------------------------------------------------------------------
      end module radcool_module
