#include "dims.h"
      subroutine radcool
      use cons_module,only: len1,len3,kmax,kmaxp1,t0,expz,p0,boltz,
     |  avo,expzmid,gask,grav,rmassinv_o2,rmassinv_o,rmassinv_n2
!
! Radiative cooling, generalized for two vertical resolutions,
! dz=0.5 and dz=0.25. See also radcool.h, vicool.F, and recur.F.
!
      implicit none
!
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
#include "lowbnd.h"
#include "mwt.h"
#include "radcool.h"
!
! Local:
      integer :: kstep, kzlb, kzpm9, kzpm5, kzpm1925
      real :: TT(zimxp),Y(zimxp),ZN2(zimxp),ZO2(zimxp),ZZ(zimxp),
     |  ALAM(zimxp,ZKMXP)
      real :: TVLB(ntndown+1),O3LB(ntndown+1)
      real :: TVV(ZIMXP,ZKMXP),UNO2(ZIMXP,ZKMXP),UNO(ZIMXP,ZKMXP),
     | UNO3(ZIMXP,ZKMXP),UNN2(ZIMXP,ZKMXP),UNCO2(ZIMXP,ZKMXP),
     | COOL(ZIMXP,ZKMXP),AL(ZIMXP),COOLDD(ZIMXP,ZKMXP),
     | CPPF(ZIMXP),UNAM(ZIMXP,ZKMXP),UNDEN(ZIMXP,ZKMXP),
     | COLCO2(ZIMXP,ZKMXP)
      real :: UTOP(zimxp),ur,a10,factor,const
      integer :: i,k,kk
      integer :: ntk,nps1k,nps2k,npo1k,npo3k,npco2k,nmsk,ncpk,nwtik,
     |  nwtek,npnok
!
! (RKO moved to radcool.h)
      UR=8.3144E+7
      A10=1.5988 
      CONST=2.63187E11
C     ****
C     **** DATA FOR TEMPERATURE AND OZONE MIXING RATIO(PPM) FROM
C     **** THE GROUND TO X=4.5 (19 VALUES AT 0.25 RESOLUTION)
! tvlb not used (use tndown(ntndown))
C     ****
      TVLB=(/288.7,277.5,266.2,254.5,242.7,231.4,221.5,213.9,209.3,
     1  207.9,208.9,211.2,213.8,216.3,218.6,220.7,222.9,225.7,229.5/)
      O3LB=(/0.028,0.034,0.061,0.106,0.166,0.239,0.321,0.428,0.523,
     1  0.661,1.051,1.606,2.335,3.223,4.201,5.190,6.113,6.925,7.646/)
!
      kzpm9 = int((zpm9-zsb)/dz)+1     ! col index at zp -9 from model lb
      kzpm5 = int((zpm5-zsb)/dz)+1     ! col index at zp -5 from model lb
      kzlb  = int((zsb-zpsrf)/dzcool)+1 ! col index to model lb from surface
      kstep = int((dz+.01)/dzcool)     ! 2 for dz 0.5, 1 for dz 0.25
      kzpm1925 = int((zpm1925-zpsrf)/dzcool)+1 ! k at zp -19.25 from surface
!
!     if (j==1)
!    |  write(6,"('radcool: j=',i2,' dz=',f6.2,' kstep=',i2,' kzpm9=',
!    |    i3,' kzpm5=',i3,' kzlb=',i3,' kzpm1925=',i3)") 
!    |    j,dz,kstep,kzpm9,kzpm5,kzlb,kzpm1925
!
      NTK = NJ+NT-1
      NPS1K = NJ+NPS-1
      NPS2K = NJ+NPS2-1
      NPO1K = NJ+NPO1-1
      NPO3K = NJ+NPO3-1
      NPCO2K = NJ+NPCO2-1
      NMSK = NJ+NMS-1
      DO K = 1,KMAX
        NTK = NTK+1
        NPS1K = NPS1K+1
        NPS2K = NPS2K+1
        NPO1K = NPO1K+1
        NPO3K = NPO3K+1
        NPCO2K = NPCO2K+1
        NMSK = NMSK+1
        DO I=1,LEN1
