#include "dims.h"
      SUBROUTINE CPKTKM
      use cons_module,only: len1,len3,kmax,kmaxp1,rmass,gask,
     |  rmassinv_o2,rmassinv_o,rmassinv_n2
      implicit none
C     ****     CALCULATES CP, KT AND KM
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
!
! Local:
      integer :: i,k 
      integer :: nps1k,nps2k,ncpk,nktk,nkmk,ntk
      real :: fmin,fmax
!
      NPS1K=NJ+NPS
      NPS2K=NPS1K+KMAXP1
      NCPK=NCP
      NKTK=NKT
      NKMK=NKM
      DO 1 I=1,LEN3
C     ****     S1=MBAR
      S1(I,1)=1./(F(I,NPS1K)*rmassinv_o2+F(I,NPS2K)*rmassinv_o+
     1(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2)
C     ****     S2=PO2
      S2(I,1)=S1(I,1)*F(I,NPS1K)*rmassinv_o2
C     ****     S3=PO
      S3(I,1)=S1(I,1)*F(I,NPS2K)*rmassinv_o
C     ****     S4=PN2
      S4(I,1)=1.-S2(I,1)-S3(I,1)
!     S4(I,1) = merge(S4(I,1),1.E-6,S4(I,1)-1.E-6>=0.)
      if (s4(i,1) < 1.e-6) s4(i,1) = 1.e-6
C     ****     CP=R/2.*(PO2*7./32.+PN2*7./28.+PO*5./16.)
      F(I,NCPK)=gask*.5*(S2(I,1)*7./32.+S4(I,1)*7./28.+S3(I,1)*5./16.)
C     ****     KM=(PO2*4.03+PN2*3.42+PO*3.9)
      F(I,NKMK)=S2(I,1)*4.03+S4(I,1)*3.42+S3(I,1)*3.9
C     ****     KT=(PO2+PN2)*56.+PO*75.9
      F(I,NKTK)=(S2(I,1)+S4(I,1))*56.+S3(I,1)*75.9
    1 CONTINUE
      NTK=NJ+NT-1
      NKTK=NKTK-1
      NKMK=NKMK-1
      DO 2 K=1,KMAX
      NKTK=NKTK+1
      NKMK=NKMK+1
      NTK=NTK+1
      DO 2 I=1,LEN1
!     F(I,NTK) = merge(F(I,NTK),100.,F(I,NTK)-100.>=0.)
      if (f(i,ntk) < 100.) f(i,ntk) = 100.
      S5(I,K) = F(I,NTK)**0.69
      F(I,NKMK) = F(I,NKMK)*S5(I,K)*1.E-6
      F(I,NKTK) = F(I,NKTK)*S5(I,K)
    2 CONTINUE
      NKMK=NKMK+1
      NKTK=NKTK+1
      DO 3 I=1,LEN1
      F(I,NKMK)=1.E-6*F(I,NKMK)*F(I,NTK)**0.69
      F(I,NKTK)=F(I,NKTK)*F(I,NTK)**0.69
    3 CONTINUE
      NKMK=NKMK+1
      NKTK=NKTK+1
      NCPK=NCPK+KMAXP1
      DO 4 K=1,KMAX
      NKMK=NKMK-1
      NKTK=NKTK-1
      NCPK=NCPK-1
      DO 4 I=1,LEN1
      F(I,NCPK)=.5*(F(I,NCPK)+F(I,NCPK-1))
      F(I,NKMK)=.5*(F(I,NKMK)+F(I,NKMK-1))
      F(I,NKTK)=.5*(F(I,NKTK)+F(I,NKTK-1))
    4 CONTINUE
      DO 5 I=1,LEN1
      F(I,NCPK-1)=2.*F(I,NCPK-1)-F(I,NCPK)
      F(I,NKMK-1)=2.*F(I,NKMK-1)-F(I,NKMK)
      F(I,NKTK-1)=2.*F(I,NKTK-1)-F(I,NKTK)
    5 CONTINUE
!     call fminmax(f(1,ncp),zimxp*zkmxp,fmin,fmax)
!     write(6,"('cpktkm j=',i2,' f(1,ncp) min,max=',2e12.4)") 
!    |  j,fmin,fmax
!     call fminmax(f(1,nkt),zimxp*zkmxp,fmin,fmax)
!     write(6,"('cpktkm j=',i2,' f(1,nkt) min,max=',2e12.4)") 
!    |  j,fmin,fmax
!     call fminmax(f(1,nkm),zimxp*zkmxp,fmin,fmax)
!     write(6,"('cpktkm j=',i2,' f(1,nkm) min,max=',2e12.4)") 
!    |  j,fmin,fmax
      RETURN
      END