!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE KFINTERFC                  ######
!######                                                      ######
!######                     Developed by                     ######
!######          Coastal Meteorology Research Program        ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE kfinterfc(nx,ny,nz,u,v,w,pprt,ptprt,qv,pbar,ptbar,zp,        & 1,1
           ptcumsrc,qcumsrc,rainc,prcrate,nca,raincv,                   &
           ptop,psfc,psb,                                               &
           sigma,hfsig,ub,vb,w0,tb,qvb,tem1,tem2)

!
!-----------------------------------------------------------------------
!   PURPOSE:
!
!    To calculate the necessary inputs for the Kain-Fritsch
!    cumulus parameterization scheme. Note that the inputs
!    to the k-f scheme are arranged in such a way that the
!    k index is counted from top down, i.e., k=1 is top and
!    k=kx is the bottom.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Zonghui Huo
!  08/01/97
!
!  MODIFICATION HISTORY:
!
!  1/6/98 (Z. Huo and J. Kain)
!  Fixed a bug in the calculation of hfsig. Optimize the code the
!  call to kfpara.f.
!
!  2/3/98 Zonghui Huo
!  A bug fix in loop 563. tem2(i,j,k) = 0.5*(tem1(i,j,k)-tem1(i,j,k-1))
!  should have been       tem2(i,j,k) = 0.5*(tem1(i,j,k)+tem1(i,j,k-1))
!
!  2/9/1998 (M. Xue)
!  Made a few additional optimization.
!
!  4/15/1998 (Donghai Wang)
!  Added the source terms for qc,qr,qi and qs equations.
!  Added the running average vertical velocity instead of the 2-level
!  time average.
!  Passed two options mphyopt and kffbfct from arps.input to K-F scheme.
!
!  4/18/2002 (Zuwen He)
!  Passed sub-saturation option kfsubsattrig from arps.input to K-F scheme.
!
!----------------------------------------------------------------------

  IMPLICIT NONE

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        time-average vertical velocity (m/s)
!    pprt     Perturbation pressure (Pascal)
!    ptprt    Perturbation potential temperature (K)
!    qv       Water vapor specific humidity (kg/kg)
!    pbar     Base state pressure (Pascal)
!    ptbar    Base state potential temperature (K)
!    zp       The physical height coordinate defined at w-point
!
!  OUTPUT:
!
!    ptcumsrc Potential temperature source term.
!    qcumsrc Water specific humidity source term
!    rainc    Accumulated precipitation by cumulus convection (mm)
!    prcrate  precipitatioon rate (mm/s)
!    raincv   K-F convective rainfall (cm)
!    nca      K-F counter for CAPE release
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  INTEGER :: nx,ny,nz

  REAL :: u      (nx,ny,nz) ! Total u-velocity (m/s)
  REAL :: v      (nx,ny,nz) ! Total v-velocity (m/s)
  REAL :: w      (nx,ny,nz) ! Total w-velocity (m/s)

  REAL :: pprt   (nx,ny,nz) ! Perturbation pressure (Pascal)
  REAL :: ptprt  (nx,ny,nz) ! Perturbation potential temperature (K)
  REAL :: qv     (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
  REAL :: pbar   (nx,ny,nz) ! Base state pressure (Pascal)
  REAL :: ptbar  (nx,ny,nz) ! Base state potential temperature (K)

  REAL :: zp     (nx,ny,nz) ! The physical height coordinate defined
                            ! at w-point
!
! Output array
!

  REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt equation.
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: rainc(nx,ny)       ! Accumulated precipitation by convec(mm)
  REAL :: prcrate(nx,ny)     ! precipitation rate (kg/(m**2*s))
  REAL :: raincv  (nx,ny)    ! K-F output (cm)
  INTEGER :: nca     (nx,ny)    ! counter for CAPE release

!
! Local array
!

  REAL :: ptop (nx,ny)    ! Pressure at the model top (cb)
  REAL :: psfc (nx,ny)    ! Surface pressure (cb=1000Pa)
  REAL :: psb  (nx,ny)    ! Psfc - Ptop, used by mm5
  REAL :: sigma(nx,ny,nz) ! sigma levels where w is defined
  REAL :: hfsig(nx,ny,nz) ! half sigma levels where all variables defined
  REAL :: ub   (nx,ny,nz) ! u at the mass points (times pstr)
  REAL :: vb   (nx,ny,nz) ! v at the mass points
  REAL :: w0   (nx,ny,nz) ! w at the mass points
  REAL :: tb   (nx,ny,nz) ! Temperature (K) at the mass points
  REAL :: qvb  (nx,ny,nz) ! Water vapor specific humidity (kg/kg) at
                          ! mass points of the physical domain
!
! More local array related to the optimization to call the kfpatra.f
!
  INTEGER :: ncuyes,k300,kmax,kmin,nzmax,nxmax
  REAL :: rovg,psfck,p300,ttmp,qtmp,es,qs,semmx,semmin
  REAL :: cp,xlv,g,aliq,bliq,cliq,dliq
  PARAMETER(nzmax=100,nxmax=500)
  PARAMETER(cp=1005.7,xlv=2.5E6,g=9.8,rovg=29.29)
  PARAMETER(aliq=613.3,bliq=17.502,cliq=4780.8,dliq=32.19)
  INTEGER :: icuyes(nxmax)
  INTEGER :: lsb   (nxmax)
  REAL :: sem   (nzmax)
  REAL :: sems  (nzmax)
  REAL :: prs   (nzmax)
  REAL :: tv0   (nzmax)
  REAL :: z0    (nzmax)

  REAL :: tem1   (nx,ny,nz) ! temporary
  REAL :: tem2   (nx,ny,nz) ! temporary

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ix,jx,kx,nk
  REAL :: delz,tv, dtbig10
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Set the loop index for the physical space mass point since
!  the K-F scheme is only performed in physical space on mass
!  points
!
!-----------------------------------------------------------------------
!

  ix=nx-1
  jx=ny-1
  kx=nz-3
!
!-----------------------------------------------------------------------
!
! Calculate the pressure at model top (ptop) and model bottom (psfc)
! The temprature (K) is writen in tv.
!
!-----------------------------------------------------------------------
!

! Surface pressure

  DO j=1,ny-1
    DO i=1,nx-1
      delz = (zp(i,j,3)-zp(i,j,2))/2.0
      tv = (ptbar(i,j,2)+ptprt(i,j,2))                                  &
          * (1.0E-5*(pbar(i,j,2)+pprt(i,j,2)))**0.286
      psfc(i,j) = (1.0E-3)*(pbar(i,j,2)+pprt(i,j,2))                    &
                  *EXP(9.8*delz/(287.0*tv))

    END DO
  END DO

! Model top pressure

  DO j=1,ny-1
    DO i=1,nx-1
      delz = (zp(i,j,nz-1)-zp(i,j,nz-2))/2.0
      tv = (ptbar(i,j,nz-2)+ptprt(i,j,nz-2))                            &
          * (1.0E-5*(pbar(i,j,nz-2)+pprt(i,j,nz-2)))**0.286
      ptop(i,j) = (1.0E-3)*(pbar(i,j,nz-2)+pprt(i,j,nz-2))              &
                  *EXP(-9.8*delz/(287.0*tv))

    END DO
  END DO

  DO j=1,jx
    DO i=1,ix
      psb(i,j) = psfc(i,j) -  ptop(i,j)
    END DO
  END DO

!
! Calculate the sigma levels corresponding to the Gal-Chen levels
! and the half sigma levels where all other variables are defined
!

  DO k=1,kx
    DO j=1,jx
      DO i=1,ix
        tem1(i,j,k)=((pbar(i,j,k+1)+pprt(i,j,k+1))*(1.0E-3)             &
                  -ptop(i,j))/(psfc(i,j)-ptop(i,j))
      END DO
    END DO
  END DO

  DO j=1,jx
    DO i=1,ix
      tem2(i,j,1) = 1.0
      tem2(i,j,kx+1) = 0.0
    END DO
  END DO

  DO k=2,kx
    DO j=1,jx
      DO i=1,ix
        tem2(i,j,k) = 0.5*(tem1(i,j,k)+tem1(i,j,k-1))
      END DO
    END DO
  END DO

  DO k=1,kx
    nk=kx-k+1
    DO j=1,jx
      DO i=1,ix
        hfsig(i,j,nk) = tem1(i,j,k)
      END DO
    END DO
  END DO

  DO k=1,kx+1
    nk=kx-k+2
    DO j=1,jx
      DO i=1,ix
        sigma(i,j,nk) = tem2(i,j,k)
      END DO
    END DO
  END DO

!
! Get u and v winds on mass points (Arakawa c grid in ARPS)
!

  DO k=1,kx
    nk=kx-k+1
    DO j=1,jx
      DO i=1,ix
        ub(i,j,nk) =  (u(i,j,k+1)+u(i+1,j,k+1))*0.5
        vb(i,j,nk) =  (v(i,j,k+1)+v(i,j+1,k+1))*0.5
      END DO
    END DO
  END DO

!
! Calculate the w on the mass points (vertical half levels)
!

  DO k=1,kx
    nk=kx-k+1
    DO j=1,jx
      DO i=1,ix
        w0(i,j,nk) =  (w(i,j,k+1)+w(i,j,k+2))*0.5
      END DO
    END DO
  END DO

!
! Calculate tb (in K) and put qvb on mass levels, Rd/Cp=2858566
!

  DO k=1,kx
    nk=kx-k+1
    DO j=1,jx
      DO i=1,ix
        tb(i,j,nk) = (ptbar(i,j,k+1)+ptprt(i,j,k+1))*                   &
                     ((pbar(i,j,k+1) +pprt(i,j,k+1))*                   &
                      1.0E-5)**0.2858566
      END DO
    END DO
  END DO

  DO k=1,kx
    nk=kx-k+1
    DO j=1,jx
      DO i=1,ix
        qvb(i,j,nk) =  qv(i,j,k+1)
      END DO
    END DO
  END DO


!-----------------------------------------------------------------------
!
  DO k=nz-3,1,-1
    nk=nz-k
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = qcumsrc(i,j,nk-1,2)
        tem2(i,j,k) = qcumsrc(i,j,nk-1,3)
      END DO
    END DO
  END DO

  DO k=nz-3,1,-1
    DO j=1,ny-1
      DO i=1,nx-1
        qcumsrc(i,j,k,2) = tem1(i,j,k)
        qcumsrc(i,j,k,3) = tem2(i,j,k)
      END DO
    END DO
  END DO

  DO k=nz-3,1,-1
    nk=nz-k
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = qcumsrc(i,j,nk-1,4)
        tem2(i,j,k) = qcumsrc(i,j,nk-1,5)
      END DO
    END DO
  END DO

  DO k=nz-3,1,-1
    DO j=1,ny-1
      DO i=1,nx-1
        qcumsrc(i,j,k,4) = tem1(i,j,k)
        qcumsrc(i,j,k,5) = tem2(i,j,k)
      END DO
    END DO
  END DO


!
!-----------------------------------------------------------------------
!
! Call the kain-Fritsch scheme. It is checked by east-west i slice.
!
! Initialize temperary arrays, which are used to hold dt/dt and dq/dt.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem1(i,j,k) = 0.0
        tem2(i,j,k) = 0.0
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!...make a quick for check for these things that can eliminate
!...grid points from the possibility of convective initiation:
!
!     1.) convection already active at this point
!     2.) point close to grid-domain boundary
!     3.) downward motion at all levels in lowest 300 mb
!     4.) no CAPE
!
!-----------------------------------------------------------------------
!
 !!!! Peter, the variable lists in the multitasking directives below
 !!!! (in SGI compiler format) should be looked at carefully since they
 !!!! were used with earlier versions of the code.  But, the idea that
 !!!! they put into effect, which is one of sending the KF call for individual
 !!!! j-slices to different processors, should be fairly easy to activate.