C     ****
C     ****     S3 = N*MBAR  (K+1/2)
C     ****
          S3(I,K)= expz(K)*p0/(boltz*F(I,NTK)*(F(I,NPS1K)*rmassinv_o2+
     1             F(I,NPS2K)*rmassinv_o+(1.-F(I,NPS1K)-F(I,NPS2K))
     2             *rmassinv_n2))
          TVV(I,K) = F(I,NTK)
          UNO2(I,K) = F(I,NPS1K)*F(I,NMSK)*rmassinv_o2
          UNO(I,K) = F(I,NPO1K)*F(I,NMSK)*rmassinv_o
          UNO3(I,K) = F(I,NPO3K)*F(I,NMSK)/RMO3
          UNN2(I,K) = (1.-F(I,NPS1K)-F(I,NPS2K))*F(I,NMSK)*rmassinv_n2
          UNCO2(I,K) = F(I,NPCO2K)*F(I,NMSK)/RMCO2
          UNAM(I,K) = F(I,NMSK)
          UNDEN(I,K) = S3(I,K)/F(I,NMSK)
        enddo
      enddo
C     ****
C     ****   CALCULATE CO2 VERTICAL COLUMN NUMBER DENSITY FOR
C     ****   RECCURANCE RELATIONS
C     ****
      NPCO2K = NJ+NPCO2
      DO I=1,LEN3
        S7(I,1) = F(I,NPCO2K)
      enddo
C     ****      S7=NO. DENSITY INTEGRAL
C     ****      TOP OF MODEL
C     ****      S7(KMAXP1)=N0*P0*PSI(N)*expz*MBAR/(RMASS(N)**2*G)
      FACTOR=avo*p0*expz(KMAX)*expzmid/(RMCO2**2*grav)
      NMSK=NJ+NMS+KMAX
      DO I=1,LEN1
        S7(I,KMAXP1)=FACTOR*.5*(S7(I,KMAX)+S7(I,KMAXP1))*F(I,NMSK)
        COLCO2(I,KMAXP1) = S7(I,KMAXP1)
      enddo
C     ****      INTEGRATE DOWN TO LEVEL 1
      do kk=1,KMAX
        K=KMAXP1-KK
        FACTOR=avo*p0*expz(K)/(RMCO2*grav)*dz
        DO I=1,LEN1
          S7(I,K)=S7(I,K+1)+FACTOR*S7(I,K)
          COLCO2(I,K) = S7(I,K)
        enddo
      enddo
!
! ntndown = 18 (params.h): number of column grid points in tndown
! from ground up to one level below the model bottom boundary.
! This corresponds to zp -21.5 to -17.25 by dzcool (0.25).
! Tndown(ntndown) is from the Solomon-Garcia model 
! (see solgar_mod.F, dlowbnd.h and lowbnd.h)
!
      do kk=1,ntndown               ! zp -21.5 to -17.25 by 0.25 (dzcool)
        DO I=1,LEN1
          O3(I,KK)=O3LB(KK)
          tv(i,kk) = tndown(j,1)    ! constant in vertical
        enddo
      enddo
      do kk=kzpm1925,ntndown        ! zp -19.25 to -17.25
        DO I=1,LEN1
          TV(I,KK)=tndown(J,KK)     ! rewrite tv with tndown from zp -19.25
        enddo
      enddo
      DO K=1,kzpm5                  ! zp -17 (model lb) to zp -5
        kk = kzlb+(k-1)*kstep       ! 19-67,2 for dz0.5, 19-67,1 for dz0.25
        DO I=1,LEN1
          TV(I,KK)=TVV(I,K)
          O3(I,KK)=UNO3(I,K)
        enddo
!
! Interpolate to midpoints for dz=0.50 model:
        if (kstep > 1 .and. k > 1) then ! true for dz0.5 (kstep==2)
          do i=1,len1
            tv(i,kk-1) = ( tvv(i,k) +  tvv(i,k-1))*0.5
            o3(i,kk-1) = (uno3(i,k) + uno3(i,k-1))*0.5
          enddo
        endif
      enddo
