#include "dims.h"
      SUBROUTINE NEWTO3P
      use cons_module,only: len1,len2,expz,kmax,rmassinv,avo
      implicit none
C     ****
C     ****     E.C.RIDLEY 7/9/87
C     ****
C     ****     ADDS IMPLICIT AND EXPLICIT CONTRIBUTIONS TO COOLING TO
C     ****     NWTI AND NWTE RESPECTIVELY.  S11 AND S12 ARE WORK SPACE.
C     ****
C     ****     COMMON BLOCKS
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
!
! Local:
      real :: AN(3),BN(3)
      integer :: i,k,ntk,nps2k,ncpk,nwtik,nwtek
C     ****
C     ****     CHANGE O(3P) COOLING FACTOR FROM 0.5 TO 1.0 OF THE BATES
C     ****       EXPRESSION
C     ****
      DATA AN/0.835E-18,0.6,0.2/,BN/228.,228.,325./
!
! Need to rewrite base case with below conditionals on xfac before
! comparing with cons_mod.
!
#if (NLEV==28)
! NLEV==28 for -7 to +7 by 0.50
!     real :: xfac(zkmxp)
!     DATA XFAC/3*.01, .05, .1, .2, .4, .55, .7, .75, 19*.8/
      real :: xfac(zkmxp) = 
     |(/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)
! NLEV==56 for -7 to +7 by 0.25
      real :: xfac(zkmxp) = 
     |(/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/)
#endif
      NTK=NJ+NT
      NPS2K=NJ+NPS2
      NCPK=NCP
      DO 1 K=1,KMAX
      DO 1 I=1,LEN1
C     ****     S11=AN(1)*XFAC*N0*PSI2/M2*EXP(-BN(1)/T)
      S11(I,K)=AN(1)*.5*(XFAC(K)+XFAC(K+1))*avo*rmassinv(2)
    1 CONTINUE
      NWTIK=NWTI
      DO 2 I=1,LEN2
      S11(I,1)=S11(I,1)*F(I,NPS2K)*EXP(-BN(1)/F(I,NTK))
C     ****     S12=1.+AN(2)*EXP(-BN(2)/T)+AN(3)*EXP(-BN(3)/T)
      S12(I,1)=1.+AN(2)*EXP(-BN(2)/F(I,NTK))+AN(3)*EXP(-BN(3)/F(I,NTK))
C     ****     S11=S11/S12= F(T(N))
      S11(I,1)=S11(I,1)/S12(I,1)
C     ****     S12=DF/DT=S11/S12/T**2*(BN(1)+(BN(1)-BN(3))*AN(3)*
C     ****          EXP(-BN(3)/T))
      S12(I,1)=S11(I,1)/S12(I,1)/F(I,NTK)**2*(BN(1)+(BN(1)-BN(3))*AN(3)*
     1EXP(-BN(3)/F(I,NTK)))
C     ****     ADD IMPLICIT CONTRIBUTIONS TO NWTI
      F(I,NWTIK)=F(I,NWTIK)+S12(I,1)/(.5*(F(I,NCPK)+F(I,NCPK+1)))
C     ****     ADD EXPLICIT CONTRIBUTIONS TO S10
      S11(I,1)=S11(I,1)-F(I,NTK)*S12(I,1)
    2 CONTINUE
      NWTEK=NWTE-1
      DO 3 K=1,KMAX
      NWTEK=NWTEK+1
      DO 3 I=1,LEN1
      F(I,NWTEK)=F(I,NWTEK)+S11(I,K)*expz(K)
    3 CONTINUE

!     call addfsech('COOL_IMP' ,' ',' ',f(1,nwti),zimxp,zkmxp,zkmxp,j)
!     call addfsech('COOL_EXP' ,' ',' ',f(1,nwte),zimxp,zkmxp,zkmxp,j)

      RETURN
      END
C