!!>****h* SpatialCell_GPU/makeKs_utility
!! SYNOPSIS:  
!!   
!! CREATOR: 
!!      Hoang-Trong Minh Tuan 
!! AFFILIATION: 
!!      George Mason University 
!! 
!! HISTORY
!!   * Rev: 0.1
!!   * 0.2: Incorporate Leak code utility (Aug,18-2010:Tuan)
!!   * 0.7: add LCC 7-state model (Aug, 12, 2010 - Tuan)
!!   * 1.0: add new algorithm setup_LCC, setup_RyR that work with any number of states (March, 2, 2012 - Tuan)
!!_****
MODULE makeKs_utility
  USE precisioninfo
  USE math_utility
  USE matrix_utility
  USE m_inspar
  USE error_codes
  !USE data_global
  IMPLICIT NONE
#include "macro.h"  
  ! DEBUG = 1: test akLv1, akLv2, akLm1, akLp1
  !         2: test bigK_Lv1
  !         3: test ballbin
#define _DEBUG_ 0

   INTERFACE read_RyR_Markov
     module procedure read_RyR_Markov_2
     module procedure read_RyR_Markov_4
   END INTERFACE
   INTERFACE setup_RyR
     module procedure setup_RyR_2
     module procedure setup_RyR_4 
   END INTERFACE
  