!
      DO K=kzpm9,kzpm5    ! 17-25 for dz=0.25 model, 33-49 for dz=0.5 model
        kk = k-kzpm9+1    ! 1-17 for dz=0.25
        if (kstep > 1) 
     |    kk = kk*kstep-1 ! 1,3,5...17 for dz=0.5 model
        DO I=1,LEN1
          SN2(I,KK)=UNN2(I,K)
          SO2(I,KK)=UNO2(I,K)
          O(I,KK)=UNO(I,K)
          CO2(I,KK)=UNCO2(I,K)
          UCO2(I,KK) = COLCO2(I,K)
          AM(I,KK)=UNAM(I,K)
          DEN(I,KK)=UNDEN(I,K)
        enddo
!
! Interpolate to midpoints for dz=0.50 model:
        if (kstep > 1 .and. k > kzpm9) then ! true for dz=0.5 model
          do i=1,len1
            sn2 (i,kk-1) = (unn2  (i,k) + unn2  (i,k-1)) *0.5
            so2 (i,kk-1) = (uno2  (i,k) + uno2  (i,k-1)) *0.5
            o   (i,kk-1) = (uno   (i,k) + uno   (i,k-1)) *0.5
            co2 (i,kk-1) = (unco2 (i,k) + unco2 (i,k-1)) *0.5
            uco2(i,kk-1) = (colco2(i,k) + colco2(i,k-1)) *0.5
            am  (i,kk-1) = (unam  (i,k) + unam  (i,k-1)) *0.5
            den (i,kk-1) = (unden (i,k) + unden (i,k-1)) *0.5
          enddo
        endif
      enddo
!
      DO K=1,nzpsrf_zpm5
        XPZ(K)=0.+(K-1)*0.25 ! not used
      enddo
      CALL RECUR
      CALL VICOOL   ! returns hco2 and ho3
      DO K=1,kzpm5  ! zp -17 to -5: 1->25 for dz=0.5, 1->49 for dz=0.25
        kk = kzpm1925+k            ! 11,12,13...59 for dz=0.25 model
        if (kstep > 1) 
     |    kk = kzpm1925+kstep*k-1  ! 11,13,15...59 for dz=0.5 model
        DO I=1,LEN1
          COOL(I,K)=-HCO2(I,KK)-HO3(I,KK)
          S12(I,K) = -HCO2(I,KK)
          S13(I,K) = -HO3(I,KK)
        enddo
      enddo
C     ****
C     ****   CO2 COOL-TO-SPACE APPROXIMATION
C     ****
      DO K = kzpm5+1,KMAX
        DO I=1,LEN1
          TT(I) = TVV(I,K)
          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))
C
C     ****  COLLISIONAL DEACTIVATION RATE:
! (RKO is in radcool.h)
C
          ZZ(I) = (UNN2(I,K)*ZN2(I)+UNO2(I,K)*ZO2(I)+UNO(I,K)*RKO)
     1            *UNDEN(I,K)
          ALAM(I,K) = A10/(A10+ZZ(I))
          COOL(I,K) = -CONST/UNAM(I,K)*UNCO2(I,K)*(1.-ALAM(I,K))*
     1                (FLUX(I)-EXP(-960.217/TT(I)))
          S12(I,K) = COOL(I,K)
          S13(I,K) = 0.
        enddo
      enddo
      DO K=1,KMAX
        DO I=1,LEN1
          CPPF(I)=86400./(0.5*UR*(7./32.*UNO2(I,K)+7./28.*UNN2(I,K)
     1            +5./16.*UNO(I,K)))
          COOLDD(I,K) = COOL(I,K)*CPPF(I)
        enddo
      enddo