!c$doacross
!c$& share(ix,jx,kx,w0,nca,psfc,a,psb,ptop,
!c$&       tb,qvb,rovg),
!c$& local(i,j,k,nk,ncuyes,icuyes,psfck,p300,k300,prs,ttmp,qtmp,
!c$&       tv0,z0,es,qs,sem,sems,semmx,semmin,kmax,kmin)
!
!
!-----------------------------------------------------------------------
!
  DO j=1,jx
    ncuyes=0
    DO i = 1,ix
      icuyes(i)=0
      lsb(i) = 0
    END DO
    DO i = 1,ix
!
!...if parameterization already active at this point, go on to next point...
!
      IF (nca(i,j) > 0) CYCLE
      psfck = psfc(i,j)*1.e3
!
!...find highest model level k-value in lowest 300 mb...
!
      p300 = psfck-3.e4
      k300 = 1
      DO k = 2,kx
        nk = kx-k+1
        prs(k)=1.e3*(hfsig(i,j,nk)*psfc(i,j)+ptop(i,j))
        IF(prs(k) < p300) EXIT
        k300 = k
      END DO
!      15     CONTINUE
!
!...vertical velocity must be upward at some level in the
!...lowest 300 mb...
!
      DO k = 1,k300
        IF(w0(i,j,k) > 0.) GO TO 25
      END DO
      CYCLE
      25     CONTINUE