CONTAINS
  !!>****f* 
  !! SYNOPSIS  
  !!    Work for L-type and RyR (Full ECC)
  !!    Matrix dimension: compK(Mbig,0:maxnk)
  !! INPUTS 
  !! OUTPUTS 
  !!   (for L-type)
  !!   compKv1, compKv2, idxKv1, idxKv2
  !!   compKm1, idxKm1
  !!   compKp1, idxKp1
  !!   maxnkv1, maxnkv2
  !!   maxnkm1
  !!   maxnkp1
  !!   (for RyR)
  !!   compKm2, idxKm2
  !!   compKp2, idxKp2
  !!   maxnkm2, maxnkp2
  !!   (common)
  !!   U_Lo, U_Ro        :  vector keep number of open channels at each release-site's state
  !!   irow_L, irow_R    :  vector size
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
!!$  SUBROUTINE makecompKs_2(compKv1, compKv2, compKm1, compKm2, compKp1, &
!!$       compKp2, idxKv1, idxKv2, idxKm1, idxKm2, idxKp1, &
!!$       idxKp2, U_Lo, U_Ro, N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, &
!!$       maxnkm2, maxnkp1, maxnkp2)
!!$
!!$    IMPLICIT NONE
!!$
!!$    !NOTE: these compK matrices have zero-th columns
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv1, compKv2, compKm1, &
!!$         compKm2, compKp1, compKp2
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxKv1, idxKv2, idxKm1, idxKm2, &
!!$         idxKp1, idxKp2
!!$    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro, U_Lo
!!$    INTEGER :: N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, maxnkm2, maxnkp1, maxnkp2
!!$
!!$                                !local
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1  
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv2  
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lm1 
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lp1 
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rm
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rp
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Rm, idxK_Rp
!!$    INTEGER :: maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1, &
!!$         maxnkRm, maxnkRp, Mbig, ii
!!$
!!$                                ! L-TYPE channel
!!$    CALL setupLCC_6state(N_L, &
!!$         compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$         idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$         bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$                                ! RyR channel
!!$    CALL setupRyR_2state(N_R, &
!!$         compK_Rm, compK_Rp, &
!!$         idxK_Rm, idxK_Rp, &
!!$         bigU_R, irow_R, maxnkRm, maxnkRp)
!!$
!!$    Mbig = irow_L*irow_R
!!$    CALL getclustercompKleft(compK_Lv1, idxK_Lv1, &
!!$         irow_L, irow_R, maxnkLv1, compKv1, idxKv1, maxnkv1)
!!$    CALL getclustercompKleft(compK_Lv2, idxK_Lv2, &
!!$         irow_L, irow_R, maxnkLv2, compKv2, idxKv2, maxnkv2)
!!$    CALL getclustercompKleft(compK_Lm1, idxK_Lm1, &
!!$         irow_L, irow_R, maxnkLm1, compKm1, idxKm1, maxnkm1)
!!$    CALL getclustercompKright(compK_Rm, idxK_Rm, &
!!$         irow_L, irow_R, maxnkRm, compKm2, idxKm2, maxnkm2)
!!$    CALL getclustercompKleft(compK_Lp1, idxK_Lp1, &
!!$         irow_L, irow_R, maxnkLp1, compKp1, idxKp1, maxnkp1)
!!$    CALL getclustercompKright(compK_Rp, idxK_Rp, &
!!$         irow_L, irow_R, maxnkRp, compKp2, idxKp2, maxnkp2)
!!$
!!$    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
!!$    ALLOCATE(U_Ro(Mbig))
!!$
!!$    IF ( ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
!!$    ALLOCATE(U_Lo(Mbig))
!!$
!!$
!!$    DO ii = 1, Mbig, irow_R
!!$       U_Ro(ii:ii+irow_R-1) = bigU_R(1:irow_R)
!!$       U_Lo(ii:ii+irow_R-1) = bigU_L((ii-1)/irow_R + 1)
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "U_Lo"
!!$    WRITE(*,*) (U_Lo(ii), ii=1,Mbig)
!!$#endif
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    IF (ALLOCATED(bigU_R)) DEALLOCATE(bigU_R)
!!$
!!$    RETURN
!!$  END SUBROUTINE makecompKs_2

  !!>****f* 
  !! SYNOPSIS  
  !!    Work for L-type and RyR (Full ECC)
  !!    Transpose version (suppose to replace non-transpose version
  !!    to run on GPU more effectively) 
  !!    Matrix dimension: compK(0:maxnk, Mbig)
  !! INPUTS 
  !! OUTPUTS 
  !!   (for L-type)
  !!   compKv1, compKv2, idxKv1, idxKv2
  !!   compKm1, idxKm1
  !!   compKp1, idxKp1
  !!   maxnkv1, maxnkv2
  !!   maxnkm1
  !!   maxnkp1
  !!   (for RyR)
  !!   compKm2, idxKm2
  !!   compKp2, idxKp2
  !!   maxnkm2, maxnkp2
  !!   (common)
  !!   U_Lo, U_Ro        :  vector keep number of open channels at each release-site's state
  !!   irow_L, irow_R    :  vector size
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO
  !!  makecompKs_2
  !!_****
!!$  SUBROUTINE makecompKs_2_T(compKv1, compKv2, compKm1, compKm2, compKp1, &
!!$       compKp2, idxKv1, idxKv2, idxKm1, idxKm2, idxKp1, &
!!$       idxKp2, U_Lo, U_Ro, N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, &
!!$       maxnkm2, maxnkp1, maxnkp2)
!!$
!!$    IMPLICIT NONE
!!$
!!$                                !NOTE: these compK matrices have zero-th columns
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv1, compKv2, compKm1, &
!!$         compKm2, compKp1, compKp2
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxKv1, idxKv2, idxKm1, idxKm2, &
!!$         idxKp1, idxKp2
!!$    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro, U_Lo
!!$    INTEGER :: N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, maxnkm2, maxnkp1, maxnkp2
!!$
!!$                                !local
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1  
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv2  
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lm1 
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lp1 
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rm
!!$    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rp
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Rm, idxK_Rp
!!$    INTEGER :: maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1, &
!!$         maxnkRm, maxnkRp, Mbig, ii
!!$
!!$                                ! L-TYPE channel
!!$    CALL setupLCC_6state_T(N_L, &
!!$         compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$         idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$         bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$                                ! RyR channel
!!$    CALL setupRyR_2state_T(N_R, &
!!$         compK_Rm, compK_Rp, &
!!$         idxK_Rm, idxK_Rp, &
!!$         bigU_R, irow_R, maxnkRm, maxnkRp)
!!$
!!$    Mbig = irow_L*irow_R
!!$    CALL getclustercompKleft_T(compK_Lv1, idxK_Lv1, &
!!$         irow_L, irow_R, maxnkLv1, compKv1, idxKv1, maxnkv1)
!!$    CALL getclustercompKleft_T(compK_Lv2, idxK_Lv2, &
!!$         irow_L, irow_R, maxnkLv2, compKv2, idxKv2, maxnkv2)
!!$    CALL getclustercompKleft_T(compK_Lm1, idxK_Lm1, &
!!$         irow_L, irow_R, maxnkLm1, compKm1, idxKm1, maxnkm1)
!!$    CALL getclustercompKright_T(compK_Rm, idxK_Rm, &
!!$         irow_L, irow_R, maxnkRm, compKm2, idxKm2, maxnkm2)
!!$    CALL getclustercompKleft_T(compK_Lp1, idxK_Lp1, &
!!$         irow_L, irow_R, maxnkLp1, compKp1, idxKp1, maxnkp1)
!!$    CALL getclustercompKright_T(compK_Rp, idxK_Rp, &
!!$         irow_L, irow_R, maxnkRp, compKp2, idxKp2, maxnkp2)
!!$
!!$    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
!!$    ALLOCATE(U_Ro(Mbig))
!!$
!!$    IF ( ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
!!$    ALLOCATE(U_Lo(Mbig))
!!$
!!$
!!$    DO ii = 1, Mbig, irow_R
!!$       U_Ro(ii:ii+irow_R-1) = bigU_R(1:irow_R)
!!$       U_Lo(ii:ii+irow_R-1) = bigU_L((ii-1)/irow_R + 1)
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "U_Lo"
!!$    WRITE(*,*) (U_Lo(ii), ii=1,Mbig)
!!$#endif
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    IF (ALLOCATED(bigU_R)) DEALLOCATE(bigU_R)
!!$
!!$    RETURN
!!$  END SUBROUTINE makecompKs_2_T


  !! Version that read procedure as argument
  ! These procedure do call appropriate LCC and RyR models
  ! This function is expected to work with any LCC/RyR models
  !   as long as LCC has v1,v2,m1, p1
  !              RyR has m2,p2
  SUBROUTINE makecompKs_3_T(compKv1, compKv2, compKm1, compKm2, compKp1, &
       compKp2, idxKv1, idxKv2, idxKm1, idxKm2, idxKp1, &
       idxKp2, U_Lo, U_Ro, N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, &
       maxnkm2, maxnkp1, maxnkp2, setupLCC, setupRyR)

    IMPLICIT NONE

    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv1, compKv2, compKm1, &
         compKm2, compKp1, compKp2
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxKv1, idxKv2, idxKm1, idxKm2, &
         idxKp1, idxKp2
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro, U_Lo
    INTEGER :: N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, maxnkm2, maxnkp1, maxnkp2
    INTERFACE
       SUBROUTINE setupLCC(N_L, &
            compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
            idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
            bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
         USE precisioninfo
         INTEGER, INTENT(in) :: N_L        ! # of LCC 
         REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
              idxK_Lp1
         REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
              compK_Lm1, compK_Lp1
         INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
       END SUBROUTINE setupLCC
       
       SUBROUTINE setupRyR(N_R, &
            compK_Rm, compK_Rp, &
            indexK_Rm, indexK_Rp, &
            bigU_R, irow_R, maxnkRm, maxnkRp, model)
         USE precisioninfo
         INTEGER, INTENT(IN) :: N_R   ! # of RyR channels
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: indexK_Rm, &
              indexK_Rp
         REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Rm, compK_Rp
         REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
         INTEGER, INTENT(OUT) :: irow_R, maxnkRm, maxnkRp
         CHARACTER(*), OPTIONAL:: model
       END SUBROUTINE setupRyR
    END INTERFACE
    !local
    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L
    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1  
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv2  
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lm1 
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lp1 
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
         idxK_Lp1
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rm
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rp
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Rm, idxK_Rp
    INTEGER :: maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1, &
         maxnkRm, maxnkRp, Mbig, ii

    ! L-TYPE channel
    CALL setupLCC(N_L, &
         compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
         idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
         bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)

    ! RyR channel
    CALL setupRyR(N_R, &
         compK_Rm, compK_Rp, &
         idxK_Rm, idxK_Rp, &
         bigU_R, irow_R, maxnkRm, maxnkRp)

    Mbig = irow_L*irow_R
    CALL getclustercompKleft_T(compK_Lv1, idxK_Lv1, &
         irow_L, irow_R, maxnkLv1, compKv1, idxKv1, maxnkv1)
    CALL getclustercompKleft_T(compK_Lv2, idxK_Lv2, &
         irow_L, irow_R, maxnkLv2, compKv2, idxKv2, maxnkv2)
    CALL getclustercompKleft_T(compK_Lm1, idxK_Lm1, &
         irow_L, irow_R, maxnkLm1, compKm1, idxKm1, maxnkm1)
    CALL getclustercompKright_T(compK_Rm, idxK_Rm, &
         irow_L, irow_R, maxnkRm, compKm2, idxKm2, maxnkm2)
    CALL getclustercompKleft_T(compK_Lp1, idxK_Lp1, &
         irow_L, irow_R, maxnkLp1, compKp1, idxKp1, maxnkp1)
    CALL getclustercompKright_T(compK_Rp, idxK_Rp, &
         irow_L, irow_R, maxnkRp, compKp2, idxKp2, maxnkp2)

    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(Mbig))

    IF ( ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
    ALLOCATE(U_Lo(Mbig))


    DO ii = 1, Mbig, irow_R
       U_Ro(ii:ii+irow_R-1) = bigU_R(1:irow_R)
       U_Lo(ii:ii+irow_R-1) = bigU_L((ii-1)/irow_R + 1)
    ENDDO
#if _DEBUG_ == 3
    WRITE(*,*) "U_Lo"
    WRITE(*,*) (U_Lo(ii), ii=1,Mbig)
#endif
    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
    IF (ALLOCATED(bigU_R)) DEALLOCATE(bigU_R)

    RETURN
  END SUBROUTINE makecompKs_3_T


  SUBROUTINE makecompKs_3(compKv1, compKv2, compKm1, compKm2, compKp1, &
       compKp2, idxKv1, idxKv2, idxKm1, idxKm2, idxKp1, &
       idxKp2, U_Lo, U_Ro, N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, &
       maxnkm2, maxnkp1, maxnkp2, setupLCC, setupRyR)

    IMPLICIT NONE

    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv1, compKv2, compKm1, &
         compKm2, compKp1, compKp2
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxKv1, idxKv2, idxKm1, idxKm2, &
         idxKp1, idxKp2
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro, U_Lo
    INTEGER :: N_L, N_R, irow_L, irow_R, maxnkv1, maxnkv2, maxnkm1, maxnkm2, maxnkp1, maxnkp2
    INTERFACE
       SUBROUTINE setupLCC(N_L, &
            compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
            idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
            bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
         USE precisioninfo
         INTEGER, INTENT(in) :: N_L        ! # of LCC 
         REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
              idxK_Lp1
         REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
              compK_Lm1, compK_Lp1
         INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
       END SUBROUTINE setupLCC
       
       SUBROUTINE setupRyR(N_R, &
            compK_Rm, compK_Rp, &
            indexK_Rm, indexK_Rp, &
            bigU_R, irow_R, maxnkRm, maxnkRp, model)
         USE precisioninfo
         INTEGER, INTENT(IN) :: N_R   ! # of RyR channels
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: indexK_Rm, &
              indexK_Rp
         REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Rm, compK_Rp
         REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
         INTEGER, INTENT(OUT) :: irow_R, maxnkRm, maxnkRp
         CHARACTER(*), OPTIONAL:: model
       END SUBROUTINE setupRyR
    END INTERFACE
    !local
    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L
    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_R
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1  
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lv2  
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lm1 
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Lp1 
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
         idxK_Lp1
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rm
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK_Rp
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Rm, idxK_Rp
    INTEGER :: maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1, &
         maxnkRm, maxnkRp, Mbig, ii

    ! L-TYPE channel
    CALL setupLCC(N_L, &
         compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
         idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
         bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)

    ! RyR channel
    CALL setupRyR(N_R, &
         compK_Rm, compK_Rp, &
         idxK_Rm, idxK_Rp, &
         bigU_R, irow_R, maxnkRm, maxnkRp)

    Mbig = irow_L*irow_R
    CALL getclustercompKleft(compK_Lv1, idxK_Lv1, &
         irow_L, irow_R, maxnkLv1, compKv1, idxKv1, maxnkv1)
    CALL getclustercompKleft(compK_Lv2, idxK_Lv2, &
         irow_L, irow_R, maxnkLv2, compKv2, idxKv2, maxnkv2)
    CALL getclustercompKleft(compK_Lm1, idxK_Lm1, &
         irow_L, irow_R, maxnkLm1, compKm1, idxKm1, maxnkm1)
    CALL getclustercompKright(compK_Rm, idxK_Rm, &
         irow_L, irow_R, maxnkRm, compKm2, idxKm2, maxnkm2)
    CALL getclustercompKleft(compK_Lp1, idxK_Lp1, &
         irow_L, irow_R, maxnkLp1, compKp1, idxKp1, maxnkp1)
    CALL getclustercompKright(compK_Rp, idxK_Rp, &
         irow_L, irow_R, maxnkRp, compKp2, idxKp2, maxnkp2)

    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(Mbig))

    IF ( ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
    ALLOCATE(U_Lo(Mbig))


    DO ii = 1, Mbig, irow_R
       U_Ro(ii:ii+irow_R-1) = bigU_R(1:irow_R)
       U_Lo(ii:ii+irow_R-1) = bigU_L((ii-1)/irow_R + 1)
    ENDDO
#if _DEBUG_ == 3
    WRITE(*,*) "U_Lo"
    WRITE(*,*) (U_Lo(ii), ii=1,Mbig)
#endif
    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
    IF (ALLOCATED(bigU_R)) DEALLOCATE(bigU_R)

    RETURN
  END SUBROUTINE makecompKs_3

  !========================================================
  ! Work with 6-state L-type
  ! return
  !   compK_Lv1, compK_Lv2,  idxK_Lv1, idxK_Lv2
  !   compK_Lm1, idxK_Lm1
  !   compK_Lp1,   idxK_Lp1
  !   bigU_L
  !   irow_L
  !   maxnkLv1, maxnkLv2
  !   maxnkLm1
  !   maxnkLp1
  !========================================================
!!$  SUBROUTINE setupLCC_6state(N_L, &
!!$       compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$       idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$       bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$    IMPLICIT NONE
!!$    !constant
!!$    INTEGER, PARAMETER:: mL=6 !mL = # of states for L-type channel
!!$    INTEGER, INTENT(in) :: N_L        ! # of LCC 
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
!!$         compK_Lm1, compK_Lp1
!!$    INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
!!$
!!$    !local
!!$    REAL(KIND=dp), DIMENSION(mL,mL) :: akLv1, akLv2, akLm1, akLp1
!!$    INTEGER, DIMENSION(mL) :: iu_L
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
!!$    INTEGER :: balls, bins, irows, icols, ii, jj
!!$    REAL(KIND=dp) :: akL25, akL16, akL34, akL21, akL56, akL23, akL54, akL52 &
!!$         , akL12, akL32, akL65, akL45, akL61, akL43 &
!!$         , rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1
!!$
!!$
!!$    !num_state_LCC = mL
!!$
!!$    ! Call ballbin to get possible state configurations      
!!$    balls=N_L
!!$    bins=mL
!!$
!!$    CALL count_ball2bin(balls, bins, iUi, irows, icols)
!!$    irow_L = irows ! # of L-type cluster-state
!!$                                ! NOTE: icols == bins
!!$                                ! iUi - matrix with each row denote a L-type group's state
!!$                                !       a row contains the number of L-type channels
!!$                                !     + in each L-type state
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "balls, bins: ", balls, bins
!!$    WRITE(*,*) "irows, icols: ", irows, icols
!!$    CALL print_mat(iUi, irows, icols)
!!$#endif
!!$
!!$                                ! with 6 state, set-up which state is open
!!$                                ! 0 = close
!!$                                ! 1 = open
!!$    iu_L = +0
!!$    iu_L(1) = +0
!!$    iu_L(2) = +0
!!$    iu_L(3) = +0
!!$    iu_L(4) = +0
!!$    iu_L(5) = 1
!!$    iu_L(6) = +0
!!$                                ! rate-transition parameters
!!$#if _USE_LCC_DATA_FROM_CONFIG_ == _YES
!!$                                !%v1-dependent
!!$    akL25 = kL25
!!$    akL16 = kL16
!!$    akL34 = kL34
!!$                                !%v2-dependent
!!$    akL21 = kL21
!!$    akL56 = kL56
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = kL23
!!$    akL54 = kL54
!!$                                !%m (Ca_ds-independent)
!!$    akL52 = kL52
!!$    akL12 = kL12
!!$    akL65 = kL65
!!$
!!$    akL32 = kL32
!!$    akL45 = kL45
!!$#else
!!$                                !%v1-dependent
!!$    akL25 = 350d0
!!$    akL16 = 350d0
!!$    akL34 = 350d0
!!$                                !%v2-dependent
!!$    akL21 = 10.d0
!!$    akL56 = 10.d0
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = 0.8d0
!!$    akL54 = 0.8d0
!!$                                !%m (Ca_ds-independent)
!!$    akL52 = 2000d0
!!$    akL12 = 10.d0
!!$    akL65 = akL12/10
!!$
!!$    akL32 = 10.d0
!!$    akL45 = akL32/10
!!$#endif
!!$    !!thermo
!!$    akL61 = (akL16*akL65*akL52*akL21)/(akL12*akL25*akL56)
!!$    akL43 = (akL34*akL45*akL52*akL23)/(akL32*akL25*akL54)
!!$                                ! zero-matrices
!!$    akLv1 = +0.0d0
!!$    akLv2 = +0.0d0
!!$    akLm1 = +0.0d0
!!$    akLp1 = +0.0d0
!!$                                ! akLv1, akLv2: voltage-dependent (for use in different functions)
!!$    akLv1(1,6) = akL16
!!$    akLv1(2,5) = akL25
!!$    akLv1(3,4) = akL34
!!$
!!$    akLv2(2,1) = akL21
!!$    akLv2(5,6) = akL56
!!$                                ! akLm1 : c_ds -independent
!!$                                ! akLp1 : c_ds -dependent
!!$    akLm1(6,1) = akL61
!!$    akLm1(1,2) = akL12
!!$    akLm1(3,2) = akL32
!!$    akLm1(5,2) = akL52
!!$    akLm1(4,3) = akL43
!!$    akLm1(4,5) = akL45
!!$    akLm1(6,5) = akL65
!!$
!!$    akLp1(2,3) = akL23
!!$    akLp1(5,4) = akL54
!!$                                !NOTE: For L-type, we don't use c_jsr dependent/independent, i.e. akLp2, akLm2
!!$
!!$
!!$    DO ii = 1,mL
!!$       rowsum_v1 = +0.d0
!!$       rowsum_v2 = +0.d0
!!$       rowsum_m1 = +0.d0
!!$       rowsum_p1 = +0.d0
!!$       DO jj = 1,mL
!!$          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
!!$          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
!!$          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
!!$          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
!!$       ENDDO
!!$       akLv1(ii,ii) = -rowsum_v1
!!$       akLv2(ii,ii) = -rowsum_v2
!!$       akLm1(ii,ii) = -rowsum_m1
!!$       akLp1(ii,ii) = -rowsum_p1
!!$    ENDDO
!!$#if _DEBUG_ == 1
!!$    PRINT *, "Here debug compK"
!!$    DO ii = 1, mL
!!$       DO jj = 1, mL
!!$          WRITE(*, "(F8.2)", ADVANCE='NO') akLv1(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(ii,jj)
!!$       END DO
!!$       WRITE(*,*)
!!$    END DO
!!$#endif
!!$
!!$    !size (m,m): akLv1, akLv2, akLm1, akLp1
!!$    !size (irow_L,m): iUi
!!$    !size (irow_L,0:maxnkL..) :: compK_L..
!!$    !size (irow_L,maxnkL..) :: idxK_L..
!!$    ! CALL getcompK_3(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1)
!!$    ! CALL getcompK_3(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2)
!!$    ! CALL getcompK_3(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1)
!!$    ! CALL getcompK_3(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1)
!!$    CALL getcompK_4(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1, N_L)
!!$    CALL getcompK_4(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2, N_L)
!!$    CALL getcompK_4(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1, N_L)
!!$    CALL getcompK_4(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1, N_L)
!!$
!!$                                ! Define bigU_L
!!$                                ! + vector that tells how many active L-types in each L-type cluster-state
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    ALLOCATE(bigU_L(irow_L))
!!$
!!$    DO ii = 1, irow_L
!!$       bigU_L(ii) = +0.d0
!!$       DO jj = 1, mL
!!$          IF (iu_L(jj).EQ.1) THEN
!!$             bigU_L(ii) = bigU_L(ii) + dfloat(iUi(ii,jj))
!!$          ENDIF
!!$       ENDDO
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE (*,*) "vector of # open state: bigU_L"
!!$    WRITE(*,*), (bigU_L(ii), ii = 1,irow_L)
!!$#endif
!!$                                !free
!!$    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)
!!$
!!$    RETURN
!!$  END SUBROUTINE setupLCC_6state
  !! Transpose version
!!$  SUBROUTINE setupLCC_6state_T(N_L, &
!!$       compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$       idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$       bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$    IMPLICIT NONE
!!$    !constant
!!$    INTEGER, PARAMETER:: mL=6 !mL = # of states for L-type channel
!!$    ! arguments
!!$    INTEGER, INTENT(in) :: N_L        ! # of LCC 
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
!!$         compK_Lm1, compK_Lp1
!!$    INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
!!$
!!$    !local
!!$    REAL(KIND=dp), DIMENSION(mL,mL) :: akLv1, akLv2, akLm1, akLp1
!!$    INTEGER, DIMENSION(mL) :: iu_L
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
!!$    INTEGER :: balls, bins, irows, icols, ii, jj
!!$    REAL(KIND=dp) :: akL25, akL16, akL34, akL21, akL56, akL23, akL54, akL52 &
!!$         , akL12, akL32, akL65, akL45, akL61, akL43 &
!!$         , rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1
!!$
!!$    !num_state_LCC = mL
!!$    
!!$    ! Call ballbin to get possible state configurations      
!!$    balls=N_L
!!$    bins=mL
!!$
!!$    CALL count_ball2bin(balls, bins, iUi, irows, icols)
!!$    irow_L = irows ! # of L-type cluster-state
!!$                                ! NOTE: icols == bins
!!$                                ! iUi - matrix with each row denote a L-type group's state
!!$                                !       a row contains the number of L-type channels
!!$                                !     + in each L-type state
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "balls, bins: ", balls, bins
!!$    WRITE(*,*) "irows, icols: ", irows, icols
!!$    CALL print_mat(iUi, irows, icols)
!!$#endif
!!$
!!$                                ! with 6 state, set-up which state is open
!!$                                ! 0 = close
!!$                                ! 1 = open
!!$    iu_L = +0
!!$    iu_L(1) = +0
!!$    iu_L(2) = +0
!!$    iu_L(3) = +0
!!$    iu_L(4) = +0
!!$    iu_L(5) = 1
!!$    iu_L(6) = +0
!!$                                ! rate-transition parameters
!!$#if _USE_LCC_DATA_FROM_CONFIG_ == _YES
!!$                                !%v1-dependent
!!$    akL25 = kL25
!!$    akL16 = kL16
!!$    akL34 = kL34
!!$                                !%v2-dependent
!!$    akL21 = kL21
!!$    akL56 = kL56
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = kL23
!!$    akL54 = kL54
!!$                                !%m (Ca_ds-independent)
!!$    akL52 = kL52
!!$    akL12 = kL12
!!$    akL65 = kL65
!!$
!!$    akL32 = kL32
!!$    akL45 = kL45
!!$#else
!!$                                !%v1-dependent
!!$    akL25 = 350d0
!!$    akL16 = 350d0
!!$    akL34 = 350d0
!!$                                !%v2-dependent
!!$    akL21 = 10.d0
!!$    akL56 = 10.d0
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = 0.8d0
!!$    akL54 = 0.8d0
!!$                                !%m (Ca_ds-independent)
!!$    akL52 = 2000d0
!!$    akL12 = 10.d0
!!$    akL65 = akL12/10
!!$
!!$    akL32 = 10.d0
!!$    akL45 = akL32/10
!!$#endif
!!$    !!thermo
!!$    akL61 = (akL16*akL65*akL52*akL21)/(akL12*akL25*akL56)
!!$    akL43 = (akL34*akL45*akL52*akL23)/(akL32*akL25*akL54)
!!$                                ! zero-matrices
!!$    akLv1 = +0.0d0
!!$    akLv2 = +0.0d0
!!$    akLm1 = +0.0d0
!!$    akLp1 = +0.0d0
!!$                                ! akLv1, akLv2: voltage-dependent (for use in different functions)
!!$    akLv1(1,6) = akL16
!!$    akLv1(2,5) = akL25
!!$    akLv1(3,4) = akL34
!!$
!!$    akLv2(2,1) = akL21
!!$    akLv2(5,6) = akL56
!!$                                ! akLm1 : c_ds -independent
!!$                                ! akLp1 : c_ds -dependent
!!$    akLm1(6,1) = akL61
!!$    akLm1(1,2) = akL12
!!$    akLm1(3,2) = akL32
!!$    akLm1(5,2) = akL52
!!$    akLm1(4,3) = akL43
!!$    akLm1(4,5) = akL45
!!$    akLm1(6,5) = akL65
!!$
!!$    akLp1(2,3) = akL23
!!$    akLp1(5,4) = akL54
!!$                                !NOTE: For L-type, we don't use c_jsr dependent/independent, i.e. akLp2, akLm2
!!$
!!$
!!$    DO ii = 1,mL
!!$       rowsum_v1 = +0.d0
!!$       rowsum_v2 = +0.d0
!!$       rowsum_m1 = +0.d0
!!$       rowsum_p1 = +0.d0
!!$       DO jj = 1,mL
!!$          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
!!$          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
!!$          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
!!$          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
!!$       ENDDO
!!$       akLv1(ii,ii) = -rowsum_v1
!!$       akLv2(ii,ii) = -rowsum_v2
!!$       akLm1(ii,ii) = -rowsum_m1
!!$       akLp1(ii,ii) = -rowsum_p1
!!$    ENDDO
!!$#if _DEBUG_ == 1
!!$    PRINT *, "Here debug compK"
!!$    DO ii = 1, mL
!!$       DO jj = 1, mL
!!$          WRITE(*, "(F8.2)", ADVANCE='NO') akLv1(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(ii,jj)
!!$                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(ii,jj)
!!$       END DO
!!$       WRITE(*,*)
!!$    END DO
!!$#endif
!!$
!!$                                !size (m,m): akLv1, akLv2, akLm1, akLp1
!!$                                !size (irow_L,m): iUi
!!$                                !size (0:maxnkL.., irow_L) :: compK_L..
!!$                                !size (0:maxnkL.., irow_L) :: idxK_L..
!!$    ! CALL getcompK_3(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1)
!!$    ! CALL getcompK_3(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2)
!!$    ! CALL getcompK_3(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1)
!!$    ! CALL getcompK_3(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1)
!!$    CALL getcompK_4_T(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1, N_L)
!!$    CALL getcompK_4_T(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2, N_L)
!!$    CALL getcompK_4_T(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1, N_L)
!!$    CALL getcompK_4_T(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1, N_L)
!!$
!!$                                ! Define bigU_L
!!$                                ! + vector that tells how many active L-types in each L-type cluster-state
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    ALLOCATE(bigU_L(irow_L))
!!$
!!$    DO ii = 1, irow_L
!!$       bigU_L(ii) = +0.d0
!!$       DO jj = 1, mL
!!$          IF (iu_L(jj).EQ.1) THEN
!!$             bigU_L(ii) = bigU_L(ii) + dfloat(iUi(ii,jj))
!!$          ENDIF
!!$       ENDDO
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE (*,*) "vector of # open state: bigU_L"
!!$    WRITE(*,*), (bigU_L(ii), ii = 1,irow_L)
!!$#endif
!!$                                !free
!!$    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)
!!$
!!$    RETURN
!!$  END SUBROUTINE setupLCC_6state_T

!!$  SUBROUTINE setupLCC_7state(N_L, &
!!$       compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$       idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$       bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$    IMPLICIT NONE
!!$    !constant
!!$    INTEGER, PARAMETER:: mL=7 !mL = # of states for L-type channel
!!$
!!$    !arguments
!!$    INTEGER, INTENT(in) :: N_L        ! # of LCC 
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
!!$         compK_Lm1, compK_Lp1
!!$    INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
!!$
!!$    !local
!!$    REAL(KIND=dp), DIMENSION(mL,mL) :: akLv1, akLv2, akLm1, akLp1
!!$    INTEGER, DIMENSION(mL) :: iu_L
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
!!$    INTEGER :: balls, bins, irows, icols, ii, jj
!!$    REAL(KIND=dp) :: akL25, akL16, akL34, akL21, akL56, akL23, akL54, akL52 &
!!$         , akL12, akL32, akL65, akL45, akL61, akL43 &
!!$         , rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1, &
!!$         akL76, akL74, akL75, akL57, akL67, akL47
!!$
!!$    !num_state_LCC = mL
!!$    
!!$    ! Call ballbin to get possible state configurations      
!!$    balls=N_L
!!$    bins=mL
!!$
!!$    CALL count_ball2bin(balls, bins, iUi, irows, icols)
!!$    irow_L = irows ! # of L-type cluster-state
!!$    ! NOTE: icols == bins
!!$    ! iUi - matrix with each row denote a L-type group's state
!!$    !       a row contains the number of L-type channels
!!$    !     + in each L-type state
!!$
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "balls, bins: ", balls, bins
!!$    WRITE(*,*) "irows, icols: ", irows, icols
!!$    CALL print_mat(iUi, irows, icols)
!!$    STOP
!!$#endif
!!$
!!$    ! with 7 state, set-up which state is open
!!$    ! 0 = close, inactive, ...
!!$    ! 1 = open
!!$    iu_L = +0
!!$    iu_L(7) = 1
!!$    ! rate-transition parameters
!!$#if _USE_LCC_DATA_FROM_CONFIG_ == _YES
!!$    !%v1-dependent
!!$    akL25 = kL25
!!$    akL16 = kL16
!!$    akL34 = kL34
!!$                                !%v2-dependent
!!$    akL21 = kL21
!!$    akL76 = kL76
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = kL23
!!$    akL74 = kL74
!!$                                !%m (Ca_ds-independent)
!!$    akL12 = kL12
!!$    akL32 = kL32
!!$    akL75 = kL75
!!$    akL61 = kL61
!!$    akL52 = kL52
!!$    akL43 = kL43
!!$    
!!$    ! opening
!!$    akL57 = kL57
!!$#else    
!!$    !%v1-dependent
!!$    akL25 = 700d0
!!$    akL16 = 700d0
!!$    akL34 = 700d0
!!$                                !%v2-dependent
!!$    akL21 = 25.d0
!!$    akL76 = 25.d0
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = 0.75d0
!!$    akL74 = 0.75d0
!!$                                !%m (Ca_ds-independent)
!!$    akL12  =  10d0
!!$    akL32  =  10d0
!!$    akL75  =  1300d0
!!$    akL61  =  500d0
!!$    akL52  =  500d0
!!$    akL43  =  500d0
!!$    
!!$    ! opening
!!$    akL57 = 500
!!$
!!$#endif
!!$    !!thermo
!!$    akL67 = (akL76*akL61*akL12*akL25*akL57)/(akL75*akL52*akL21*akL16)
!!$    akL47 = (akL74*akL43*akL32*akL25*akL57)/(akL75*akL52*akL23*akL34)
!!$
!!$    WRITE(log_string, *) "akL67, akL47 : ", akL67, akL47
!!$    CALL write2log(log_string)
!!$    
!!$    ! zero-matrices
!!$    akLv1 = +0.0d0
!!$    akLv2 = +0.0d0
!!$    akLm1 = +0.0d0
!!$    akLp1 = +0.0d0
!!$    ! akLv1: Vm-dependent activation 
!!$    ! akLv2: Vm-dependent inactivation
!!$    akLv1(2,5) = akL25
!!$    akLv1(1,6) = akL16
!!$    akLv1(3,4) = akL34
!!$
!!$    akLv2(2,1) = akL21
!!$    akLv2(7,6) = akL76
!!$
!!$    ! akLm1 : c_ds -independent
!!$    ! akLp1 : c_ds -dependent
!!$    akLm1(1,2) = akL12
!!$    akLm1(3,2) = akL32
!!$    akLm1(4,3) = akL43
!!$    akLm1(4,7) = akL47
!!$    akLm1(5,2) = akL52
!!$    akLm1(5,7) = akL57
!!$    akLm1(6,1) = akL61
!!$    akLm1(6,7) = akL67
!!$    akLm1(7,5) = akL75
!!$    
!!$    akLp1(2,3) = akL23
!!$    akLp1(7,4) = akL74
!!$    !NOTE: For L-type, we don't use c_jsr dependent/independent, i.e. akLp2, akLm2
!!$
!!$
!!$    DO ii = 1,mL
!!$       rowsum_v1 = +0.d0
!!$       rowsum_v2 = +0.d0
!!$       rowsum_m1 = +0.d0
!!$       rowsum_p1 = +0.d0
!!$       DO jj = 1,mL
!!$          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
!!$          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
!!$          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
!!$          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
!!$       ENDDO
!!$       akLv1(ii,ii) = -rowsum_v1
!!$       akLv2(ii,ii) = -rowsum_v2
!!$       akLm1(ii,ii) = -rowsum_m1
!!$       akLp1(ii,ii) = -rowsum_p1
!!$    ENDDO
!!$#if _DEBUG_ == 1
!!$    PRINT *, "Here debug compK"
!!$    DO ii = 1, mL
!!$       DO jj = 1, mL
!!$          WRITE(*, "(F8.2)", ADVANCE='NO') akLv1(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(ii,jj)
!!$       END DO
!!$       WRITE(*,*)
!!$    END DO
!!$#endif
!!$
!!$    !size (m,m): akLv1, akLv2, akLm1, akLp1
!!$    !size (irow_L,m): iUi
!!$    !size (irow_L,0:maxnkL..) :: compK_L..
!!$    !size (irow_L,maxnkL..) :: idxK_L..
!!$    CALL getcompK_4(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1, N_L)
!!$    CALL getcompK_4(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2, N_L)
!!$    CALL getcompK_4(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1, N_L)
!!$    CALL getcompK_4(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1, N_L)
!!$
!!$    ! Define bigU_L
!!$    ! + vector that tells how many active L-types in each L-type cluster-state
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    ALLOCATE(bigU_L(irow_L))
!!$
!!$    DO ii = 1, irow_L
!!$       bigU_L(ii) = +0.d0
!!$       DO jj = 1, mL
!!$          IF (iu_L(jj).EQ.1) THEN
!!$             bigU_L(ii) = bigU_L(ii) + dfloat(iUi(ii,jj))
!!$          ENDIF
!!$       ENDDO
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE (*,*) "vector of # open state: bigU_L"
!!$    WRITE(*,*), (bigU_L(ii), ii = 1,irow_L)
!!$#endif
!!$                                !free
!!$    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)
!!$
!!$    RETURN
!!$  END SUBROUTINE setupLCC_7state
  
  !! Transpose version
!!$  SUBROUTINE setupLCC_7state_T(N_L, &
!!$       compK_Lv1, compK_Lv2, compK_Lm1, compK_Lp1, &
!!$       idxK_Lv1, idxK_Lv2, idxK_Lm1, idxK_Lp1, &
!!$       bigU_L, irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1)
!!$
!!$    IMPLICIT NONE
!!$    !constant
!!$    INTEGER, PARAMETER:: mL=7 !mL = # of states for L-type channel
!!$
!!$    !arguments
!!$    INTEGER, INTENT(in) :: N_L        ! # of LCC 
!!$    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bigU_L  ! # of open LCC at a cluster-state
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK_Lv1, idxK_Lv2, idxK_Lm1, &
!!$         idxK_Lp1
!!$    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK_Lv1, compK_Lv2, &
!!$         compK_Lm1, compK_Lp1
!!$    INTEGER, INTENT(OUT) :: irow_L, maxnkLv1, maxnkLv2, maxnkLm1, maxnkLp1
!!$
!!$    !local
!!$    REAL(KIND=dp), DIMENSION(mL,mL) :: akLv1, akLv2, akLm1, akLp1
!!$    INTEGER, DIMENSION(mL) :: iu_L
!!$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
!!$    INTEGER :: balls, bins, irows, icols, ii, jj
!!$    REAL(KIND=dp) :: akL25, akL16, akL34, akL21, akL56, akL23, akL54, akL52 &
!!$         , akL12, akL32, akL65, akL45, akL61, akL43 &
!!$         , rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1, &
!!$         akL76, akL74, akL75, akL57, akL67, akL47
!!$
!!$    !num_state_LCC = mL
!!$    
!!$    ! Call ballbin to get possible state configurations      
!!$    balls=N_L
!!$    bins=mL
!!$
!!$    CALL count_ball2bin(balls, bins, iUi, irows, icols)
!!$    irow_L = irows ! # of L-type cluster-state
!!$    ! NOTE: icols == bins
!!$    ! iUi - matrix with each row denote a L-type group's state
!!$    !       a row contains the number of L-type channels
!!$    !     + in each L-type state
!!$
!!$#if _DEBUG_ == 3
!!$    WRITE(*,*) "balls, bins: ", balls, bins
!!$    WRITE(*,*) "irows, icols: ", irows, icols
!!$    CALL print_mat(iUi, irows, icols)
!!$    STOP
!!$#endif
!!$
!!$    ! with 7 state, set-up which state is open
!!$    ! 0 = close, inactive, ...
!!$    ! 1 = open
!!$    iu_L = +0
!!$    iu_L(7) = 1
!!$    ! rate-transition parameters
!!$#if _USE_LCC_DATA_FROM_CONFIG_ == _YES
!!$    !%v1-dependent
!!$    akL25 = kL25
!!$    akL16 = kL16
!!$    akL34 = kL34
!!$                                !%v2-dependent
!!$    akL21 = kL21
!!$    akL76 = kL76
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = kL23
!!$    akL74 = kL74
!!$                                !%m (Ca_ds-independent)
!!$    akL12 = kL12
!!$    akL32 = kL32
!!$    akL75 = kL75
!!$    akL61 = kL61
!!$    akL52 = kL52
!!$    akL43 = kL43
!!$    
!!$    ! opening
!!$    akL57 = kL57
!!$#else    
!!$    !%v1-dependent
!!$    akL25 = 700d0
!!$    akL16 = 700d0
!!$    akL34 = 700d0
!!$                                !%v2-dependent
!!$    akL21 = 25.d0
!!$    akL76 = 25.d0
!!$                                !%c1 (Ca_ds-dependent)
!!$    akL23 = 0.75d0
!!$    akL74 = 0.75d0
!!$                                !%m (Ca_ds-independent)
!!$    akL12 = 10d0
!!$    akL32 = 10d0
!!$    akL75 = 1300d0
!!$    akL61 = 500d0
!!$    akL52 = 500d0
!!$    akL43 = 500d0
!!$    
!!$    ! opening
!!$    akL57 = 500
!!$
!!$#endif
!!$    !!thermo
!!$    akL67 = (akL76*akL61*akL12*akL25*akL57)/(akL75*akL52*akL21*akL16)
!!$    akL47 = (akL74*akL43*akL32*akL25*akL57)/(akL75*akL52*akL23*akL34)
!!$    
!!$    WRITE(log_string, *) "akL67, akL47 : ", akL67, akL47
!!$    CALL write2log(log_string)
!!$
!!$    ! zero-matrices
!!$    akLv1 = +0.0d0
!!$    akLv2 = +0.0d0
!!$    akLm1 = +0.0d0
!!$    akLp1 = +0.0d0
!!$    ! akLv1: Vm-dependent activation 
!!$    ! akLv2: Vm-dependent inactivation
!!$    akLv1(2,5) = akL25
!!$    akLv1(1,6) = akL16
!!$    akLv1(3,4) = akL34
!!$
!!$    akLv2(2,1) = akL21
!!$    akLv2(7,6) = akL76
!!$
!!$    ! akLm1 : c_ds -independent
!!$    ! akLp1 : c_ds -dependent
!!$    akLm1(1,2) = akL12
!!$    akLm1(3,2) = akL32
!!$    akLm1(4,3) = akL43
!!$    akLm1(4,7) = akL47
!!$    akLm1(5,2) = akL52
!!$    akLm1(5,7) = akL57
!!$    akLm1(6,1) = akL61
!!$    akLm1(6,7) = akL67
!!$    akLm1(7,5) = akL75
!!$    
!!$    akLp1(2,3) = akL23
!!$    akLp1(7,4) = akL74
!!$    !NOTE: For L-type, we don't use c_jsr dependent/independent, i.e. akLp2, akLm2
!!$
!!$
!!$    DO ii = 1,mL
!!$       rowsum_v1 = +0.d0
!!$       rowsum_v2 = +0.d0
!!$       rowsum_m1 = +0.d0
!!$       rowsum_p1 = +0.d0
!!$       DO jj = 1,mL
!!$          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
!!$          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
!!$          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
!!$          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
!!$       ENDDO
!!$       akLv1(ii,ii) = -rowsum_v1
!!$       akLv2(ii,ii) = -rowsum_v2
!!$       akLm1(ii,ii) = -rowsum_m1
!!$       akLp1(ii,ii) = -rowsum_p1
!!$    ENDDO
!!$#if _DEBUG_ == 1
!!$    PRINT *, "Here debug compK"
!!$    DO ii = 1, mL
!!$       DO jj = 1, mL
!!$          WRITE(*, "(F8.2)", ADVANCE='NO') akLv1(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(ii,jj)
!!$          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(ii,jj)
!!$       END DO
!!$       WRITE(*,*)
!!$    END DO
!!$#endif
!!$
!!$    !size (m,m): akLv1, akLv2, akLm1, akLp1
!!$    !size (irow_L,m): iUi
!!$    !size (0:maxnkL.., irow_L) :: compK_L..
!!$    !size (0:maxnkL.., irow_L) :: idxK_L..
!!$    ! CALL getcompK_3(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1)
!!$    ! CALL getcompK_3(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2)
!!$    ! CALL getcompK_3(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1)
!!$    ! CALL getcompK_3(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1)
!!$    CALL getcompK_4_T(akLv1, mL, iUi, irow_L, compK_Lv1, idxK_Lv1, maxnkLv1, N_L)
!!$    CALL getcompK_4_T(akLv2, mL, iUi, irow_L, compK_Lv2, idxK_Lv2, maxnkLv2, N_L)
!!$    CALL getcompK_4_T(akLm1, mL, iUi, irow_L, compK_Lm1, idxK_Lm1, maxnkLm1, N_L)
!!$    CALL getcompK_4_T(akLp1, mL, iUi, irow_L, compK_Lp1, idxK_Lp1, maxnkLp1, N_L)
!!$
!!$    ! Define bigU_L
!!$    ! + vector that tells how many active L-types in each L-type cluster-state
!!$    IF (ALLOCATED(bigU_L)) DEALLOCATE(bigU_L)
!!$    ALLOCATE(bigU_L(irow_L))
!!$
!!$    DO ii = 1, irow_L
!!$       bigU_L(ii) = +0.d0
!!$       DO jj = 1, mL
!!$          IF (iu_L(jj).EQ.1) THEN
!!$             bigU_L(ii) = bigU_L(ii) + dfloat(iUi(ii,jj))
!!$          ENDIF
!!$       ENDDO
!!$    ENDDO
!!$#if _DEBUG_ == 3
!!$    WRITE (*,*) "vector of # open state: bigU_L"
!!$    WRITE(*,*), (bigU_L(ii), ii = 1,irow_L)
!!$#endif
!!$                                !free
!!$    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)
!!$
!!$    RETURN
!!$  END SUBROUTINE setupLCC_7state_T

  !==============================================
  ! version 3 - work with full ECC model
  ! return
  !   compK, idxK
  !   maxnkl
  !   iBound
  !==============================================      
  SUBROUTINE getcompK_3(ak, m, iUi, irows, compK, idxK, maxnkl)
    IMPLICIT NONE

    REAL(KIND=dp), DIMENSION(m,m) :: ak ! rate transition from one channel state to another
    INTEGER, DIMENSION(irows,m) :: iUi 
    INTEGER, INTENT(IN) :: m, & ! # of channel states
         irows  ! # of cluster states
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK ! matrix of rate transition from one cluster-state
                                ! to another (keep only non-zero values)
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  idxK          ! keep index of the state
    INTEGER, INTENT(OUT) :: maxnkl
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    INTEGER :: icount, ii, jj, initial_state, final_state
    REAL(kind=dp) :: asum

                                !  Find maximum # of possible neighbors
                                !+ maximum number of state transition is the # non-zero non-diagonal elements
    icount = +0
    DO ii = 1, m 
       mask =  ak(ii,:) .GT. EPS
       icount = icount + COUNT(mask)
    ENDDO
    !! The above value work for N_L > m
    !! If N_L < m, a better boundary would be
    !!+ count the number non-zero non-diagonal elements from
    !!+ N_L lines with maximum number of non-zero values
                                !(to be implemented)
    !! Finally
    maxnkl = icount + 1     ! +1 for state unchanged

    !!     Gather information from big K matrices
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
                                ! the zero-th column is for consistency with idxK matrix
    ALLOCATE(compK(irows, 0:maxnkl))       

    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
                                ! As different state has different number of state transition
                                !+ the zero-th column keep the value indicating the true # possible state transition,
                                !+ excluding the current state
    ALLOCATE(idxK(irows, 0:maxnkl))

    !! default values
    compK = +0.d0
    idxK = +0

                                ! check all pairs of cluster state, see which pair is valid for a neighbor in state transition
    DO ii = 1, irows
       icount = 1
       asum = +0.d0
       DO jj = 1,irows
          tmpvec = iUi(ii,1:m) - iUi(jj,1:m)
                                ! valid for state change = yes
                                !+ when only one element = 1, one element = -1, the others = 0
          IF (SUM(ABS(tmpvec)) .EQ. 2) THEN 
             initial_state = MAXLOC(tmpvec, 1)  ! location of element with value 1
             final_state = MINLOC(tmpvec,1)     ! location of element with value -1
             IF (ak(initial_state, final_state) .GT. EPS) THEN ! positive rate transition 
                icount = icount + 1
                compK(ii,icount) = dfloat(iUi(ii,initial_state)) &
                     *ak(initial_state,final_state)
                idxK(ii,icount) = jj
                asum = asum + compK(ii, icount)
             ENDIF
          ENDIF
       ENDDO
       idxK(ii, 1) = ii
       ! Find the negative value at each row of aK
       !+ and keep it at the first column
       !+ one minus to this value would the probability that the state doesn't change
       IF (asum .GT. EPS) THEN
          compK(ii, 1) = -asum
       ENDIF
       ! keep the true number of neighbors at state ii (including itself)
       IF (icount > 1) THEN
          idxK(ii, 0) = icount
       ELSE
          idxK(ii, 0) = +0
       ENDIF
    ENDDO
    RETURN
  END SUBROUTINE getcompK_3
                                ! Transpose version 
  SUBROUTINE getcompK_3_T(ak, m, iUi, irows, compK, idxK, maxnkl)
    IMPLICIT NONE

    REAL(KIND=dp), DIMENSION(m,m) :: ak ! rate transition from one channel state to another
    INTEGER, DIMENSION(irows,m) :: iUi 
    INTEGER, INTENT(IN) :: m, & ! # of channel states
         irows  ! # of cluster states
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK ! matrix of rate transition from one cluster-state
                                ! to another (keep only non-zero values)
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  idxK          ! keep index of the state
    INTEGER, INTENT(OUT) :: maxnkl
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    INTEGER :: icount, ii, jj, initial_state, final_state
    REAL(kind=dp) :: asum

    !  Find maximum # of possible neighbors
    !+ maximum number of state transition is the # non-zero non-diagonal elements
    icount = +0
    DO ii = 1, m 
       mask =  ak(ii,:) .GT. EPS
       icount = icount + COUNT(mask)
    ENDDO
    !! The above value work for N_L > m
    !! If N_L < m, a better boundary would be
    !!+ count the number non-zero non-diagonal elements from
    !!+ N_L lines with maximum number of non-zero values
                                !(to be implemented)
    !! Finally
    maxnkl = icount + 1     ! +1 for state unchanged

    !!     Gather information from big K matrices
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
                                ! the zero-th column is for consistency with idxK matrix
    ALLOCATE(compK(0:maxnkl, irows))       

    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
    ! As different state has different number of state transition
    !+ the zero-th column keep the value indicating the true # possible state transition,
    !+ excluding the current state
    ALLOCATE(idxK(0:maxnkl, irows))

    !! default values
    compK = +0.d0
    idxK = +0

                                ! check all pairs of cluster state, see which pair is valid for a neighbor in state transition
    DO ii = 1, irows
       icount = 1
       asum = +0.d0
       DO jj = 1,irows
          tmpvec = iUi(ii,1:m) - iUi(jj,1:m)
                                ! valid for state change = yes
                                !+ when only one element = 1, one element = -1, the others = 0
          IF (SUM(ABS(tmpvec)) .EQ. 2) THEN 
             initial_state = MAXLOC(tmpvec, 1)  ! location of element with value 1
             final_state = MINLOC(tmpvec,1)     ! location of element with value -1
             IF (ak(initial_state, final_state) .GT. EPS) THEN ! positive rate transition 
                icount = icount + 1
                compK(icount, ii) = dfloat(iUi(ii,initial_state)) &
                     *ak(initial_state,final_state)
                idxK(icount, ii) = jj
                asum = asum + compK(icount, ii)
             ENDIF
          ENDIF
       ENDDO
       idxK(1, ii) = ii
       ! Find the negative value at each row of aK
       !+ and keep it at the first column
       !+ one minus to this value would the probability that the state doesn't change
       IF (asum .GT. EPS) THEN
          compK(1, ii) = -asum
       ENDIF
                                ! keep the true number of neighbors at state ii (including itself)
       IF (icount > 1) THEN
          idxK(0, ii) = icount
       ELSE
          idxK(0, ii) = +0
       ENDIF
    ENDDO
    RETURN
  END SUBROUTINE getcompK_3_T

  !!======================================
  !! Find the Kronecker:
  !!   compK = compK_L (X) Identity
  !! return
  !!    compK, idxK
  !!    maxnk
  !!======================================
  SUBROUTINE getclustercompKleft(compK_L, idxK_L, &
       irow_L, irow_R, maxnkL, compK, idxK, maxnk)
    IMPLICIT NONE

    REAL(KIND=DP), DIMENSION(irow_L,0:maxnkL) :: compK_L
    INTEGER, DIMENSION(irow_L,0:maxnkL) :: idxK_L
    INTEGER, INTENT(IN) :: irow_L, irow_R, maxnkL
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK
    INTEGER, INTENT(OUT) :: maxnk

                                ! local
    INTEGER :: Mbig
    INTEGER :: ii, jj, kk

    Mbig = irow_L * irow_R
    maxnk = maxnkL

                                ! (X) = kronecker
                                ! compK_L (X) Identity
    ALLOCATE(compK(Mbig,0:maxnk))
    ALLOCATE(idxK(Mbig,0:maxnk))
    compK = +0.0d0
    idxK = +0
    DO ii = 1, irow_L
       DO kk = 1, maxnk
          IF (ABS(compK_L(ii,kk)) .GT. EPS) THEN
             compK((ii-1)*irow_R + 1 : ii*irow_R, kk) = compK_L(ii,kk)
             idxK((ii-1)*irow_R + 1 : ii*irow_R, 0) = idxK_L(ii,0)
             DO jj=1, irow_R
                idxK((ii-1)*irow_R + jj, kk) = (idxK_L(ii,kk)-1)*irow_R + jj
             ENDDO
          ENDIF
       ENDDO
    ENDDO
  END SUBROUTINE getclustercompKleft
  SUBROUTINE getclustercompKleft_T(compK_L, idxK_L, &
       irow_L, irow_R, maxnkL, compK, idxK, maxnk)
    IMPLICIT NONE

    REAL(KIND=DP), DIMENSION(0:maxnkL, irow_L) :: compK_L
    INTEGER, DIMENSION(0:maxnkL, irow_L) :: idxK_L
    INTEGER, INTENT(IN) :: irow_L, irow_R, maxnkL
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK
    INTEGER, INTENT(OUT) :: maxnk

                                ! local
    INTEGER :: Mbig
    INTEGER :: ii, jj, kk

    Mbig = irow_L * irow_R
    maxnk = maxnkL

                                ! (X) = kronecker
                                ! compK_L (X) Identity
                                ! NOTE: Here, we transpose the matrix
    ALLOCATE(compK(0:maxnk, Mbig))
    ALLOCATE(idxK(0:maxnk, Mbig))
    compK = +0.0d0
    idxK = +0
    DO ii = 1, irow_L
       DO kk = 1, maxnk
          IF (ABS(compK_L(kk, ii)) .GT. EPS) THEN
             compK(kk, (ii-1)*irow_R + 1 : ii*irow_R) = compK_L(kk, ii)
             idxK(0, (ii-1)*irow_R + 1 : ii*irow_R) = idxK_L(0, ii)
             DO jj=1, irow_R
                idxK(kk, (ii-1)*irow_R + jj) = (idxK_L(kk, ii)-1)*irow_R + jj
             ENDDO
          ENDIF
       ENDDO
    ENDDO
  END SUBROUTINE getclustercompKleft_T

  !!======================================
  !! Find the Kronecker:
  !!   compK = Identity (X) compK_R
  !! return
  !!    compK, idxK
  !!    maxnk
  !!======================================
  SUBROUTINE getclustercompKright(compK_R, idxK_R, &
       irow_L, irow_R, maxnkR, compK, idxK, maxnk)
    IMPLICIT NONE

    REAL(KIND=DP), DIMENSION(irow_R,0:maxnkR) :: compK_R
    INTEGER, DIMENSION(irow_R,0:maxnkR) :: idxK_R
    INTEGER, INTENT(IN) :: irow_L, irow_R, maxnkR
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK
    INTEGER, INTENT(OUT) :: maxnk
                                ! local
    INTEGER :: Mbig
    INTEGER :: ii, jj, kk

    Mbig = irow_L * irow_R
    maxnk = maxnkR

                                ! (X) = kronecker 
                                ! Identity (X) compK_R
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
    ALLOCATE(compK(Mbig,0:maxnk))
    ALLOCATE(idxK(Mbig,0:maxnk))
    compK = +0.0d0
    idxK = +0
    DO ii = 1, Mbig, irow_R
       compK(ii:ii+irow_R-1, 1:maxnkR) = compK_R(1:irow_R, 1:maxnkR)
       idxK(ii:ii+irow_R-1, 0) = idxK_R(1:irow_R, 0)
    ENDDO
    DO ii = 1, irow_R
       DO jj = 1, maxnkR
          IF (ABS(compK_R(ii,jj)) .GT. EPS) THEN
             DO kk = 1, Mbig, irow_R
                idxK(kk+ii-1, jj) = kk-1 + idxK_R(ii,jj)
             ENDDO
          ENDIF
       ENDDO
    ENDDO
  END SUBROUTINE getclustercompKright
  SUBROUTINE getclustercompKright_T(compK_R, idxK_R, &
       irow_L, irow_R, maxnkR, compK, idxK, maxnk)
    IMPLICIT NONE

    REAL(KIND=DP), DIMENSION(0:maxnkR, irow_R) :: compK_R
    INTEGER, DIMENSION(0:maxnkR, irow_R) :: idxK_R
    INTEGER, INTENT(IN) :: irow_L, irow_R, maxnkR
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: idxK
    INTEGER, INTENT(OUT) :: maxnk
                                ! local
    INTEGER :: Mbig
    INTEGER :: ii, jj, kk

    Mbig = irow_L * irow_R
    maxnk = maxnkR

                                ! (X) = kronecker 
                                ! Identity (X) compK_R
                                !NOTE: Here transpose the matrix
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
    ALLOCATE(compK(0:maxnk, Mbig))
    ALLOCATE(idxK(0:maxnk, Mbig))
    compK = +0.0d0
    idxK = +0
    DO ii = 1, Mbig, irow_R
       compK(1:maxnkR, ii:ii+irow_R-1) = compK_R(1:maxnkR, 1:irow_R)
       idxK(0, ii:ii+irow_R-1) = idxK_R(0, 1:irow_R)
    ENDDO
    DO ii = 1, irow_R
       DO jj = 1, maxnkR
          IF (ABS(compK_R(jj, ii)) .GT. EPS) THEN
             DO kk = 1, Mbig, irow_R
                idxK(jj, kk+ii-1) = kk-1 + idxK_R(jj, ii)
             ENDDO
          ENDIF
       ENDDO
    ENDDO
  END SUBROUTINE getclustercompKright_T

  !!==============================
  !! used in makeKs
  !! return the transition matrix of a cluster of ion channels
  !!  given the transition matrix of a single channel
  !!==============================
  SUBROUTINE getbigK(ak,bigK,m,irows,icols,iUi)
    IMPLICIT NONE
    INTEGER, INTENT(IN) ::  m
    INTEGER, INTENT(INOUT) :: irows,icols
    REAL(KIND=dp), DIMENSION(m,m) ::  ak
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: bigK
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    REAL(KIND=dp) :: asum
    INTEGER :: i, j, ii, jj, initial_state, final_state

    IF (ALLOCATED(bigK)) DEALLOCATE(bigK)
    ALLOCATE(bigK(irows, irows))
    icols = irows

    DO i = 1,irows
       DO j = 1,irows
          bigk(i,j) = 0.0d0
       ENDDO
    ENDDO

    DO ii = 1,irows
       asum  =  0.0d0
       DO jj = 1,irows
          tmpvec = iUi(ii,1:m) - iUi(jj,1:m)
                                ! valid for state change = yes
                                !+ when only one element = 1, one element = -1, the others = 0
          IF (SUM(ABS(tmpvec)) .EQ. 2) THEN 
             initial_state = MAXLOC(tmpvec, 1)  ! location of element with value 1
             final_state = MINLOC(tmpvec,1)     ! location of element with value -1
             IF (ak(initial_state, final_state) .GT. EPS) THEN ! positive rate transition 
                bigK(ii,jj) = dfloat(iUi(ii,initial_state))*ak(initial_state,final_state)
                asum = asum + bigK(ii,jj)
             ENDIF
          ENDIF
       ENDDO
       ! Find the negative value at each row of aK
       !+ and keep it at the first column
       !+ one minus to this value would the probability that the state doesn't change
       IF (asum .GT. EPS) THEN
          bigK(ii, ii) = -asum
       ENDIF
    ENDDO
    RETURN
  END SUBROUTINE getbigK

  !!=========================================
  !! use for testing purpose (not efficient)
  !!   - not update data
  !!=========================================
  SUBROUTINE makeKs(aKv1,aKv2,aKm1,aKp1,aKm2,aKp2,U_Lo,U_Ro,N_L,N_R,irow_L,irow_R)
    use math_utility
    USE matrix_utility
                                !      implicit double precision (a-h,o-z)
    IMPLICIT NONE
    !! constant
    INTEGER, PARAMETER:: mL=6 !mL = # of states for L-type channel
    INTEGER, PARAMETER:: mR=2 !mR = # of states for RyR channel
    !! arguments
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKv1  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKv2  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKm1 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKm2 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKp1  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aKp2  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: U_Lo
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: U_Ro
    INTEGER :: N_L, N_R, irow_L, irow_R
    !! local
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: e_L
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: e_R
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bigU_L
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bigU_R
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: UI_L
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: UI_R
    INTEGER balls,bins,irows,icols
    INTEGER :: iu_L, iU_R
    DIMENSION iu_L(mL),iU_R(mR)
    DOUBLE PRECISION :: akLv1, akLv2, akLm1, akLp1, akRm2, akRp2
    DIMENSION akLv1(mL,mL),akLv2(mL,mL),akLm1(mL,mL),akLp1(mL,mL)
    DIMENSION akRm2(mR,mR),akRp2(mR,mR)
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Lv1  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Lv2  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Lm1 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Lp1 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Rm2 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bigK_Rp2 
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: eye  
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aI_L 
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: aI_R 
    INTEGER :: ii, jj, i, j
    REAL(KIND=dp) :: akL25, akL16, akL34, akL21, akL56, akL23, akL54, akL52 &
         , akL12, akL32, akL65, akL45, akL61, akL43, akR12, akR21 &
         , rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1

    !! L-type channel

    !! Call ballbin to get possible state configurations
    balls = N_L
    bins = mL

    CALL count_ball2bin(balls, bins, iUi, irows, icols)
    irow_L = irows ! # of L-type cluster-state
    ! NOTE: icols == bins
    ! iUi(irows x icols) - matrix with each row denote a L-type group's state
    !       by telling the number of L-type channels
    !     + in each L-type channel state
    
    ! with 6 state, set-up which state is open
    ! 0 = close
    ! 1 = open
    iu_L = +0
    iu_L(1) = +0
    iu_L(2) = +0
    iu_L(3) = +0
    iu_L(4) = +0
    iu_L(5) = 1
    iu_L(6) = +0
                                ! rate-transition parameters
                                !%v1
    akL25 = 350d0
    akL16 = 350d0
    akL34 = 350d0
                                !%v2
    akL21 = 10.d0
    akL56 = 10.d0
                                !%c1
    akL23 = 0.8d0
    akL54 = 0.8d0
                                !%m
    akL52 = 2000d0
    akL12 = 10.d0
    akL65 = akL12/10

    akL32 = 10.d0
    akL45 = akL32/10

    akL61 = (akL16*akL65*akL52*akL21)/(akL12*akL25*akL56)
    akL43 = (akL34*akL45*akL52*akL23)/(akL32*akL25*akL54)
                                ! zero-matrices
    akLv1 = +0.0d0
    akLv2 = +0.0d0
    akLm1 = +0.0d0
    akLp1 = +0.0d0
                                ! akLv1, akLv2: voltage-dependent (for use in different functions)
    akLv1(1,6) = akL16
    akLv1(2,5) = akL25
    akLv1(3,4) = akL34

    akLv2(2,1) = akL21
    akLv2(5,6) = akL56
                                ! akLm1 : c_ds -independent
                                ! akLp1 : c_ds -dependent
    akLm1(6,1) = akL61
    akLm1(1,2) = akL12
    akLm1(3,2) = akL32
    akLm1(5,2) = akL52
    akLm1(4,3) = akL43
    akLm1(4,5) = akL45
    akLm1(6,5) = akL65

    akLp1(2,3) = akL23
    akLp1(5,4) = akL54
                                !NOTE: For L-type, we don't use c_jsr dependent/independent, i.e. akLp2, akLm2

    DO ii = 1,mL
       rowsum_v1 = +0.d0
       rowsum_v2 = +0.d0
       rowsum_m1 = +0.d0
       rowsum_p1 = +0.d0
       DO jj = 1,mL
          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
       ENDDO
       akLv1(ii,ii) = -rowsum_v1
       akLv2(ii,ii) = -rowsum_v2
       akLm1(ii,ii) = -rowsum_m1
       akLp1(ii,ii) = -rowsum_p1
    ENDDO
#if _DEBUG_ == 1
    PRINT *, "Here debug makeKs"
    DO I = 1, mL
       DO J = 1, mL
          WRITE(*, "(F8.2)", ADVANCE='NO') akLv1(I,J)
          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(I,J)
          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(I,J)
          !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(I,J)
       END DO
       WRITE(*,*)
    END DO
    return
#endif

    !! For L-type cluster
    CALL getbigK(akLv1,bigK_Lv1,mL,irows,icols,iUi)
    CALL getbigK(akLv2,bigK_Lv2,mL,irows,icols,iUi)
    CALL getbigK(akLm1,bigK_Lm1,mL,irows,icols,iUi)
    CALL getbigK(akLp1,bigK_Lp1,mL,irows,icols,iUi)
    irow_L = irows
#if _DEBUG_ == 2
    PRINT *, "Here debug compK_Lv"
    DO I = 1, irows
       DO J = 1, icols
          WRITE(*, "(F8.2)", ADVANCE='NO') bigK_Lv1(I,J)
                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLv2(I,J)
                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLm1(I,J)
                                !           WRITE(*, "(F8.2)", ADVANCE='NO') akLp1(I,J)
       END DO
       WRITE(*,*)
    END DO
    return
#endif

    !! Define bigU_L
    ALLOCATE(bigU_L(irow_L))
    bigU_L = 0.0d0
    DO i = 1,irow_L
       DO j = 1,mL
          IF (iu_L(j).EQ.1) THEN
             bigU_L(i) = bigU_L(i) + dfloat(iUi(i,j))
          ENDIF
       ENDDO
    ENDDO
    !! free memory
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    !! RyR channel
    iU_R(1) = 0
    iU_R(2) = 1
    !! Call ballbin to get possible state configurations
    balls = N_R
    bins = mR
    CALL count_ball2bin(balls, bins, iUi, irows, icols)

    !! rate transition
    akR12 = 3d0!2d5
    akR21 = 400d0
                                ! zeros matrices
    akRm2 = 0.0d0
    akRp2 = 0.0d0
                                ! NOTE:
                                !  + With RyR, we don't have voltage-dependent/independent, i.e. akRv1, akRv2    
                                !  + With RyR, we have (c_ds+c_jsr)-dependent/independent
    akRm2(2,1) = akR21
    akRm2(2,2) = -akR21

    akRp2(1,2) = akR12
    akRp2(1,1) = -akR12

    !! For RyR cluster
    CALL getbigK(akRm2,bigK_Rm2,mR,irows,icols,iUi)
    CALL getbigK(akRp2,bigK_Rp2,mR,irows,icols,iUi)
    irow_R = irows

    !! Define bigU_R
    ALLOCATE(bigU_R(irow_R))
    bigU_R = 0.0
    DO i = 1,irow_R
       DO j = 1,mR
          IF (iu_R(j).EQ.1) THEN
             bigU_R(i) = bigU_R(i) + dfloat(iUi(i,j))
          ENDIF
       ENDDO
    ENDDO
    !! free memory
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)


    !! Set up identity matrix
    ALLOCATE(aI_L(irow_L,irow_L))
    aI_L = 0.0
    DO i = 1,irow_L
       aI_L(i,i) = 1.0
    ENDDO
    ALLOCATE(aI_R(irow_R,irow_R))
    aI_R = 0.0
    DO i = 1,irow_R
       aI_R(i,i) = 1.0
    ENDDO
    !! Kronecker things
    ALLOCATE(aKv1(irow_L*irow_R,irow_L*irow_R))
    ALLOCATE(aKv2(irow_L*irow_R,irow_L*irow_R))
    CALL kron(bigK_Lv1,aI_R,aKv1,irow_L,irow_L,irow_R,irow_R)
    CALL kron(bigK_Lv2,aI_R,aKv2,irow_L,irow_L,irow_R,irow_R)
    IF (ALLOCATED(bigK_Lv1)) DEALLOCATE(bigK_Lv1)
    IF (ALLOCATED(bigK_Lv2)) DEALLOCATE(bigK_Lv2)

    ALLOCATE(aKm1(irow_L*irow_R,irow_L*irow_R))
    ALLOCATE(aKm2(irow_L*irow_R,irow_L*irow_R))
    CALL kron(bigK_Lm1,aI_R,aKm1,irow_L,irow_L,irow_R,irow_R)
    CALL kron(aI_L,bigK_Rm2,aKm2,irow_L,irow_L,irow_R,irow_R)
    IF (ALLOCATED(bigK_Lm1)) DEALLOCATE(bigK_Lm1)
    IF (ALLOCATED(bigK_Rm2)) DEALLOCATE(bigK_Rm2)

    ALLOCATE(aKp1(irow_L*irow_R,irow_L*irow_R))
    ALLOCATE(aKp2(irow_L*irow_R,irow_L*irow_R))
    CALL kron(bigK_Lp1,aI_R,aKp1,irow_L,irow_L,irow_R,irow_R)
    CALL kron(aI_L,bigK_Rp2,aKp2,irow_L,irow_L,irow_R,irow_R)
    IF (ALLOCATED(bigK_Lp1)) DEALLOCATE(bigK_Lp1)
    IF (ALLOCATED(bigK_Rp2)) DEALLOCATE(bigK_Rp2)

    IF (ALLOCATED(aI_L)) DEALLOCATE(aI_L)
    IF (ALLOCATED(aI_R)) DEALLOCATE(aI_R)


    ALLOCATE(e_L(1,irow_L))
    e_L = 1.0
    ALLOCATE(U_Ro(1,irow_L*irow_R))
    CALL kron(e_L,bigU_R,U_Ro,1,irow_L,1,irow_R)
    IF (ALLOCATED(e_L)) DEALLOCATE(e_L)

    ALLOCATE(e_R(1,irow_R))
    e_R = 1.0
    ALLOCATE(U_Lo(1,irow_L*irow_R))
    CALL kron(bigU_L,e_R,U_Lo,1,irow_L,1,irow_R)
    IF (ALLOCATED(e_R)) DEALLOCATE(e_R)

    ALLOCATE(eye(irow_L*irow_R,irow_L*irow_R))
    eye = 0.0
    DO i = 1,irow_L*irow_R
       eye(i,i) = 1.0
    ENDDO

    RETURN
  END SUBROUTINE makeKs

  !==============================================
  ! version 4 - work with full ECC model
  ! return
  !   compK, idxK
  !   maxnkl
  !   iBound
  !==============================================      
  SUBROUTINE getcompK_4(ak, m, iUi, irows, compK, idxK, maxnkl, N_chan)
    IMPLICIT NONE

    REAL(KIND=dp), DIMENSION(m,m) :: ak ! rate transition from one channel state to another
    INTEGER, DIMENSION(irows,m) :: iUi 
    INTEGER, INTENT(IN) :: m, & ! # of channel state
         irows, &  ! # of cluster states
         N_chan    ! # of channels
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK ! matrix of rate transition from one cluster-state
                                ! to another (keep only non-zero values)
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  idxK          ! keep index of the state
    INTEGER, INTENT(OUT) :: maxnkl
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    REAL(KIND=dp), DIMENSION(m) :: vcount
    INTEGER :: icount, ii, jj, initial_state, final_state
    REAL(kind=dp) :: asum

    !  Find maximum # of possible neighbors
    !+ maximum number of state transition is the # non-zero non-diagonal elements
    icount = +0
!!$    IF (N_chan .GE. m) THEN
!!$       DO ii = 1, m 
!!$          mask =  ak(ii,:) .GT. EPS
!!$          icount = icount + COUNT(mask)
!!$       ENDDO
!!$    ELSE
    !! The above value work for N_chan > m
    !! If N_chan < m, a better boundary would be
    !!+ count the number non-zero non-diagonal elements from
    !!+ N_chan lines of maximum number of non-zero values
!!$       DO ii = 1, m 
!!$          mask =  ak(ii,:) .GT. EPS
!!$          vcount(ii) = -COUNT(mask)
!!$       ENDDO
!!$       CALL inspar(vcount, N_chan)
!!$       icount = SUM(-vcount(1:N_chan))
!!$    END IF
    DO ii = 1, m 
       mask =  ak(ii,:) .GT. EPS
       vcount(ii) = -COUNT(mask)
    ENDDO
    CALL inspar(vcount, MIN(N_chan,m))
    icount = SUM(-vcount(1:MIN(N_chan,m)))
#if _DEBUG_ == 3
    WRITE(*,*) icount
#endif
    !! Finally
    maxnkl = icount + 1     ! +1 column to keep state unchanged
#if _DEBUG_ == 3
    WRITE(*,*) "maxnkl", maxnkl
#endif
    
    !!     Gather information from big K matrices
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
                                ! the zero-th column is for consistency with idxK matrix
    ALLOCATE(compK(irows, 0:maxnkl))

    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
    ! As different state has different number of state transition
    !+ the zero-th column keep the value indicating the true # possible state transition,
    !+ excluding the current state
    ALLOCATE(idxK(irows, 0:maxnkl))
    !! default values
    compK = +0.d0
    idxK = +0

    ! check all pairs of cluster state, see which pair is valid for a neighbor in state transition
    DO ii = 1, irows
       icount = 1
       asum = +0.d0
       DO jj = 1,irows
          tmpvec = iUi(ii,1:m) - iUi(jj,1:m)
          ! valid for state change = yes
          !+ when only one element = 1, one element = -1, the others = 0
          IF (SUM(ABS(tmpvec)) .EQ. 2) THEN 
             initial_state = MAXLOC(tmpvec, 1)  ! location of element with value 1
             final_state = MINLOC(tmpvec,1)     ! location of element with value -1
             IF (ak(initial_state, final_state) .GT. EPS) THEN ! positive rate transition 
                icount = icount + 1
                compK(ii,icount) = dfloat(iUi(ii,initial_state)) &
                     *ak(initial_state,final_state)
                idxK(ii,icount) = jj
                asum = asum + compK(ii, icount)
             ENDIF
          ENDIF
       ENDDO
       idxK(ii, 1) = ii
       ! Find the negative value at each row of aK
       !+ and keep it at the first column
       !+ one minus to this value would the probability that the state doesn't change
       IF (asum .GT. EPS) THEN
          compK(ii, 1) = -asum
       ENDIF
       ! keep the true number of neighbors at state ii (including itself)
       IF (icount > 1) THEN
          idxK(ii, 0) = icount
       ELSE
          idxK(ii, 0) = +0
       ENDIF
    ENDDO
    
    RETURN
  END SUBROUTINE getcompK_4
  !! Transpose version
  SUBROUTINE getcompK_4_T(ak, m, iUi, irows, compK, idxK, maxnkl, N_chan)
    IMPLICIT NONE

    REAL(KIND=dp), DIMENSION(m,m) :: ak ! rate transition from one channel state to another
    INTEGER, DIMENSION(irows,m) :: iUi 
    INTEGER, INTENT(IN) :: m, & ! # of channel state
         irows, &  ! # of cluster states
         N_chan    ! # of channels
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: compK ! matrix of rate transition from one cluster-state
                                ! to another (keep only non-zero values)
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  idxK          ! keep index of the state
    INTEGER, INTENT(OUT) :: maxnkl
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    REAL(KIND=dp), DIMENSION(m) :: vcount
    INTEGER :: icount, ii, jj, initial_state, final_state
    REAL(kind=dp) :: asum

                                !  Find maximum # of possible neighbors
                                !+ maximum number of state transition is the # non-zero non-diagonal elements
    icount = +0
!!$    IF (N_chan .GE. m) THEN
!!$       DO ii = 1, m 
!!$          mask =  ak(ii,:) .GT. EPS
!!$          icount = icount + COUNT(mask)
!!$       ENDDO
!!$    ELSE
    !! The above value work for N_chan > m
    !! If N_chan < m, a better boundary would be
    !!+ count the number non-zero non-diagonal elements from
    !!+ N_chan lines of maximum number of non-zero values
!!$       DO ii = 1, m 
!!$          mask =  ak(ii,:) .GT. EPS
!!$          vcount(ii) = -COUNT(mask)
!!$       ENDDO
!!$       CALL inspar(vcount, N_chan)
!!$       icount = SUM(-vcount(1:N_chan))
!!$    END IF
    DO ii = 1, m 
       mask =  ak(ii,:) .GT. EPS
       vcount(ii) = -COUNT(mask)
    ENDDO
    CALL inspar(vcount, MIN(N_chan,m))
    icount = SUM(-vcount(1:MIN(N_chan,m)))
#if _DEBUG_ == 3
    WRITE(*,*) icount
#endif
    !! Finally
    maxnkl = icount + 1     ! +1 column to keep state unchanged
#if _DEBUG_ == 3
    WRITE(*,*) "maxnkl", maxnkl
#endif

    !!     Gather information from big K matrices
    IF (ALLOCATED(compK)) DEALLOCATE(compK)
                                ! the zero-th row is for consistency with idxK matrix
    ALLOCATE(compK(0:maxnkl, irows))

    IF (ALLOCATED(idxK)) DEALLOCATE(idxK)
    ! As different state has different number of state transition
    !+ the zero-th column keep the value indicating the true # possible state transition,
    !+ excluding the current state
    ALLOCATE(idxK(0:maxnkl, irows))

    !! default values
    compK = +0.d0
    idxK = +0

    ! check all pairs of cluster state, see which pair is valid for a neighbor in state transition
    DO ii = 1, irows
       icount = 1
       asum = +0.d0
       DO jj = 1,irows
          tmpvec = iUi(ii,1:m) - iUi(jj,1:m)
          ! valid for state change = yes
          !+ when only one element = 1, one element = -1, the others = 0
          IF (SUM(ABS(tmpvec)) .EQ. 2) THEN 
             initial_state = MAXLOC(tmpvec, 1)  ! location of element with value 1
             final_state = MINLOC(tmpvec,1)     ! location of element with value -1
             IF (ak(initial_state, final_state) .GT. EPS) THEN ! positive rate transition 
                icount = icount + 1
                compK(icount, ii) = dfloat(iUi(ii,initial_state)) &
                     *ak(initial_state,final_state)
                idxK(icount, ii) = jj
                asum = asum + compK(icount, ii)
             ENDIF
          ENDIF
       ENDDO
       idxK(1, ii) = ii
       ! Find the negative value at each row of aK
       !+ and keep it at the first column
       !+ one minus to this value would the probability that the state doesn't change
       IF (asum .GT. EPS) THEN
          compK(1, ii) = -asum
       ENDIF
       ! keep the true number of neighbors at state ii (including itself)
       IF (icount > 1) THEN
          idxK(0, ii) = icount
       ELSE
          idxK(0, ii) = +0
       ENDIF
    ENDDO
    RETURN
  END SUBROUTINE getcompK_4_T

  !==============================================
  ! return
  ! compK, indexK
  ! maxnkl
  ! NOTE: This version works with 2 state only
  !       A general version should be used, defined as getcompK2
  !==============================================      
  SUBROUTINE getcompK(ak,compK, indexK,m,irows,maxnkl,iUi)
    IMPLICIT NONE
    REAL(KIND=dp), DIMENSION(m,m):: ak
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  indexK
    INTEGER, INTENT(IN) :: m, irows
    INTEGER, INTENT(OUT) :: maxnkl
    INTEGER, DIMENSION(irows,m) :: iUi  
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    INTEGER :: ii, jj, initial_state, ifinal_state, icount, istate_change
    REAL(KIND=dp) :: asum

    maxnkl = m*(m-1)
                                !     Gather information from big K matrices
    if (.not. allocated(compK)) then
       allocate(compK(irows, 0:maxnkl))
    endif
    if (.not. allocated(indexK)) then
       allocate(indexK(irows, maxnkl))
    endif

                                ! default values       
    compK = 0d0
    indexK = 0
    istate_change = 0

                                ! check all pairs of cluster state
    do ii = 1,irows
       icount = 0
       asum = 0d0
       do jj = 1,irows
          tmpvec = iUI(ii,:) - iUI(jj,:)
          if (SUM(ABS(tmpvec)) .eq. 2) then
             initial_state = maxloc(tmpvec, 1)
             ifinal_state = minloc(tmpvec,1)
             if (abs(ak(initial_state, ifinal_state)) .gt. 0d0) then
                icount = icount + 1
                compK(ii,icount) = dfloat(iUI(ii,initial_state)) &
                     *ak(initial_state,ifinal_state)
                indexK(ii,icount) = jj
                asum = asum + compK(ii, icount)
             endif
          endif
       enddo
                                ! find the negative value at each row of aK
       if (asum .gt. +0d0) then
!!$          icount = icount + 1
!!$          compK(ii, icount) = -asum
!!$          indexK(ii, icount) = ii
          compK(ii, maxnkl) = -asum
          indexK(ii, maxnkl) = ii
       endif
    enddo
    return
  END SUBROUTINE getcompK

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                ! Work with Leak model
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE getcompK2(ak,compK, indexK,m,irows,maxnkl,iUi)
                                !    IMPLICIT REAL(KIND=dp) (a-h,o-z)
    IMPLICIT NONE
    REAL(KIND=dp), DIMENSION(m,m):: ak
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  indexK
    INTEGER, INTENT(IN) :: m, irows
    INTEGER, INTENT(OUT) :: maxnkl
    INTEGER, DIMENSION(irows,m) :: iUi  
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    INTEGER :: ii, jj, icount, istate_change, initial_state, ifinal_state
    REAL(KIND=dp) :: asum


                                !old formula:    maxnkl = m*(m-1)
                                !better with:
                                !find number of rows with non-zero elements
    icount = 0
    DO ii = 1, m
       DO jj=1,m
          IF (ak(ii,jj) .GT. EPS) THEN
             icount = icount + 1
          ENDIF
       ENDDO
    ENDDO
    maxnkl = icount + 1

                                !     Gather information from big K matrices
    if (.not. allocated(compK)) then
       allocate(compK(irows, 0:maxnkl))
    endif
    if (.not. allocated(indexK)) then
       allocate(indexK(irows, maxnkl))
    endif

                                ! default values       
    compK = 0d0
    indexK = 0
    istate_change = 0

                                ! check all pairs of cluster state
    do ii = 1,irows
       icount = 0
       asum = 0d0
       do jj = 1,irows
          tmpvec = iUI(ii,:) - iUI(jj,:)
          if (SUM(ABS(tmpvec)) .eq. 2) then
             initial_state = maxloc(tmpvec, 1)
             ifinal_state = minloc(tmpvec,1)
             IF (ABS(ak(initial_state, ifinal_state)) .GT. EPS) THEN
                icount = icount + 1
                compK(ii,icount) = dfloat(iUI(ii,initial_state)) &
                     *ak(initial_state,ifinal_state)
                indexK(ii,icount) = jj
                asum = asum + compK(ii, icount)
             endif
          endif
       enddo
                                ! find the negative value at each row of aK
       if (asum .gt. EPS) then
          compK(ii, maxnkl) = -asum
          indexK(ii, maxnkl) = ii
       endif
    enddo
    return
  END SUBROUTINE getcompK2
                                ! The transpose version
                                ! NOTE: Beside transposing the matrix, another change compare to the non-transpose version
                                !  is that
                                !       the value to keep the negative value is now at row 1, instead of the last row
  SUBROUTINE getcompK2_T(ak,compK, indexK,m,irows,maxnkl,iUi)
    IMPLICIT NONE
    REAL(KIND=dp), DIMENSION(m,m):: ak
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compK
    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  indexK
    INTEGER, INTENT(IN) :: m, irows
    INTEGER, INTENT(OUT) :: maxnkl
    INTEGER, DIMENSION(irows,m) :: iUi  
                                !local
    INTEGER, DIMENSION(m) :: tmpvec
    LOGICAL, DIMENSION(m) :: mask
    INTEGER :: ii, jj, icount, istate_change, initial_state, ifinal_state
    REAL(KIND=dp) :: asum


                                !old formula:    maxnkl = m*(m-1)
                                !better with:
                                !find number of rows with non-zero elements
    icount = 0
    DO ii = 1, m
       DO jj=1,m
          IF (ak(ii,jj) .GT. EPS) THEN
             icount = icount + 1
          ENDIF
       ENDDO
    ENDDO
    maxnkl = icount + 1

                                !     Gather information from big K matrices
    if (.not. allocated(compK)) then
       ALLOCATE(compK(0:maxnkl, irows))
    endif
    if (.not. allocated(indexK)) then
       ALLOCATE(indexK(maxnkl, irows))
    endif

                                ! default values       
    compK = 0d0
    indexK = 0
    istate_change = 0

                                ! check all pairs of cluster state
    do ii = 1,irows
       icount = 1
       asum = 0d0
       do jj = 1,irows
          tmpvec = iUI(ii,:) - iUI(jj,:)
          if (SUM(ABS(tmpvec)) .eq. 2) then
             initial_state = maxloc(tmpvec, 1)
             ifinal_state = minloc(tmpvec,1)
             IF (ABS(ak(initial_state, ifinal_state)) .GT. EPS) THEN
                icount = icount + 1
                compK(icount, ii) = dfloat(iUI(ii,initial_state)) &
                     *ak(initial_state,ifinal_state)
                indexK(icount, ii) = jj
                asum = asum + compK(icount, ii)
             endif
          endif
       enddo
                                ! find the negative value at each row of aK
       if (asum .gt. EPS) then
          compK(1, ii) = -asum
          indexK(1, ii) = ii
       endif
    enddo
    return
  END SUBROUTINE getcompK2_T



  SUBROUTINE export_Umat(numChan, numState, filename)
    INTEGER :: numChan, numState
    CHARACTER(*), OPTIONAL:: filename
    !local
    INTEGER :: balls, bins
    INTEGER :: irows, icols
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    
    balls=numChan
    bins=numState

    CALL count_ball2bin(balls, bins, iUi, irows, icols)
    ! NOTE: icols == bins
    ! iUi - matrix with each row denote a L-type group's state
    !       a row contains the number of L-type channels
    !     + in each L-type state

    WRITE(*,*) "balls, bins: ", balls, bins
    WRITE(*,*) "irows, icols: ", irows, icols
    IF (PRESENT(filename)) THEN
       CALL print_mat(iUi, irows, icols, filename)
    ELSE
       CALL print_mat(iUi, irows, icols)
    END IF
    STOP

  END SUBROUTINE export_Umat

  !!FUNCTION read_RyR_Markov_4(....) RESULT(res)
  FUNCTION read_RyR_Markov_4(filename, akRp_mode1, akRm_mode1, akRp_mode2, akRm_mode2, akRm, &
                             akRos, mR, iU_R, initial_state) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    CHARACTER(FILENAME_LENGTH) :: filename  ! name of the configuration file
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: akRm_mode1, akRp_mode1, akRm_mode2, akRp_mode2, akRm, akRos
    INTEGER, INTENT(OUT) :: mR
    INTEGER, DIMENSION(:), ALLOCATABLE :: iU_R
    INTEGER, INTENT(OUT) :: initial_state
    !local
    INTEGER :: fileID, error_code
    LOGICAL :: not_found
    CHARACTER :: dummy_char
    CHARACTER(40) :: line
    REAL(KIND=dp) :: rowsum_p_mode1, rowsum_m_mode1,rowsum_p_mode2, rowsum_m_mode2, rowsum_m, rowsum_Ros
    INTEGER :: ii, jj

    res = SUCCESS
    error_code = 0
    
    fileID = get_file_ID()
    OPEN(unit=fileID, file=filename, ACTION='READ')

    !! Read-in the number of state mR
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) mR
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !! Read-in the open states
    ! set-up which state is open
    not_found = .TRUE.
    ALLOCATE(iU_R(mR))
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,*) iU_R
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !! Read-in which state should be the initial state
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) initial_state
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    ALLOCATE(akRp_mode1(mR,mR), akRp_mode2(mR,mR), akRos(mR,mR))    
    ALLOCATE(akRm_mode1(mR,mR), akRm_mode2(mR,mR), akRm(mR,mR))    
    ! zero-matrices
    akRp_mode1 = +0.0d0
    akRm_mode1 = +0.0d0
    akRp_mode2 = +0.0d0
    akRm_mode2 = +0.0d0
    akRm = +0.0d0
    akRos = +0.0d0
    
    !! Read-in the transition-rates
    DO WHILE (error_code == 0)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) line
          SELECT CASE (TRIM(line)) 
          CASE ('!Ca_act_mode1')
             !PRINT *, 'Caact'
             CALL read_transition(akRp_mode1, fileID)
          CASE ('!Ca_act_mode2')
             CALL read_transition(akRp_mode2, fileID)
          CASE ('!free_mode1')
             CALL read_transition(akRm_mode1, fileID)
          CASE ('!free_mode2')
             CALL read_transition(akRm_mode2, fileID)
          CASE ('!free')
             !PRINT *, 'free'
             CALL read_transition(akRm, fileID)
          CASE ('!Ros_dep')
             !PRINT *, 'free'
             CALL read_transition(akRos, fileID)
          END SELECT
       END IF
    END DO

    CLOSE(unit = fileID)

    DO ii = 1,mR
       rowsum_p_mode1 = +0.d0
       rowsum_p_mode2 = +0.d0
       rowsum_m_mode1 = +0.d0
       rowsum_m_mode2 = +0.d0
       rowsum_m = +0.d0
       rowsum_Ros = +0.d0
       DO jj = 1,mR
          rowsum_p_mode1 = rowsum_p_mode1 + akRp_mode1(ii,jj)
          rowsum_p_mode2 = rowsum_p_mode2 + akRp_mode2(ii,jj)
          rowsum_m_mode1 = rowsum_m_mode1 + akRm_mode1(ii,jj)
          rowsum_m_mode2 = rowsum_m_mode2 + akRm_mode2(ii,jj)          
          rowsum_m = rowsum_m + akRm(ii,jj)
          rowsum_Ros = rowsum_Ros + akRos(ii,jj)
       ENDDO
       akRp_mode1(ii,ii) = -rowsum_p_mode1
       akRp_mode2(ii,ii) = -rowsum_p_mode2
       akRm_mode1(ii,ii) = -rowsum_m_mode1
       akRm_mode2(ii,ii) = -rowsum_m_mode2
       akRm(ii,ii) = -rowsum_m
       akRos(ii,ii) = -rowsum_Ros
    ENDDO


  END FUNCTION read_RyR_Markov_4

  FUNCTION read_RyR_Markov_2(filename, akRp, akRm, mR, iU_R, initial_state) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    CHARACTER(FILENAME_LENGTH) :: filename  ! name of the configuration file
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: akRm, akRp
    INTEGER, INTENT(OUT) :: mR
    INTEGER, DIMENSION(:), ALLOCATABLE :: iU_R
    INTEGER, INTENT(OUT) :: initial_state
    !local
    INTEGER :: fileID, error_code
    LOGICAL :: not_found
    CHARACTER :: dummy_char
    CHARACTER(40) :: line
    REAL(KIND=dp) :: rowsum_m, rowsum_p
    INTEGER :: ii, jj

    res = SUCCESS
    error_code = 0
    
    fileID = get_file_ID()
    OPEN(unit=fileID, file=filename, ACTION='READ')

    !! Read-in the number of state mR
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) mR
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !! Read-in the open states
    ! set-up which state is open
    not_found = .TRUE.
    ALLOCATE(iU_R(mR))
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,*) iU_R
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !! Read-in which state should be the initial state
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) initial_state
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    ALLOCATE(akRm(mR,mR), akRp(mR,mR))    
    ! zero-matrices
    akRm = +0.0d0
    akRp = +0.0d0
    !! Read-in the transition-rates
    DO WHILE (error_code == 0)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) line
          SELECT CASE (TRIM(line)) 
          CASE ('!Ca_act')
             !PRINT *, 'Caact'
             CALL read_transition(akRp, fileID)
          CASE ('!free')
             !PRINT *, 'free'
             CALL read_transition(akRm, fileID)
          END SELECT
       END IF
    END DO

    CLOSE(unit = fileID)

    DO ii = 1,mR
       rowsum_m = +0.d0
       rowsum_p = +0.d0
       DO jj = 1,mR
          rowsum_m = rowsum_m + akRm(ii,jj)
          rowsum_p = rowsum_p + akRp(ii,jj)
       ENDDO
       akRm(ii,ii) = -rowsum_m
       akRp(ii,ii) = -rowsum_p
    ENDDO
    
  END FUNCTION  read_RyR_Markov_2


  FUNCTION read_LCC_Markov(filename, akLv1, akLv2, akLp1, akLm1, mL, iU_L, initial_state) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: akLv1, akLv2, akLm1, akLp1
    INTEGER, INTENT(OUT) :: mL
    CHARACTER(FILENAME_LENGTH) :: filename  ! name of the configuration file
    INTEGER, DIMENSION(:), ALLOCATABLE :: iU_L
    INTEGER, INTENT(OUT) :: initial_state
    !local
    INTEGER :: fileID, error_code
    LOGICAL :: not_found
    CHARACTER :: dummy_char
    CHARACTER(40) :: line
    REAL(KIND=dp) :: rowsum_v1, rowsum_v2, rowsum_m1, rowsum_p1
    INTEGER :: ii, jj

    res = SUCCESS
    error_code = 0
    
    fileID = get_file_ID()
    OPEN(unit=fileID, file=filename, ACTION='READ')
    !WRITE(strFormat, '("(",I2,"(I2,1x))")') mL
    !READ(fileID,strFormat, IOSTAT=error_code) LCC_states

    !! Read-in the number of state mL
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) mL
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !WRITE(strFormat, '("(",I2,"(I2,1x))")') mL

    !! Read-in the open states
    ! set-up which state is open
    not_found = .TRUE.
    ALLOCATE(iU_L(mL))
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,*) iU_L
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    !! Read-in which state should be the initial state
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) initial_state
          not_found = .FALSE.
       END IF
    END DO
    IF (error_code .NE. 0) THEN 
       res = ERR_FILE_READING_NOT_COMPLETED
       WRITE(error_message, *) "Wrong Markov model file-format"
       RETURN 
    ENDIF

    ALLOCATE(akLv1(mL,mL),akLv2(mL,mL),akLm1(mL,mL),akLp1(mL,mL))
    ! zero-matrices
    akLv1 = +0.0d0
    akLv2 = +0.0d0
    akLm1 = +0.0d0
    akLp1 = +0.0d0
    !! Read-in the transition-rates
    DO WHILE (error_code == 0)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID, *) line
          SELECT CASE (TRIM(line)) 
          CASE ('!Vm_act')
             CALL read_transition(akLv1, fileID)
          CASE ('!Vm_inact')
             CALL read_transition(akLv2, fileID)
          CASE ('!Ca_inact')
             CALL read_transition(akLp1, fileID)
          CASE ('!free')
             CALL read_transition(akLm1, fileID)
          END SELECT
       END IF
    END DO

    CLOSE(unit = fileID)

    DO ii = 1,mL
       rowsum_v1 = +0.d0
       rowsum_v2 = +0.d0
       rowsum_m1 = +0.d0
       rowsum_p1 = +0.d0
       DO jj = 1,mL
          rowsum_v1 = rowsum_v1 + akLv1(ii,jj)
          rowsum_v2 = rowsum_v2 + akLv2(ii,jj)
          rowsum_m1 = rowsum_m1 + akLm1(ii,jj)
          rowsum_p1 = rowsum_p1 + akLp1(ii,jj)
       ENDDO
       akLv1(ii,ii) = -rowsum_v1
       akLv2(ii,ii) = -rowsum_v2
       akLm1(ii,ii) = -rowsum_m1
       akLp1(ii,ii) = -rowsum_p1
    ENDDO
  END FUNCTION read_LCC_Markov
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Deal with each homogeneous cluster separately
  !! These functions are expected to work with any LCC/RyR models
  !! as long as LCC has v1,v2,m1, p1
  !!            RyR has m2,p2
  !! IMPORTANT: We use non-transpose form as default, i.e. compK(irow,maxnkp)
  !!  mR = number of RyR channels in the cluster
  !!  U_Ro = vector that keep # of open RyRs in each homogeneous cluster state
  !!  N_R = total number of RyR
  !!  Ncol_... = number of column... in compact form matrices
  !!  compK... = compact form matrices
  !!  indxK... = index ...
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  FUNCTION setup_RyR_2(akRp, akRm, N_R, mR, iU_R, compKcal_act, compKfree, indxKcal_act, &
                indxKfree, U_Ro, Ncol_Kcal_act, Ncol_Kfree, irow_R) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:) :: akRp, akRm
    INTEGER, INTENT(IN) :: N_R !# RyRs
    INTEGER, INTENT(IN) :: mR !# of state in Markov-model of single channel
    INTEGER, DIMENSION(:) :: iU_R
    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKcal_act, compKfree 
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: indxKcal_act, indxKfree
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro
    INTEGER, INTENT(OUT) :: Ncol_Kcal_act, Ncol_Kfree
    INTEGER, INTENT(OUT) :: irow_R !# of RyR cluster-state
		
    !local
    INTEGER :: balls, bins
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    INTEGER :: icols, ii, jj

    
    res = SUCCESS

    balls=N_R
    bins=mR

    CALL count_ball2bin(balls, bins, iUi, irow_R, icols)

    !size (mR,mR): akRm, akRp
    !size (irow_R,mR): iUi
    !size (irow_R,0:Ncol_...) :: compK_R.
    !size (irow_R,0:Ncol_...) :: indxK_R.
    CALL getcompK_4(akRm, mR, iUi, irow_R, compKfree, indxKfree, Ncol_Kfree, N_R)
    CALL getcompK_4(akRp, mR, iUi, irow_R, compKcal_act, indxKcal_act, Ncol_Kcal_act, N_R)

    
    ! Define bigU_R
    ! + vector that tells how many active RyRs in each RyR cluster-state
    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(irow_R))

    U_Ro = +0.d0
    DO ii = 1, irow_R
       DO jj = 1,mR
          IF (iU_R(jj).EQ.1) THEN
             U_Ro(ii) = U_Ro(ii) + dfloat(iUi(ii,jj))
          ENDIF
       ENDDO
    ENDDO

    !free
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    PRINT *, "+end setup_RyR()"
    RETURN  
  END FUNCTION   setup_RyR_2
  FUNCTION setup_RyR_4(akRp_mode1, akRm_mode1,akRp_mode2,akRm_mode2,akRm,akRos, N_R, mR, iU_R, compKcal_act_mode1, &
                compKfree_mode1, compKcal_act_mode2,compKfree_mode2, compKfree, compKRos, indxKcal_act_mode1, &
                indxKfree_mode1, indxKcal_act_mode2, indxKfree_mode2, indxKfree, indxKRos, U_Ro, U_RyR_State, &
                Ncol_Kcal_act_mode1, Ncol_Kfree_mode1, Ncol_Kcal_act_mode2, Ncol_Kfree_mode2, &
                Ncol_Kfree, Ncol_KRos, irow_R) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:) :: akRp_mode1, akRm_mode1, akRp_mode2, akRm_mode2, akRm, akRos
    INTEGER, INTENT(IN) :: N_R !# RyRs
    INTEGER, INTENT(IN) :: mR !# of state in Markov-model of single channel
    INTEGER, DIMENSION(:) :: iU_R
    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKcal_act_mode1, compKfree_mode1, compKcal_act_mode2, &
                                                  compKfree_mode2, compKfree, compKRos 
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: indxKcal_act_mode1, indxKfree_mode1, indxKcal_act_mode2, indxKfree_mode2, &
                                            indxKfree, indxKRos
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: U_RyR_State
    INTEGER, INTENT(OUT) :: Ncol_Kcal_act_mode1, Ncol_Kfree_mode1, Ncol_Kcal_act_mode2, Ncol_Kfree_mode2, Ncol_Kfree, Ncol_KRos
    INTEGER, INTENT(OUT) :: irow_R !# of RyR cluster-state
		
    !local
    INTEGER :: balls, bins
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    INTEGER :: icols, ii, jj

    
    res = SUCCESS

    balls=N_R
    bins=mR
    CALL count_ball2bin(balls, bins, iUi, irow_R, icols)

    !size (mR,mR): akRm, akRp
    !size (irow_R,mR): iUi
    !size (irow_R,0:Ncol_...) :: compK_R.
    !size (irow_R,0:Ncol_...) :: indxK_R.
    CALL getcompK_4(akRm_mode1, mR, iUi, irow_R, compKfree_mode1, indxKfree_mode1, Ncol_Kfree_mode1, N_R)
    PRINT *," after akRm_mode1 getcompK_4()" 
    CALL getcompK_4(akRm_mode2, mR, iUi, irow_R, compKfree_mode2, indxKfree_mode2, Ncol_Kfree_mode2, N_R)
    PRINT *," after akRm_mode2 getcompK_4()" 
    CALL getcompK_4(akRp_mode1, mR, iUi, irow_R, compKcal_act_mode1, indxKcal_act_mode1, Ncol_Kcal_act_mode1, N_R)
    PRINT *," after akRp_mode1 getcompK_4()" 
    CALL getcompK_4(akRp_mode2, mR, iUi, irow_R, compKcal_act_mode2, indxKcal_act_mode2, Ncol_Kcal_act_mode2, N_R)
    PRINT *," after akRp_mode2 getcompK_4()" 
    CALL getcompK_4(akRm, mR, iUi, irow_R, compKfree, indxKfree, Ncol_Kfree, N_R)
    PRINT *," after akRm getcompK_4()" 
    CALL getcompK_4(akRos, mR, iUi, irow_R, compKRos, indxKRos, Ncol_KRos, N_R)

    
    ! Define bigU_R
    ! + vector that tells how many active RyRs in each RyR cluster-state
    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(irow_R))
    IF (ALLOCATED(U_RyR_State)) DEALLOCATE(U_RyR_State)
    ALLOCATE(U_RyR_State(irow_R, mR))

    U_Ro = +0.d0
    U_RyR_State = +0.d0
    DO ii = 1, irow_R
       DO jj = 1,mR
          IF (iU_R(jj).EQ.1) THEN
             U_Ro(ii) = U_Ro(ii) + dfloat(iUi(ii,jj))
          ENDIF
          U_RyR_State(ii,jj) = U_RyR_State(ii,jj) + dfloat(iUi(ii,jj))
       ENDDO
    ENDDO

    !free
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    PRINT *, "+end setup_RyR()"
    RETURN  
  END FUNCTION   setup_RyR_4
  FUNCTION setup_RyR_T(akRp, akRm, N_R, mR, iU_R, compKcal_act, compKfree, &
       indxKcal_act, indxKfree, &
       U_Ro, Ncol_Kcal_act, Ncol_Kfree, irow_R) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:) :: akRp, akRm
    INTEGER, INTENT(IN) :: N_R !# RyRs
    INTEGER, INTENT(IN) :: mR !# of state in Markov-model of single channel
    INTEGER, DIMENSION(:) :: iU_R
    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKcal_act, compKfree
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: indxKcal_act, indxKfree
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Ro
    INTEGER, INTENT(OUT) :: Ncol_Kcal_act, Ncol_Kfree
    INTEGER, INTENT(OUT) :: irow_R !# of RyR cluster-state
		
    !local
    INTEGER :: balls, bins
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi

    INTEGER :: icols, ii, jj

    res = SUCCESS

    balls=N_R
    bins=mR

    CALL count_ball2bin(balls, bins, iUi, irow_R, icols)
    
    !size (mR,mR): akRm, akRp
    !size (irow_R,mR): iUi
    !size (0:Ncol_...,irow_R) :: compK_R.
    !size (0:Ncol_...,irow_R) :: indxK_R.
    CALL getcompK_4_T(akRm, mR, iUi, irow_R, compKfree, indxKfree, Ncol_Kfree, N_R)
    CALL getcompK_4_T(akRp, mR, iUi, irow_R, compKcal_act, indxKcal_act, Ncol_Kcal_act, N_R)

    ! Define bigU_R
    ! + vector that tells how many active RyRs in each RyR cluster-state
    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(irow_R))

    U_Ro = +0.d0
    DO ii = 1, irow_R
       DO jj = 1,mR
          IF (iU_R(jj).EQ.1) THEN
             U_Ro(ii) = U_Ro(ii) + dfloat(iUi(ii,jj))
          ENDIF
       ENDDO
    ENDDO

    !free
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    PRINT *, "+end setup_RyR_T()"
    RETURN
  
  END FUNCTION   setup_RyR_T
 
  !!!!!!!!!!!!!!!!!!!!!
  !! Vm-dependent activation/inactivation
  !! [Ca]ds-dependent inactivation
  !! Free parameters
  FUNCTION setup_LCC(akLv1, akLv2, akLp1, akLm1,  N_L, mL, iU_L, compKv_act, compKv_inact, compKcal_inact, compKfree, &
            indxKv_act, indxKv_inact, indxKcal_inact, indxKfree, &
            U_Lo, Ncol_Kv_act, Ncol_Kv_inact, Ncol_Kcal_inact, Ncol_Kfree, irow_L) RESULT(res)
    IMPLICIT NONE

    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:) :: akLv1, akLv2, akLm1, akLp1
    INTEGER, INTENT(IN) :: N_L !# LCCs
    INTEGER, INTENT(IN) :: mL !# of state in Markov-model of single channel
    INTEGER, DIMENSION(:) :: iU_L
    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv_act, compKv_inact, compKcal_inact, compKfree
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: indxKv_act, indxKv_inact, indxKcal_inact, indxKfree
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Lo
    INTEGER, INTENT(OUT) :: Ncol_Kv_act, Ncol_Kv_inact, Ncol_Kcal_inact, Ncol_Kfree
    INTEGER, INTENT(OUT) :: irow_L  !# of L-type cluster-state
    !local
    INTEGER :: balls, bins
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    INTEGER :: icols, ii, jj

    res = SUCCESS

    ! Call ballbin to get possible state configurations      
    balls=N_L
    bins=mL

    CALL count_ball2bin(balls, bins, iUi, irow_L, icols)
    !irow_L = irows ! 
    ! NOTE: icols == bins
    ! iUi - matrix with each row denote a L-type group's state
    !       a row contains the number of L-type channels
    !     + in each L-type state

    !size (m,m): akLv1, akLv2, akLm1, akLp1
    !size (irow_L,m): iUi
    !size (irow_L,0:maxnkL..) :: compK_L..
    !size (irow_L,0:maxnkL..) :: idxK_L..
    CALL getcompK_4(akLv1, mL, iUi, irow_L, compKv_act, indxKv_act, Ncol_Kv_act, N_L)
    CALL getcompK_4(akLv2, mL, iUi, irow_L, compKv_inact, indxKv_inact, Ncol_Kv_inact, N_L)
    CALL getcompK_4(akLm1, mL, iUi, irow_L, compKfree, indxKfree, Ncol_Kfree, N_L)
    CALL getcompK_4(akLp1, mL, iUi, irow_L, compKcal_inact, indxKcal_inact, Ncol_Kcal_inact, N_L)

    ! Define U_Lo
    ! + vector that tells how many active L-types in each L-type cluster-state
    IF (ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
    ALLOCATE(U_Lo(irow_L))

    DO ii = 1, irow_L
       U_Lo(ii) = +0.d0
       DO jj = 1, mL
          IF (iU_L(jj).EQ.1) THEN
             U_Lo(ii) = U_Lo(ii) + dfloat(iUi(ii,jj))
          ENDIF
       ENDDO
    ENDDO

    !free
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    PRINT *, "+end setup_LCC()"
    RETURN
  END FUNCTION  setup_LCC
  FUNCTION setup_LCC_T(akLv1, akLv2, akLp1, akLm1,  N_L, mL, iU_L, compKv_act, compKv_inact, compKcal_inact, compKfree, &
            indxKv_act, indxKv_inact, indxKcal_inact, indxKfree, &
            U_Lo, Ncol_Kv_act, Ncol_Kv_inact, Ncol_Kcal_inact, Ncol_Kfree, irow_L) RESULT(res)
    IMPLICIT NONE

    INTEGER :: res
    REAL(KIND=dp), DIMENSION(:,:) :: akLv1, akLv2, akLm1, akLp1
    INTEGER, INTENT(IN) :: N_L !# LCCs
    INTEGER, INTENT(IN) :: mL !# of state in Markov-model of single channel
    INTEGER, DIMENSION(:) :: iU_L
    !NOTE: these compK matrices have zero-th columns
    REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: compKv_act, compKv_inact, compKcal_inact, compKfree
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: indxKv_act, indxKv_inact, indxKcal_inact, indxKfree
    REAL(KIND=DP), DIMENSION(:), ALLOCATABLE :: U_Lo
    INTEGER, INTENT(OUT) :: Ncol_Kv_act, Ncol_Kv_inact, Ncol_Kcal_inact, Ncol_Kfree
    INTEGER, INTENT(OUT) :: irow_L  !# of L-type cluster-state
    !local
    INTEGER :: balls, bins
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    INTEGER :: icols, ii, jj
    
    res = SUCCESS
    ! Call ballbin to get possible state configurations      
    balls=N_L
    bins=mL

    CALL count_ball2bin(balls, bins, iUi, irow_L, icols)
    !irow_L = irows ! 
    ! NOTE: icols == bins
    ! iUi - matrix with each row denote a L-type group's state
    !       a row contains the number of L-type channels
    !     + in each L-type state

    !size (m,m): akLv1, akLv2, akLm1, akLp1
    !size (irow_L,m): iUi
    !size (0:maxnkL.., irow_L) :: compK_L..
    !size (0:maxnkL.., irow_L) :: idxK_L..
    CALL getcompK_4_T(akLv1, mL, iUi, irow_L, compKv_act, indxKv_act, Ncol_Kv_act, N_L)
    CALL getcompK_4_T(akLv2, mL, iUi, irow_L, compKv_inact, indxKv_inact, Ncol_Kv_inact, N_L)
    CALL getcompK_4_T(akLm1, mL, iUi, irow_L, compKfree, indxKfree, Ncol_Kfree, N_L)
    CALL getcompK_4_T(akLp1, mL, iUi, irow_L, compKcal_inact, indxKcal_inact, Ncol_Kcal_inact, N_L)

!!$    CALL print_mat(compKv_act, 0, Ncol_Kv_act+1, 1, 5)
!!$    stop

    ! Define U_Lo
    ! + vector that tells how many active L-types in each L-type cluster-state
    IF (ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
    ALLOCATE(U_Lo(irow_L))

    DO ii = 1, irow_L
       U_Lo(ii) = +0.d0
       DO jj = 1, mL
          IF (iU_L(jj).EQ.1) THEN
             U_Lo(ii) = U_Lo(ii) + dfloat(iUi(ii,jj))
          ENDIF
       ENDDO
    ENDDO
    !free
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    PRINT *, "+end setup_LCC_T()"
    RETURN
  END FUNCTION  setup_LCC_T
  
  SUBROUTINE read_transition(matA, fileID)
    !stop until '!' is seen or end of file
    IMPLICIT NONE
    INTEGER :: fileID
    REAL(KIND=dp), DIMENSION(:,:) :: matA

    !local
    INTEGER :: error_code
    LOGICAL :: not_found
    INTEGER :: ifrom, ito
    REAL(KIND=dp) :: rate
    CHARACTER :: dummy_char
    
    error_code = 0
    
    not_found = .TRUE.
    DO WHILE (error_code == 0 .AND. not_found)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       !PRINT *, "dumm", dummy_char
       IF (dummy_char .EQ. '!') THEN
          BACKSPACE(fileID)
          not_found = .FALSE.          
       ELSEIF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,*) ifrom, ito, rate
          PRINT *, ifrom, ito, rate
          matA(ifrom, ito) = rate
       END IF
    END DO
  END SUBROUTINE read_transition

  FUNCTION makeKs_nj(akRp2, akRm2, mR, iU_R, ryr_initial_state, &
       akLv1, akLv2, akLp1, akLm1, mL, iU_L, lcc_initial_state, &
       aKv1, aKv2, aKm1, aKp1, aKm2, aKp2, U_Lo, &
       U_Ro, pL_nj, pR_nj, irow_L, irow_R)  RESULT(res)
    IMPLICIT NONE
    INTEGER, PARAMETER :: N_L_nj = 1, N_R_nj = 1 ! assembled channel (RyR, LCC)
    INTEGER :: res
    !! arguments
    REAL(KIND=dp), DIMENSION(:,:) :: akLv1, akLv2, akLm1, akLp1
    REAL(KIND=dp), DIMENSION(:,:) :: akRm2, akRp2
    INTEGER:: mL  !mL = # of states for L-type channel
    INTEGER:: mR !mR = # of states for RyR channel
    INTEGER :: lcc_initial_state, ryr_initial_state
    INTEGER, DIMENSION(:) :: iU_L, iU_R
    !CHARACTER(FILENAME_LENGTH) :: lcc_filename, ryr_filename
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKv1
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKv2
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKm1
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKm2
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKp1
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKp2
    INTEGER, DIMENSION(:), ALLOCATABLE :: U_Lo
    INTEGER, DIMENSION(:), ALLOCATABLE :: U_Ro
    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pL_nj, pR_nj
    INTEGER, INTENT(OUT) :: irow_L, irow_R
    !! local
    INTEGER :: balls, bins, irows, icols
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: iUi
    !INTEGER :: ii, jj, fileID
    !CHARACTER :: dummy_char
    !INTEGER, DIMENSION(:), ALLOCATABLE :: iU_R, iU_L
    !INTEGER :: error_code
    !CHARACTER(40) :: line
    !LOGICAL :: not_found
    

    res = SUCCESS
    !! non-junctional RyR
    balls=N_R_nj
    bins=mR
    CALL count_ball2bin(balls, bins, iUi, irows, icols)
    irow_R = irows

    IF (ALLOCATED(U_Ro)) DEALLOCATE(U_Ro)
    ALLOCATE(U_Ro(irow_R))


    CALL getbigK(akRm2,aKm2,mR,irows,icols,iUi)
    CALL getbigK(akRp2,aKp2,mR,irows,icols,iUi)
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)

    !! non-junctional DHPR (lumped channel)
    balls = N_L_nj
    bins=mL
    CALL count_ball2bin(balls, bins, iUi, irows, icols)
    irow_L = irows ! # of states of assembled LCC

    IF (ALLOCATED(U_Lo)) DEALLOCATE(U_Lo)
    ALLOCATE(U_Lo(irow_L))
    U_Lo = iU_L


    CALL getbigK(akLv1,aKv1,mL,irows,icols,iUi)
    CALL getbigK(akLv2,aKv2,mL,irows,icols,iUi)
    CALL getbigK(akLm1,aKm1,mL,irows,icols,iUi)
    CALL getbigK(akLp1,aKp1,mL,irows,icols,iUi)
    !! free memory
    IF (ALLOCATED(iUi)) DEALLOCATE(iUi)
    
    !for non-junctional probabilities, i.e. fraction of channels in each state
    ALLOCATE(pL_nj(irow_L), pR_nj(irow_R))

    ! + set fraction of channel in each state
    pL_nj = 0d0
    pR_nj = 0d0
    pL_nj(lcc_initial_state) = 1d0  ! all LCC channel is in state 2 (as state 1st is inactive)
    pR_nj(ryr_initial_state) = 1d0  ! all RyR channel is in state 1 (state 1st is close)

    ALLOCATE(U_Lo(mL), U_Ro(mR))
    U_Lo = iU_L
    U_Ro = iU_R
    
    PRINT *, '+end makeKs_nj()'
    RETURN
  END FUNCTION makeKs_nj

  
END MODULE makeKs_utility