C     ****
C     ****   NO COOL TO SPACE APPROXIMATION(SAME AS BEFORE)
C     ****
      NMSK=NJ+NMS
      NCPK=NCP-1
      NTK=NJ+NT-1
      NPS1K=NJ+NPS-1
      NPS2K=NJ+NPS2-1
      NPO1K=NJ+NPO1-1
      NMSK=NJ+NMS-1
      NPCO2K=NJ+NPCO2-1
      NWTIK=NWTI-1
      NWTEK=NWTE-1
      NPNOK=NJ+NPNO-1
      DO K=1,KMAX
        NCPK=NCPK+1
        NTK=NTK+1
        NPS1K=NPS1K+1
        NPS2K=NPS2K+1
        NPO1K=NPO1K+1
        NPCO2K=NPCO2K+1
        NMSK=NMSK+1
        NWTIK=NWTIK+1
        NWTEK=NWTEK+1
        NPNOK=NPNOK+1
        DO I=1,LEN1
C         ****
C         ****     S10 = N(NO)/RHO
C         ****
          S10(I,K)=avo*F(I,NPNOK)/RMNO
C         ****
C         ****     S11 = A(NO)
C         ****
          S11(I,K)=.5*(F(I,NMSK)+F(I,NMSK+1))*avo*p0*expz(K)/
C    1      (gask*(F(I,NTK)+.5*(T0(K)+T0(K+1))))*(6.5E-11*F(I,NPO1K)/
     1      (gask*(F(I,NTK)+.5*(T0(K)+T0(K+1))))*(2.7E-11*F(I,NPO1K)*
     2      rmassinv_o+2.4E-14*F(I,NPS1K)*rmassinv_o2)
C         ****
C         ****     FORM NO COOLING CONTRIBUTION IN S11
C         ****
          S11(I,K)=4.956E-12*S10(I,K)*(S11(I,K)/(S11(I,K)+13.3))*
     1      EXP(-2700./F(I,NTK))
C         ****
C         ****     FORM TOTAL COOLING TERMS IN NWTI AND NWTE
C         ****
          F(I,NWTIK)=(2700.*S11(I,K)+960.*S12(I,K)+1500.*S13(I,K))/
     1      (.5*(F(I,NCPK)+F(I,NCPK+1))*F(I,NTK)**2)
          F(I,NWTEK)=((1.-2700./F(I,NTK))*S11(I,K)+(1.-960./F(I,NTK))*
     1      S12(I,K)+(1.-1500./F(I,NTK))*S13(I,K))*expz(K)
        enddo
      enddo
!
!     real :: TVV(ZIMXP,ZKMXP),UNO2(ZIMXP,ZKMXP),UNO(ZIMXP,ZKMXP),
!    | UNO3(ZIMXP,ZKMXP),UNN2(ZIMXP,ZKMXP),UNCO2(ZIMXP,ZKMXP),
!    | COOL(ZIMXP,ZKMXP),AL(ZIMXP),COOLDD(ZIMXP,ZKMXP),
!    | CPPF(ZIMXP),UNAM(ZIMXP,ZKMXP),UNDEN(ZIMXP,ZKMXP),
!    | COLCO2(ZIMXP,ZKMXP)
!
!     call addfsech('COOL'  ,' ',' ',cool  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('ALAM'  ,' ',' ',alam  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('TVV'   ,' ',' ',tvv   ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNO2'  ,' ',' ',uno2  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNO'   ,' ',' ',uno   ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNO3'  ,' ',' ',uno3  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNN2'  ,' ',' ',unn2  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNCO2' ,' ',' ',unco2 ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNAM'  ,' ',' ',unam  ,zimxp,zkmxp,zkmx,j)
!     call addfsech('UNDEN' ,' ',' ',unden ,zimxp,zkmxp,zkmx,j)
!     call addfsech('COLCO2',' ',' ',colco2,zimxp,zkmxp,zkmx,j)
!     call addfsech('RADC_S10',' ',' ',s10 ,zimxp,zkmxp,zkmx,j)
!     call addfsech('RADC_S11',' ',' ',s11 ,zimxp,zkmxp,zkmx,j)
!     call addfsech('RADC_S13',' ',' ',s13 ,zimxp,zkmxp,zkmx,j)
!     call addfsech('NWTI',' ',' ',f(1,nwti),zimxp,zkmxp,zkmx,j)
!     call addfsech('NWTE',' ',' ',f(1,nwte),zimxp,zkmxp,zkmx,j)
      RETURN
      END