!
!...calculate moist static energy and saturation moist static energy
!...at each vertical level...
!
      DO k = kx,1,-1
        ttmp = tb(i,j,k)
        qtmp = qvb(i,j,k)
        prs(k)=1.e3*(hfsig(i,j,k)*psfc(i,j)+ptop(i,j))
        tv0(k) = ttmp*(1.+.608*qtmp)
        IF(k == kx) THEN
          z0(k) = 0.
        ELSE
          z0(k) = z0(k+1)-rovg*0.5*(tv0(k)+tv0(k+1))*                   &
                ALOG(prs(k)/prs(k+1))
        END IF
        es=aliq*EXP((bliq*ttmp-cliq)/(ttmp-dliq))
        qs=0.622*es/(prs(k)-es)
        sems(k)=cp*ttmp+xlv*qs+g*z0(k)
        sem(k)=cp*ttmp+xlv*qtmp+g*z0(k)
      END DO
!
!...determine level of maximum MSE in lowest 300 mb...
!
      semmx = 0.
      DO nk =1,k300
        k = kx-nk+1
        semmx=AMAX1(semmx,sem(k))
        IF(semmx == sem(k)) kmax=k
      END DO
!
!...determine level of minimum saturated MSE above kmax...
!
      semmin = 1.e6
      DO k = kmax,1,-1
        semmin=AMIN1(semmin,sems(k))
        IF(semmin == sems(k)) kmin=k
      END DO
      IF(semmx > semmin)THEN
        ncuyes=ncuyes+1
        icuyes(ncuyes)=i
        lsb(i) = kx-kmax+1
      END IF
    END DO
    IF(ncuyes > 0)THEN

      CALL kfpara(nx,ny,nz,j,dtbig,dx,mphyopt,kffbfct,                  &
                  kfsubsattrig,                                         & 
                  ptop,psb,sigma,hfsig,ub,vb,w0,tb,qvb,                 &
                  tem1,tem2,qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),          &
                  qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                    &
                  nca,raincv,ncuyes,icuyes,lsb)
    END IF
  END DO

!
!-----------------------------------------------------------------------
!
! Put the K-F results into the forcing terms (note the index change,
! and the lateral boundary values, the k-f temperature rate should be
! converted back to potential temperature rate)) Rd/Cp=2858566
!
!-----------------------------------------------------------------------
!
  DO k=2,nz-2
    nk=nz-k
    DO j=1,ny-1
      DO i=1,nx-1
        ptcumsrc(i,j,k) = ptcumsrc(i,j,k)                               &
                        + tem1(i,j,nk-1)*(1.0E5/(pbar(i,j,k)            &
                        + pprt(i,j,k)))**0.2858566
        qcumsrc(i,j,k,1) = qcumsrc(i,j,k,1) + tem2(i,j,nk-1)
      END DO
    END DO
  END DO

  DO k=2,nz-2
    nk=nz-k
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = qcumsrc(i,j,nk-1,2)
        tem2(i,j,k) = qcumsrc(i,j,nk-1,3)
      END DO
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny-1
      DO i=1,nx-1
        qcumsrc(i,j,k,2) = tem1(i,j,k)
        qcumsrc(i,j,k,3) = tem2(i,j,k)
      END DO
    END DO
  END DO

  DO k=2,nz-2
    nk=nz-k
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = qcumsrc(i,j,nk-1,4)
        tem2(i,j,k) = qcumsrc(i,j,nk-1,5)
      END DO
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny-1
      DO i=1,nx-1
        qcumsrc(i,j,k,4) = tem1(i,j,k)
        qcumsrc(i,j,k,5) = tem2(i,j,k)
      END DO
    END DO
  END DO

  dtbig10 = 10.0/dtbig

  DO j=1,ny-1
    DO i=1,nx-1
      prcrate(i,j) = raincv(i,j)*dtbig10
    END DO
  END DO
  RETURN
END SUBROUTINE kfinterfc