!!>**** 
!! SYNOPSIS:  
!!   ECC code utility on GPU
!!   
!! CREATOR: 
!!      Hoang-Trong Minh Tuan 
!! AFFILIATION: 
!!      George Mason University 
!! 
!! HISTORY
!!    Rev: 0.1 
!!    0.2: add sum reduction on GPU, min/max on GPU
!!    0.3: add random number generator on GPU using SARU (Hoang-Trong Minh Tuan - Oct, 18, 2010)
!!    0.4: rename from spark_gpu_utility to ecc_gpu_utility
!!    0.5: working version on Fermi (Tuan, Nov 30, 2010)
!!    0.6: add CaB_htrpn (Tuan, Dec 29, 2010)
!!    0.7: fix bug on I-clamp (Tuan,Feb, 01, 2010)
!!    0.8: revise the code, work with PGI 11.3 only (Tuan, Jul, 11, 2011)
!!    0.9: fix bug with PRNG working with PGI 11.7 (Jafri+Tuan, Jul 25,2011)
!!    1.0: rename data to reflect using the new algorithm (Mar 3rd, 2012 - Tuan)
!!_****
MODULE ecc_gpu_utility
  USE precisioninfo
  USE data_global
  USE data_ecc
  USE data_ecc_gpu
  USE rand_utility
  USE hdf5_utility
  USE iso_c_binding
  USE cudafor
  USE cuda_interop
  USE error_codes
  IMPLICIT NONE
#include "macro.h"
!!!!! SAME CHI !!!!!!!!!

#define _USE_CaDS_APPROX_ _NO
!#define _RUNMODE_   _MODE_V_CLAMP
!#define _RUNMODE_ _MODE_RESTING
#define _RUNMODE_ _MODE_I_CLAMP
!#define _RUNMODE_ _MODE_CLOSEDCELL
!#define _DATAMODE_  _IO_2_HDF2
#define _DATAMODE_  _IO_TEST
#define _GPU_MODEL_ _FERMI
#define _USE_PRNG_ _PRNG_SARU
#define _STRETCHMODE_ _SINGLE
!#define _STRETCHMODE_ _PACED
  
#define _BENCHMARK_ _NO
#define _SAFE_MODE_ _NO
  ! DEBUG
#define _DEBUG_   _NO

#if _TEST_NEW_ORPHAN_CRU_ == _YES 
  #define _USE_CaDS_APPROX_ _NO
#endif
  
  INTEGER :: CHUNK_SIZE   !number of records
  INTEGER :: CHUNKS_PER_FILE = 1 ! set to 1 for now
  INTEGER :: FIELDS_PER_RECORD
  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: data_buffer, data_2restartfile
  INTEGER :: num_record = 0
  INTEGER :: file_index = 0
  INTEGER :: RF_FIELDS_PER_RECORD  !in restart file
  INTEGER, PARAMETER :: copy_interval=100 ! how often to copy data back/forth GPU-CPU

  CHARACTER(350) :: fstring, file_out

  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), pinned:: Ca_ds, Ca_jsr
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), pinned :: J_efflux_total, &
       J_refill_total, I_dhpr_total, J_ryr_total, J_ryr_orphan_total,Cads_total
  INTEGER, DIMENSION(:), ALLOCATABLE, pinned :: Nquarks, Nsparks, Nblinks

  REAL(KIND=dp) :: J_dhpr_total, GAIN, accumulate_J_ryr, accumulate_J_dhpr, accumulate_J_ryr_orphan
  REAL(KIND=dp) :: temp_kon_htrpn, temp_kon_ltrpn

  !+ APD
  INTEGER, PARAMETER :: NUM_APD_TYPES = 4 !APD90, APD50, APD20, APD0mV
  REAL(KIND=dp), DIMENSION(NUM_APD_TYPES) :: Vmthreshold  
  REAL(KIND=dp), DIMENSION(NUM_APD_TYPES) :: APD, APstart, APend
  LOGICAL, DIMENSION(NUM_APD_TYPES) :: is_inside_AP  
  INTEGER, DIMENSION(NUM_APD_TYPES) :: APcount
  REAL :: time_start_DiastolicPhase
  LOGICAL :: is_in_DiastolicPhase, is_done_firstbeat
  LOGICAL :: has_not_reset_states, has_not_intro_ROS,has_not_removed_ROS
  INTEGER :: current_step_no
  REAL :: step_size,freq,letgo

  
  !+ Spark/Quark/Blink analysis
  LOGICAL, DIMENSION(:), ALLOCATABLE, device :: flag1, flag2, flag3, flag4
   ! flag1 (spark/quark start), flag2 (blink start), flag3 (spark/quark rise), flag4 (true spark)
  REAL, DIMENSION(:), ALLOCATABLE, device :: ts_spark, ts_blink
  INTEGER, DIMENSION(:), ALLOCATABLE, device :: Nquarks_dev, Nsparks_dev, Nblinks_dev
  !REAL, DIMENSION(:,:), ALLOCATABLE, device :: SparkPeaks, SparkDurations, BlinkDurations
  REAL, DIMENSION(:), ALLOCATABLE, device :: tot_SparkPeaks, tot_SparkDurations, tot_BlinkDurations !total
  !NTEGER, PARAMETER :: MAX_NUM_QUARKS_PER_SEC = 10000, MAX_NUM_SPARKS_PER_SEC = 1000
  INTEGER, DIMENSION(:), ALLOCATABLE, device :: prev_RyRgate
  REAL, DIMENSION(:), ALLOCATABLE, device :: prev_Ca_ds  

  INTEGER :: num_Sparks, num_Blinks
  REAL:: avg_SparkDuration, avg_BlinkDuration, avg_SparkPeak, &
       SFWHM, BFWHM  ! spark/blink full half-width maximum
  INTEGER:: num_Quarks
  REAL :: Spark_Freq
  REAL(KIND=dp) :: tsegment, time_prev1, time_prev2, time_prev_lessfreq, time_start_spark_analysis, time_prev_fb

  !+ device data
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, device :: Ca_ds_dev, Ca_jsr_dev
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, device :: X_r_dev ! vector keeps random-numbers

  !!+ cuda stream
  INTEGER :: stream0, stream1

  !+ device constant data
  INTEGER, PARAMETER :: SIZE_INT_SHARE = 15, SIZE_DP_SHARE = 20
  INTEGER, DIMENSION(SIZE_INT_SHARE), constant :: int_ca ! constant array
  REAL(KIND=dp), DIMENSION(SIZE_DP_SHARE), constant :: dp_ca

  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE, device :: lcc_compKv_act_dev, lcc_compKv_inact_dev, &
       lcc_compKfree_dev, ryr_compKfree_mode1_dev, lcc_compKcal_inact_dev, ryr_compKcal_act_mode1_dev, &
       ryr_compKfree_mode2_dev, ryr_compKcal_act_mode2_dev, ryr_compKfree_dev, ryr_compKRos_dev
  INTEGER, DIMENSION(:,:), ALLOCATABLE, device :: lcc_indxKv_act_dev, lcc_indxKv_inact_dev, &
       lcc_indxKfree_dev, ryr_indxKfree_mode1_dev, lcc_indxKcal_inact_dev, ryr_indxKcal_act_mode1_dev, &
       ryr_indxKfree_mode2_dev, ryr_indxKcal_act_mode2_dev, ryr_indxKfree_dev, ryr_indxKRos_dev

  !! old: size (NSFU, maxnklm+maxnklp-1)
  !! new method: size (NSFU, 1)
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: compPdt_dev
#if _DEBUG_ == _YES
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: compPdt
#endif  
  INTEGER, ALLOCATABLE, DIMENSION(:) :: Nclose_RyRs_mode1, Nclose_RyRs_mode2, Nopen_RyRs_mode1,Nopen_RyRs_mode2
  ! size (NSFU, maxnklm+maxnklp-1)
  INTEGER, ALLOCATABLE, DIMENSION(:, :), device :: indexPdt_dev
  ! size (irow_R)
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device ::  U_Ro_dev, U_Lo_dev
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:,:), device:: U_RyR_State_dev
  ! size (NSFU)  - keep current state index at each SFU
  INTEGER, ALLOCATABLE, DIMENSION(:), device ::  CRU_states_dev
  INTEGER, ALLOCATABLE, DIMENSION(:), device :: Nclose_RyRs_mode1_dev, Nclose_RyRs_mode2_dev, Nopen_RyRs_mode1_dev, Nopen_RyRs_mode2_dev
  ! size (NSFU)
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: RYRgate_dev, DHPRgate_dev
  ! size (NSFU)
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: J_efflux_dev
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: J_refill_dev
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: I_dhpr_dev
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), pinned :: I_dhpr
  !REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: J_dhpr_dev, J_ryr_dev
  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), device :: J_ryr_dev

  !+buffering
  REAL(KIND=dp) :: beta_myo
  !+voltage
  REAL(KIND=dp) :: Vm &    ! membrane voltage
       , E_Ca, E_Na, E_K  ! reversal potential
  !+ ionic currents
  REAL(KIND=dp) :: I_ncx, I_Na, I_K1, I_bK, I_bCa, I_NaK, I_pmca, &
       I_app, I_bNa, I_Ktof, I_Kss, I_Ktos, I_dhpr_nj_T
  !+ global calcium ionic fluxes
  REAL(KIND=dp) :: J_ncx, J_serca, J_ryr_nj_T, J_Ca_sarmem
  REAL(KIND=dp) ::  J_leak, J_htrpn, J_ltrpn, J_ryr_orphan
  INTEGER :: NUM_STATES_NJ
  INTEGER :: NUM_FIELDS4CURRENTS, NUM_INDIVIDUAL_FIELDS

  REAL(KIND=dp) :: CaB_htrpn, CaB_ltrpn,CaB_trpn, CaF
  
  !+ Ros
  REAL(KIND=dp) :: RosDS, RosI, baseRos_DS,baseRos_I,baseRos_prod,gsh 
  REAL(KIND=dp) :: Vefflux_Ros_T, v_RosRemoval, v_RosRemoval_i, v_RosRemoval_ds, v_RosRemoval_base
  REAL(KIND=dp) :: Vefflux_Ros, lambda_ds_Ros, inv_lambda_ds_Ros, effective_NSFU
  REAL(KIND=dp) :: RosEffluxed, RosRemoved, RosProduced  
  REAL(KIND=dp) :: t_stretch, t_unstretch, stretch_area
  REAL(KIND=dp) :: t_this_cycle,RosProd_atstart_unstr,RosProd_atstart_str
  REAL(KIND=dp) :: cycle_size, prod_rise_time,Ros_relaxTime
  REAL(KIND=dp) :: t_stretch_start, stretch_duration, t_stretch_end, t_stretch_after_I_stim
  REAL(KIND=dp) :: K_Ros, eta_Ros,bigRos,V_ds_Ros

  !+ summary data
  REAL(KIND=dp), DIMENSION(100) :: summary_data ! to store GAIN, APD...
  INTEGER :: count_summarydata
  CHARACTER(40), DIMENSION(100) :: description_summarydata
  
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: p_nj
  CHARACTER(8)  :: date
  INTEGER,DIMENSION(8) :: values
  
  TYPE(HDFFILE) :: file_struct

#if _TEST_ORPHAN_CRU_ == _YES
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pR_orphan, tmp_pR_orphan
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: chiC_orphan, chiO_orphan
#endif
  LOGICAL :: is_being_clamped

#if _USE_PRNG_ == _PRNG_SARU
  TYPE(c_ptr) :: saruObj
#elif _USE_PRNG_ == _PRNG_MT
  TYPE(c_ptr) :: d_status
  TYPE(c_ptr) :: d_texture
#elif _USE_PRNG_ == _PRNG_MTCPU
  INTEGER, DIMENSION(:), ALLOCATABLE :: mt_STATE(:)
  REAL(KIND=DP) :: mt_LB, mt_UB   !lower-bound, upper-bound mersenne twister
  INTEGER :: INFO
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: X_r
#elif _USE_PRNG_ == _PRNG_CURAND
  
#endif

  TYPE(c_ptr) :: odata_d
  REAL(KIND=dp) :: dp_arg1, dp_arg2
   
  SAVE
CONTAINS
  !!>****f* 
  !! SYNOPSIS  
  !!   Assign frequently used read-only data to device constant memory
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  SUBROUTINE set_device_constant()
    IMPLICIT NONE

    int_ca = (/ NSFU, N_R, ryr_Ncol_Kcal_act_mode1, ryr_Ncol_Kfree_mode1, lcc_Ncol_Kv_act, & !5
         lcc_Ncol_Kv_inact, lcc_Ncol_Kcal_inact, lcc_Ncol_Kfree, lcc_Ncol_Kv_act+lcc_Ncol_Kv_inact+ &
         lcc_Ncol_Kfree+ryr_Ncol_Kfree_mode1+lcc_Ncol_Kcal_inact+ryr_Ncol_Kcal_act_mode1-5, irow_R, & ! 10
         ryr_Ncol_Kcal_act_mode2, ryr_Ncol_Kfree_mode2, ryr_Ncol_Kfree,ryr_Ncol_KRos, ryr_Ncol_Kcal_act_mode1+ &
         ryr_Ncol_Kfree_mode1+lcc_Ncol_Kv_act+lcc_Ncol_Kv_inact+lcc_Ncol_Kcal_inact+lcc_Ncol_Kfree+ &
         ryr_Ncol_Kcal_act_mode2+ryr_Ncol_Kfree_mode2+ ryr_Ncol_Kfree+ryr_Ncol_Kros-9 & ! 15
         /)
    !beta_ds : 12 --> 6
    dp_ca = (/ -Ej*0.5d0, Ecc, Eoo, k_jsr0, k_jsr1, & ! 5
         beta_ds, inv_lambda_ds, inv_lambda_jsr, v_ryr, v_refill, & ! 10
         v_efflux, 0.d0, P_dhpr, eta_RyR , eta_LCC, &  ! 15
         F2zCa2_RT, inv_2FVmyo, 0.341d0*Ca_o, Km_jsr, B_jsrT &  ! 20
         /)
  END SUBROUTINE set_device_constant

  !!>****f* 
  !! SYNOPSIS  
  !!   Initialize data for the simulation
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  FUNCTION init_Simulation() RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    !local     
    LOGICAL :: plog
    INTEGER :: iserror

    res = print_mode_info()

    !! memory allocation
    !+ device
    ALLOCATE(lcc_compKv_act_dev(irow_L,0:lcc_Ncol_Kv_act), &
         lcc_compKv_inact_dev(irow_L,0:lcc_Ncol_Kv_inact), &
         lcc_compKcal_inact_dev(irow_L,0:lcc_Ncol_Kcal_inact) , &
         lcc_compKfree_dev(irow_L,0:lcc_Ncol_Kfree), &
         ryr_compKfree_mode1_dev(irow_R,0:ryr_Ncol_Kfree_mode1), &
         ryr_compKcal_act_mode1_dev(irow_R,0:ryr_Ncol_Kcal_act_mode1),&
         ryr_compKfree_mode2_dev(irow_R,0:ryr_Ncol_Kfree_mode2), &
         ryr_compKcal_act_mode2_dev(irow_R,0:ryr_Ncol_Kcal_act_mode2),&
         ryr_compKfree_dev(irow_R,0:ryr_Ncol_Kfree), &
         ryr_compKRos_dev(irow_R,0:ryr_Ncol_KRos), STAT=iserror)

    IF (iserror) THEN
       PRINT *, "irow_R: ",irow_R, irow_L
       PRINT *, "ryr_Ncol_Kfree_mode1: ",ryr_Ncol_Kfree_mode1
       PRINT *, "ryr_Ncol_Kcal_act_mode1: ",ryr_Ncol_Kcal_act_mode1
       PRINT *, "ryr_Ncol_Kfree_mode2: ",ryr_Ncol_Kfree_mode2
       PRINT *, "ryr_Ncol_Kcal_act_mode2: ",ryr_Ncol_Kcal_act_mode2
       PRINT *, "ryr_Ncol_Kfree: ",ryr_Ncol_Kfree
       PRINT *, "ryr_Ncol_KRos: ",ryr_Ncol_KRos
      


       WRITE(error_message, *) "FAILED: device memory allocation compK"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF
    ALLOCATE(lcc_indxKv_act_dev(irow_L, 0:lcc_Ncol_Kv_act), &
         lcc_indxKv_inact_dev(irow_L, 0:lcc_Ncol_Kv_inact), lcc_indxKcal_inact_dev(irow_L, 0:lcc_Ncol_Kcal_inact), &
         lcc_indxKfree_dev(irow_L, 0:lcc_Ncol_Kfree), &
         ryr_indxKfree_mode1_dev(irow_R, 0:ryr_Ncol_Kfree_mode1),  &
         ryr_indxKcal_act_mode1_dev(irow_R, 0:ryr_Ncol_Kcal_act_mode1), &
         ryr_indxKfree_mode2_dev(irow_R, 0:ryr_Ncol_Kfree_mode2),  &
         ryr_indxKcal_act_mode2_dev(irow_R, 0:ryr_Ncol_Kcal_act_mode2), &
         ryr_indxKfree_dev(irow_R, 0:ryr_Ncol_Kfree),  &
         ryr_indxKRos_dev(irow_R, 0:ryr_Ncol_KRos), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation indxK"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    ALLOCATE(U_Ro_dev(irow_R), U_Lo_dev(irow_L),U_RyR_State_dev(irow_R,mR), STAT = iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation U_Ro, U_Lo, U_RyR_State "
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    ALLOCATE(Ca_ds_dev(NSFU), Ca_jsr_dev(NSFU), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation Cajsr, Cads "
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF
    
    ALLOCATE(J_efflux_dev(NSFU), J_refill_dev(NSFU), J_ryr_dev(NSFU), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation fluxes Jrefill, Jefflux, Jryr"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    !ALLOCATE(I_dhpr_dev(NSFU), J_dhpr_dev(NSFU), J_ryr_dev(NSFU))
    ALLOCATE(I_dhpr_dev(NSFU), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation Idhpr"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    ALLOCATE(RYRgate_dev(NSFU),  DHPRgate_dev(NSFU), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation RyRgate, DHPRgate"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    ALLOCATE(Nclose_RyRs_mode1_dev(NSFU), Nclose_RyRs_mode2_dev(NSFU), Nopen_RyRs_mode1_dev(NSFU), Nopen_RyRs_mode2_dev(NSFU), STAT=iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation Nclose_RyRs, Nopen_RyRs"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    ALLOCATE(CRU_states_dev(NSFU), STAT = iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation state indices"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    
    ALLOCATE(compPdt_dev(NSFU), STAT = iserror)
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation compPdt"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF
    
    ALLOCATE(X_r_dev(NSFU*copy_interval), STAT = iserror)  ! for random numbers
    IF (iserror) THEN
       WRITE(error_message, *) "FAILED: device memory allocation X_r random numbers"
       res = ERR_CUDA_ALLOCATE_MEM
       RETURN
    END IF

    !+ pinned memory
    ALLOCATE(Ca_ds(NSFU), Ca_jsr(NSFU), PINNED=plog)
    IF (.NOT. plog) THEN
       WRITE(log_string,*) "FAILED: Ca_ds, Ca_jsr on PINNED memory"
       CALL write2log(log_string)
    END IF

    ALLOCATE(J_efflux_total(1), J_refill_total(1), I_dhpr_total(1), J_ryr_total(1), J_ryr_orphan_total(1), PINNED=plog)
    IF (.NOT. plog) THEN
       WRITE(log_string,*) "FAILED: J_efflux_total, J_refill_total, I_dhpr_total on PINNED memory"
       CALL write2log(log_string)
    END IF

    ALLOCATE(Cads_total(1), PINNED=plog)
    IF (.NOT. plog) THEN
       WRITE(log_string,*) "FAILED: Cads_total on PINNED memory"
       CALL write2log(log_string)
    END IF
    ALLOCATE(I_dhpr(NSFU), PINNED=plog)
    IF (.NOT. plog) THEN
       WRITE(log_string,*) "FAILED: I_dhpr on PINNED memory"
       CALL write2log(log_string)
    END IF
    
    !+ host memory
    ALLOCATE(CRU_states(NSFU), LCCcluster_states(NSFU), RYRcluster_states(NSFU))
    ALLOCATE(RYRgate(NSFU),  DHPRgate(NSFU))
    ALLOCATE(Nclose_RyRs_mode1(NSFU), Nclose_RyRs_mode2(NSFU), Nopen_RyRs_mode1(NSFU), Nopen_RyRs_mode2(NSFU)) 
    ALLOCATE(p_nj(irow_L_nj*irow_R_nj))

#if _DEBUG_ == _YES
    ALLOCATE(compPdt(NSFU))
#endif    

    !! spark analysis data
    ALLOCATE(flag1(NSFU), flag2(NSFU), flag3(NSFU), flag4(NSFU))
    flag1 = .FALSE.
    flag2 = .FALSE.
    flag3 = .FALSE.
    flag4 = .FALSE.
    ALLOCATE(ts_spark(NSFU), ts_blink(NSFU))
    ts_spark = +0d0
    ts_blink = +0d0
    ALLOCATE(Nquarks_dev(NSFU), Nsparks_dev(NSFU), Nblinks_dev(NSFU))
    Nquarks_dev = 0
    Nsparks_dev = 0
    Nblinks_dev = 0

    ALLOCATE(Nquarks(NSFU), Nsparks(NSFU), Nblinks(NSFU))

    ALLOCATE(tot_SparkPeaks(NSFU), tot_SparkDurations(NSFU), tot_BlinkDurations(NSFU))
    tot_SparkDurations = 0d0
    tot_SparkPeaks = 0d0
    tot_BlinkDurations = 0d0
    ALLOCATE(prev_RyRgate(NSFU), prev_Ca_ds(NSFU))

    !copy from host to device
    lcc_compKv_act_dev = lcc_compKv_act
    lcc_compKv_inact_dev = lcc_compKv_inact
    lcc_compKfree_dev = lcc_compKfree
    lcc_compKcal_inact_dev = lcc_compKcal_inact
    ryr_compKcal_act_mode1_dev = ryr_compKcal_act_mode1
    ryr_compKfree_mode1_dev = ryr_compKfree_mode1
    ryr_compKcal_act_mode2_dev = ryr_compKcal_act_mode2
    ryr_compKfree_mode2_dev = ryr_compKfree_mode2
    ryr_compKfree_dev = ryr_compKfree
    ryr_compKRos_dev = ryr_compKRos
    lcc_indxKv_act_dev = lcc_indxKv_act
    lcc_indxKv_inact_dev = lcc_indxKv_inact
    lcc_indxKfree_dev = lcc_indxKfree
    lcc_indxKcal_inact_dev = lcc_indxKcal_inact
    ryr_indxKcal_act_mode1_dev = ryr_indxKcal_act_mode1
    ryr_indxKfree_mode1_dev = ryr_indxKfree_mode1
    ryr_indxKcal_act_mode2_dev = ryr_indxKcal_act_mode2
    ryr_indxKfree_mode2_dev = ryr_indxKfree_mode2
    ryr_indxKfree_dev = ryr_indxKfree
    ryr_indxKRos_dev = ryr_indxKRos
    U_Ro_dev = U_Ro  !irow_R length
    U_Lo_dev = U_Lo  !irow_L length
    U_RyR_State_dev = U_RyR_State

    I_app = +0d0  ! uA/cm^2
    Vm = V_rest   ! mV

    is_being_clamped = .FALSE.
    is_inside_AP = .FALSE.
    is_EndSim = .FALSE.

    !!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! data setting
    CALL loaddata_from_conf()
    PRINT *,"loaddata_from_conf() SUCCESSFUL"
    IF (TRIM(restart_file) .NE. "") THEN
       ! overload the setting from parameter file
       CALL loaddata_from_restartfile()
    PRINT *,"loaddata_from_restartfile() SUCCESSFUL"
       CALL state_breakdown_all(CRU_states, irow_R, LCCcluster_states, RYRcluster_states)
    END IF
    
    ! initialize gpu data
    Ca_ds_dev = Ca_ds
    Ca_jsr_dev = Ca_jsr
      
    compPdt_dev = +0d0

    IF (TRIM(restart_file) .EQ. "") THEN
       ! initialize release site's state
       CALL init_CRU_states(lcc_initialstate, ryr_initialstate, N_L, mL, N_R, mR, &
            CRU_states, LCCcluster_states, RYRcluster_states)
       PRINT *, "init_CRU_states() SUCCESSFUL"
    ENDIF
        
    PRINT *, "after the IF block for init_CRU_states()"
    !copy from host to device
    CRU_states_dev = CRU_states
    RYRgate_dev = U_Ro(RYRcluster_states) !IMPORTANT: now we need to deal with homogeneous cluster state index
    DHPRgate_dev = U_Lo(LCCcluster_states)
    Nclose_RyRs_mode1_dev = U_RyR_State(RYRcluster_states(:),1)
    Nclose_RyRs_mode2_dev = U_RyR_State(RYRcluster_states(:),4)
    Nopen_RyRs_mode1_dev = U_RyR_State(RYRcluster_states(:),2)
    Nopen_RyRs_mode2_dev = U_RyR_State(RYRcluster_states(:),3)

    ! Pseudo-random-number-generator
    CALL init_PRNG()
    !PRINT *, "init_PRNG() SUCCESSFUL"
    ! Define simulation parameters
    t_start = t_start + real_timestart ! NOTE: read-in value for t_start is 0
    time = t_start
    t_end = t_end + real_timestart
    !TUAN: spark statistics are derived from (maximum) the last 5 seconds of each simulation
    !+ when to start spark analysis
    time_start_spark_analysis = MAX(t_start, t_end - 5.0d0)
    IF (duration_write_lastbeats .GT. t_end-t_start) THEN
       WRITE(error_message, *) "Duration to write lastbeats longer than simulation setting"
       res = ERR_FILE_INVALID_DATA
       RETURN
    ENDIF
    time_prev1 =  t_end - duration_write_lastbeats - write_interval ! a time tracker for data I/O
    time_prev_lessfreq = time_prev1
    time_prev_fb = t_start
    time_prev2 = time_start_spark_analysis ! a time tracker to find integrate flux
    
    tsegment = 1.d0/(Hz)
    IF (Hz .GT. 1) THEN
       t_V_stim_start = tsegment/2
       t_V_stim_end = t_V_stim_start + t_V_stim_duration
       t_I_stim_start = tsegment/2
       t_I_stim_end = t_I_stim_start + t_I_stim_duration
    END IF
    
    ! Vm-dependent rate transitions
    INCLUDE "Vm_dependent.inc"  !!kva,kvi
    INCLUDE "dp_args.inc"  !! dp_arg1, dp_arg2

    odata_d = initreduction_double()  !prepare for the reduction

    res = setup_KernelConfiguration()
    IF (res .NE. SUCCESS) RETURN
    
    ! copy data to constant memory
    CALL set_device_constant()
    PRINT *, "set_device_constant() SUCCESSFUL"
    !! Initial setting values
    beta_myo = 0d0;  
    J_efflux_total(1) = 0d0; J_efflux_dev = 0.0d0
    J_leak = 0d0;
    J_refill_total(1) = 0d0;     J_refill_dev = 0.0d0
    J_Ca_sarmem = 0d0;   J_ryr_nj_T = 0d0;  J_serca = 0d0; J_htrpn = 0d0; J_ltrpn = 0d0
    I_Na = 0d0; I_dhpr_nj_T = 0d0
    I_dhpr_total(1) = 0d0;     I_dhpr_dev = 0d0; I_dhpr = 0d0
    J_ryr_orphan = 0d0
    J_ryr_total(1) = 0d0;     J_ryr_dev = +0d0; 
    J_ryr_orphan_total(1) = 0d0
    Cads_total(1) = 0d0
    GAIN = 0d0
    accumulate_J_dhpr = 0.0d0
    accumulate_J_ryr = 0.0d0; accumulate_J_ryr_orphan = 0.0d0
    APD = 0d0; APstart = 0d0; APend = 0d0;
    Vmthreshold = (/ V_rest * 0.9d0, V_rest * 0.5d0, V_rest * 0.2d0, +0.0d0 /)
    APcount = +0
    is_in_DiastolicPhase = .FALSE.
    is_done_firstbeat = .FALSE.
    has_not_reset_states = .TRUE.
    I_Na = 0.0d0
    I_Ktof = 0.0d0
    I_Ktos = 0.0d0
    I_K1 = 0.0d0
    I_Kss = 0.0d0
    I_bNa = 0.0d0
    I_bK = 0.0d0
    I_bCa = 0.0d0
    I_NaK = 0.0d0
    I_pmca = 0.0d0
    I_dhpr_nj_T = 0.0d0

!!!initializing ROS parameters

    !has_not_intro_ROS = .TRUE.
    !has_not_removed_ROS = .TRUE.
    bigRos = 0.0d0
    !base_Ros = 0.0d0
    RosEffluxed = 0
    RosRemoved = 0
    RosProduced = 0
    t_this_cycle = 0.0d0
    
    eta_Ros = 4.0d-1
    !K_Ros = 2.0d-1
    Vefflux_Ros_T = 1.25d0/4.0d0!!2.6d-5
    stretch_area = 8.0d-1
    prod_rise_time = 1.2d-1  
    Ros_relaxTime = 6.0d-1 !1.2d-1
    v_RosRemoval = 2.0d0
    gsh = 1.0d0
    baseRos_DS = 1.0d-2
    baseRos_I = 1.0d-2

    effective_NSFU = NSFU*stretch_area;
    Vefflux_Ros = Vefflux_Ros_T/effective_NSFU;
    V_ds_Ros = V_ds_T/(effective_NSFU)
    lambda_ds_Ros = V_ds_Ros/V_myo_T
    inv_lambda_ds_Ros = 1.0d0/lambda_ds_Ros
    v_RosRemoval_i = v_RosRemoval*gsh
    v_RosRemoval_ds = lambda_ds_Ros*v_RosRemoval*gsh
    v_RosRemoval_base = v_RosRemoval*baseRos_I
    !baseRos_prod = Vefflux_Ros*(baseRos_DS-baseRos_I) + v_RosRemoval_ds*(baseRos_DS)
    baseRos_prod = 0.0d0
    
    RosDS = baseRos_DS
    RosI = baseRos_I
    RosProd_atstart_unstr = 0.0d0
    RosProd_atstart_str = baseRos_prod
    !inv_K_removalRos = 1.0d0/(K_removalRos)
    PRINT *,"effective_NSFU: ", effective_NSFU
    PRINT *,"inv_lambda_ds: ",inv_lambda_ds_Ros
    PRINT *,"Vefflux_Ros: ",Vefflux_Ros

   !!initialize for ROS/stretching 
    letgo = 1.0d1
    t_stretch_start = 0.0d0
    stretch_duration = 1.0d1 !2.0d-1 !1.0d1

    !! _MODE_RESTING
      !!! initialize freq for _MODE_RESTING && _STRETCHMODE_ _PACED
    freq = 1.0d0
    temp_kon_htrpn = 2.0d1*(5.5d-1) !(7.0d-1) !(6.67d-1) !(5.0d-1) !1.0d2
      !koff_htrpn = 6.60d-2
    temp_kon_ltrpn = 4.0d1*(5.5d-1)  !(7.0d-1) !(6.67d-1) !(5.0d-1) !4.00d0 !3.9d1
     
     !!_MODE_I_CLAMP
      !!!initialize starting time of stretch for _MODE_I_CLAMP 
    !t_stretch_after_I_stim = 0.0d0!9.5d-1 !6.0d-1 !!How long after I-stim should the Ros-stim be started (stretching be done)?

#if _STRETCHMODE_ == _SINGLE
    cycle_size = 2.0d0*(stretch_duration)
    !cycle_size = 1.5d2
    !t_stretch_end = stretch_duration

  !#if _RUNMODE_ == _MODE_I_CLAMP
    !t_stretch_start = t_I_stim_start + t_stretch_after_I_stim
    !t_stretch_start = 0.0d0
  !#else
   ! t_stretch_start = 0.0d0
  !#endif
#else
  #if _RUNMODE_ == _MODE_I_CLAMP
    cycle_size = 1.0d0/(Hz)    
    !t_stretch_start = t_I_stim_start + t_stretch_after_I_stim
    !t_stretch_end = stretch_duration
  #else
    cycle_size = 1.0d0/freq
    stretch_duration = cycle_size/2.0d0
    !t_stretch_end = cycle_size/2.0d0
    !t_stretch_start = 0.0d0
  #endif
#endif

    PRINT *,"cycle_size: ",cycle_size
    PRINT *,"Stretch at: ",t_stretch_start," sec"
!!#if _RUNMODE_ == _MODE_I_CLAMP
!    letgo = 1.0d1
!    t_Ros_stim_start = t_I_stim_start + t_Ros_after_I_stim
!  #if _STRETCHMODE_ == _SINGLE
!    cycle_size = 2.0d0*(t_unstretch - t_stretch)
!    t_Ros_stim_end = cycle_size/2.0d0
!  #else
!    cycle_size = 1.0d0/(Hz)
!    t_Ros_stim_end = t_Ros_stim_start + t_Ros_stim_duration
!  #endif
!    PRINT *,"cycle_size: ",cycle_size
!    PRINT *,"freq: ",Hz, "Stretch at: ",t_Ros_stim_start," sec"
!#else

!  #if _STRETCHMODE_ == _SINGLE
!    cycle_size = 2.0d0*(t_unstretch - t_stretch)
!    letgo = t_stretch
!    PRINT *,"cycle size: ",cycle_size 
!  #else
!    cycle_size = 1.0d0/freq
!    PRINT *,"cycle size: ",cycle_size
!  #endif
!    t_Ros_stim_start = 0.0d0
!    t_Ros_stim_end = cycle_size/2.0d0
!#endif

    ! I/O init
    !!!TUAN: check to make sure I/O for firstbeat working properly
    IF (write_firstbeat .NE. 0) THEN
       CALL data_io_init_fb()
       CALL data_io_write_fb()
    ELSE
       CALL data_io_init()
    ENDIF
    !PRINT *,"irow_R: ",irow_R
    !PRINT *,"mR: ",mR

    PRINT *, "+end init_Simulation()"
    !STOP
  END FUNCTION init_Simulation

  !!>****f* 
  !! SYNOPSIS  
  !!     Calculate the calcium concentration using GPU 
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  attributes(host) SUBROUTINE continue_Simulation()
    IMPLICIT NONE
    !+ other local
    INTEGER, PARAMETER :: MAX_ITER=600000000
    INTEGER :: ii, nsteps, err
            
    PRINT *, "start continue_Simulation()"
    
    nsteps = copy_interval
    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !! TIMING LOOP
    OUTERLOOP: DO ii = 1, MAX_ITER
#if _BENCHMARK_ == _YES
       PRINT *, "ii", ii
       IF (ii .EQ. 3) THEN
          EXIT OUTERLOOP
       ENDIF
#endif

       !! run a number of time steps
       err = simulate_NextSteps(nsteps, is_EndSim)
#if _SAFE_MODE_ == _YES
       IF (err .NE. SUCCESS) THEN
          WRITE(log_string, *) "Error in simulate_NextSteps()"
          CALL write2log(log_string)
          EXIT OUTERLOOP
       END IF
#endif

       IF (is_EndSim) EXIT OUTERLOOP
    END DO OUTERLOOP

    PRINT *, "++end  continue_Simulation()"

    !! writing the last time
    CALL data_io_finalize()

        
#if _BENCHMARK_ == _NO
    !!TUANHOANG
    CALL spark_summarize()
#endif
    
    RETURN
  END SUBROUTINE continue_Simulation

  SUBROUTINE finalize_Simulation
    !!clean data on GPU
    CALL finalizereduction_double_tesla(odata_d)
  #if _USE_PRNG_ == _PRNG_SARU    
    CALL finalize_saru(saruObj)
  #endif
    CALL free_device()
  END SUBROUTINE finalize_Simulation

  !!>****f* ecc_gpu_utility/simulate_NextSteps
  !! SYNOPSIS  
  !!     Perform a number of time steps (i.e. "nsteps" steps)
  !!     NOTE: make sure nsteps .LE.  copy_interval
  !!   
  !! INPUTS
  !!   nsteps
  !! OUTPUTS 
  !!   is_EndSim : .TRUE. if the simulation is complete
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  INTEGER FUNCTION simulate_NextSteps(nsteps, is_EndSim) RESULT(res)
    IMPLICIT NONE
    LOGICAL, INTENT(INOUT) :: is_EndSim
    INTEGER, INTENT(IN):: nsteps
    !!local 
    INTEGER :: iinner, cuErr, err
    REAL(kind=dp), PARAMETER :: ramp_time = 1.0d-3  ! should be less than 1ms
    REAL(KIND=dp), PARAMETER :: tau_time = 1.0d-3
    REAL(KIND=dp) :: tresidue, dt, invdt, min_Pdt, dtmax
    INTEGER :: ii
    REAL(KIND=dp) :: dummy1, dummy2,tempRos,tempRos1
    !+ parameters for gating variables
    REAL(KIND=dp) :: TCa_i2, TCa_sr2, vcycle_pump
    REAL(KIND=dp), DIMENSION(irow_L_nj) :: tmp_pL_nj
    REAL(KIND=dp), DIMENSION(irow_R_nj) :: tmp_pR_nj




    res = SUCCESS

    IF (is_EndSim) RETURN

#if _SAFE_MODE_ == _YES
    IF (nsteps .GT. copy_interval) THEN
       PRINT *, "ERROR: simulate_NextSteps\n nsteps is larger than copy_interval"
       err = ERR_INVALID_DATA_SETTING
       RETURN
    END IF
#endif
       
#if _USE_PRNG_ == _PRNG_SARU
    !!SARU pseudo-RNG
    CALL gen_rand_saru_c2050(saruObj, X_r_dev, NSFU, nsteps)
#elif _USE_PRNG_ == _PRNG_MT
    CALL generate_mt_data(d_status, d_texture, NSFU*nsteps, X_r_dev)
#elif _USE_PRNG_ ==  _PRNG_MTCPU
    CALL DRANDUNIFORM(NSFU*nsteps,mt_LB,mt_UB,mt_STATE,X_r,INFO)
    X_r_dev = X_r
#endif

    INNERLOOP: DO iinner = 1, nsteps  ! inner loop
       ! CHECK FOR COMPLETE SIMULATION
       IF (time .GT. t_end) THEN
          is_EndSim = .TRUE.
          EXIT INNERLOOP
       ENDIF
#if _BENCHMARK_ == _YES
       !PRINT *, "iinner", iinner
#endif
       
#if (_RUNMODE_ == _MODE_I_CLAMP || _RUNMODE_ == _MODE_V_CLAMP)
       tresidue = MODULO(time, tsegment)
 #if _RUNMODE_ == _MODE_V_CLAMP
       !V-clamp
       !INCLUDE "Vclamp_method1.inc"
       !INCLUDE "Vclamp_method2.inc"
       INCLUDE "Vclamp_method3.inc"
       !INCLUDE "Vclamp_method4.inc"
       !INCLUDE "Vclamp_method5.inc"
       
 #elif _RUNMODE_ == _MODE_I_CLAMP
       !I-clamp
       IF (tresidue.GT. t_I_stim_start .AND. tresidue.LT. t_I_stim_end) THEN
          IF (.NOT. is_being_clamped) THEN
             !current-clamp
             I_app = I_stim
             is_being_clamped = .TRUE.
          END IF
       ELSEIF (is_being_clamped .EQ. .TRUE.) THEN
          I_app = +0d0             
          is_being_clamped = .FALSE.
       END IF
       INCLUDE "Vm_dependent.inc"  !!kva,kvi
       INCLUDE "dp_args.inc"  !! dp_arg1, dp_arg2
 #endif

#endif

       
       !IF(time<10) THEN
       !   Ros = 0.0d0       
       !ELSE IF(time >10 .AND. has_not_intro_ROS) THEN
       !   Ros = 2.0d1
       !   has_not_intro_ROS = .FALSE.
       !ELSE IF(time >20 .AND. has_not_removed_ROS) THEN
       !   Ros = 0.0d0
       !   has_not_removed_ROS = .FALSE.
       !END IF
 

       !!find rate transition
#if _TEST_NEW_ORPHAN_CRU_ == _YES 
       CALL get_compPdt_2b <<<gridSize, blocksize>>> (kva, kvi, NSFU_normal, Ej_orphan)
#else
       CALL get_compPdt_2 <<<gridSize, blocksize>>> (kva, kvi,bigRos)
#endif

       cuErr = cudaDeviceSynchronize()

#if _SAFE_MODE_ == _YES
       IF (cuErr .NE. cudaSuccess) THEN
          WRITE(log_string,*) "Kernel Error: get_compPdt_2"
          CALL write2log(log_string)
          WRITE(error_message,*) cudaGetErrorString(cudaGetLastError())
          res = CUDA_ERROR
          RETURN
       END IF
#endif
       
       !TUAN: use min on GPU
       min_Pdt = min_reduce_double_2host(compPdt_dev(1:NSFU), NSFU, odata_d)
       
#if _DEBUG_ == _YES
       compPdt = compPdt_dev
       IF (min_Pdt .NE. MINVAL(compPdt)) THEN
          WRITE(log_string,*) "Error: min_reduce_double_2host not working"
          CALL write2log(log_string)
          STOP
       END IF
#endif
       
       dtmax =-0.1d0/min_Pdt
       dt = MAX(MIN(dtmax,dt_UB),dt_LB)
       
!!$#if _RUNMODE_ == _MODE_V_CLAMP
!!$       IF (tresidue .GT. t_V_stim_end .AND. tresidue .LE. t_V_stim_end + ramp_time) THEN
!!$          dt = 1d-8
!!$       ENDIF
!!$#endif       
       time = time + dt
       invdt = 1.0d0/dt
       
    !IF(time .LT. (letgo+t_stretch_start)) THEN
    !  RosDS = base_Ros
    !ELSE
    !   t_this_cycle = MOD((time-letgo-t_stretch_start),cycle_size)
       CALL calc_new_RosDS(dt, RosProduced, RosEffluxed)
       CALL calc_new_RosI(dt, RosEffluxed)
    !ENDIF
       tempRos = LOG(RosDS)
       !bigRos = EXP(tempRos*eta_Ros)
       tempRos1 = EXP(tempRos*eta_Ros)
       !bigRos = tempRos1/(K_Ros+tempRos1)
       bigRos = tempRos1

       !!TUANHOANG
#if _TEST_NEW_ORPHAN_CRU_ == _YES 
       !CALL updateSystem_2b <<<gridSize, blocksize>>> &
            !(X_r_dev((iinner-1)*NSFU+1:iinner*NSFU), &
            !dt, invdt, Ca_nsr, Ca_myo, Vm, kva, kvi, dp_arg1, dp_arg2, &
            !NSFU_normal, inv_lambda_ds_orphan, Km_myo, B_myoT, Ej_orphan)
#else
       CALL updateSystem_2 <<<gridSize, blocksize>>> &
            (X_r_dev((iinner-1)*NSFU+1:iinner*NSFU), &
            dt, invdt, Ca_nsr, Ca_myo, Vm, kva, kvi, dp_arg1, dp_arg2, bigRos)
#endif
          
#if _TEST_ORPHAN_CRU_ == _YES
       J_ryr_orphan = SUM(pR_orphan*U_Ro_orphan, MASK=(U_Ro_orphan .GE. 1)) * v_ryr * (Ca_nsr - Ca_myo)
       chiC_orphan = EXP(-Ej*0.5d0*(dfloat(N_R-U_Ro_orphan)*Ecc - dfloat(U_Ro_orphan-1)*Eoo))
       chiO_orphan = EXP(-Ej*0.5d0*(dfloat(U_Ro_orphan)*Eoo - dfloat(N_R-U_Ro_orphan-1)*Ecc))
       tmp_pR_orphan = pR_orphan
       DO ii = 1, irow_R_orphan
          pR_orphan(ii) = tmp_pR_orphan(ii) + dt * SUM(tmp_pR_orphan(:)* (chiC_orphan(ii) * aKm2_orphan(:,ii) + chiO_orphan(ii) * (k_jsr0*Ca_nsr+k_jsr1)*Ca_myo**eta_RyR * aKp2_orphan(:,ii)))
       END DO
#endif

#if _RUNMODE_ != _MODE_CLOSEDCELL          
       !!Resting potential
       E_Na = RT_F*LOG(Na_o/Na_i)
       E_K = RT_F*LOG(K_o/K_i)
       E_Ca = RT_FzCa*LOG(Ca_o/Ca_myo)
       
       !! Fast Na current (Luo-Rudy)
       IF (Vm .GE. -40.d0) THEN
          a_h = 1d3*0
          a_j = 1d3*0
          b_h = 1d3*1/(0.13d0*(1+EXP((Vm+10.66d0)/-11.1d0)))
          b_j = 1d3*0.3d0*(EXP(-2.535d-7*Vm)/(1d0+EXP(-0.1d0*(Vm+32d0))))
       ELSEIF (Vm .LT. -40.d0) THEN 
          a_h = 1d3*(0.135d0*EXP((80d0+Vm)/-6.8d0))
          a_j = 1d3*(-1.2714d5*EXP(0.2444d0*Vm)-3.474d-5*EXP(-0.04391d0*Vm))* &
               ( (Vm+37.78d0)/(1d0+EXP(0.311d0*(Vm+79.23d0))))
          b_h = 1d3*(3.56d0*EXP(0.079d0*Vm)+3.1d5*EXP(0.35d0*Vm))
          b_j = 1d3*(0.1212d0*EXP(-0.01052d0*Vm)/(1d0+EXP(-0.1378d0*(Vm+40.14d0))))
       ELSE ! NaN blankxo
          WRITE(log_string,*) Vm
          CALL write2log(log_string)
          WRITE(log_string,*) "Error: NaN value Vm"
          CALL write2log(log_string)
          WRITE(log_string, *) time, V_rest, Vm
          CALL write2log(log_string)
          RETURN
       END IF
       a_m = 1d3*0.32d0*(Vm+47.13d0)/(1d0-EXP(-0.1d0*(Vm+47.13d0)))
       b_m = 1d3*0.08d0*EXP(-Vm/11d0)  ![1/sec]

       I_Na = g_Na * m_na**3 * h_na * j_na * (Vm - E_Na)  ![uA/cm^2]

       !! K currents [uA/cm^2]
       !I_Ktof = g_ktof* a_ktof**3 *in_ktof*(Vm-E_K) 
       !I_Ktos = g_ktos* a_ktos**3 *in_ktos*(Vm-E_K) 
       !I_K1   = g_k1*(K_o/(K_o+210.0d0))*((Vm-E_K-1.73d0)/(1d0+EXP(0.0596d0*(Vm-E_K-1.73d0)))) 
          
       I_Ktof = g_ktof* a_ktof *in_ktof*(Vm-E_K) 
       I_Ktos = g_ktos* a_ktos *in_ktos*(Vm-E_K) 
       I_K1   = g_k1*(K_o/(K_o+210.0d0))*((Vm-E_K-1.73d0)/(1d0+EXP(0.0896d0*(Vm-E_K-1.73d0)))) 
       I_Kss  = g_kss*a_kss*in_kss*(Vm-E_K)

       !! Background currents [uA/cm^2]
       I_bNa = g_bna*(Vm-E_Na)
       I_bK  = g_bk*(Vm-E_K)
       I_bCa = g_bca*(Vm-E_Ca)
          
       !! Sarcolemma pumps and Exchangers
       I_NaK  = I_NaK_bar*(1d0/(1d0+0.1245d0*EXP(-0.1d0*(Vm*F_RT))+0.0365d0*(1/7d0*(EXP(Na_o/67300)-1d0))* &
            EXP(-(Vm*F_RT))))*(1d0/(1d0+(21000d0/Na_i)**1.5d0))*(K_o/(K_o+1500d0))
       I_pmca = I_pmca_bar*(Ca_myo**eta_pmca/(Km_pmca**eta_pmca+Ca_myo**eta_pmca))
       !I_ncx  = I_ncx_bar*(Na_i**3 * Ca_o*EXP(eta_ncx*Vm*F_RT) - Na_o**3 *Ca_myo*EXP((eta_ncx-1d0)*Vm*F_RT))/ &
       !     ((Km_ncxNa**3 + Na_o**3)*(Km_ncxCa+Ca_o)*(1d0+ksat_ncx*EXP((eta_ncx-1d0)*Vm*F_RT)))
   
      Cads_total(1) = sum_reduce_double_2host(Ca_ds_dev, NSFU, odata_d)
     I_ncx  = 0.2d0 *I_ncx_bar*(Na_i**3 * Ca_o*EXP(eta_ncx*Vm*F_RT) - Na_o**3 *((Cads_total(1)/NSFU))*EXP((eta_ncx-1d0)*Vm*F_RT))/ &
            ((Km_ncxNa**3 + Na_o**3)*(Km_ncxCa+Ca_o)*(1d0+ksat_ncx*EXP((eta_ncx-1d0)*Vm*F_RT))) + &
      (1.0d0-0.2d0) * I_ncx_bar*(Na_i**3 * Ca_o*EXP(eta_ncx*Vm*F_RT) - Na_o**3 * Ca_myo*EXP((eta_ncx-1d0)*Vm*F_RT))/ &
            ((Km_ncxNa**3 + Na_o**3)*(Km_ncxCa+Ca_o)*(1d0+ksat_ncx*EXP((eta_ncx-1d0)*Vm*F_RT)))


       !I_dhpr_nj_T = SUM(pL_nj, mask=U_Lo_nj .EQ. 1) * P_dhpr_nj *(F2zCa2_RT*Vm) * &
       !     (Ca_myo * EXP(Vm*FzCa_RT) - 0.341d0 * Ca_o)/(EXP(Vm*FzCa_RT)-1d0) ! [uA]
       !TUAN
       !I_dhpr_nj_T = 0.0d0
#endif
          
       !! Updating Fluxes
       !+ membrane fluxes inward Ca2+ fluxes
       !!J_Ca_sarmem = J_bCa - J_pmca - J_ncx + J_dhpr_nj_T
       !!TUANHOANG
       J_Ca_sarmem = (-I_bCa - I_pmca + zCa * I_ncx) * Am_2FVmyo  - I_dhpr_nj_T * inv_2FVmyo  ! [uM/s]
       
          
       !+ internal fluxes
       J_leak = v_leak *(Ca_nsr - Ca_myo)  ![uM/s]
       
       TCa_i2 = (Ca_myo*invKp_myo)**2
       TCa_sr2 = (Ca_nsr*invKp_nsr)**2
       dummy1 = (3.24873d12 * TCa_i2**2 + TCa_i2 * (9.17846d6 - 11478.2d0 * TCa_sr2) - 0.329904d0 * TCa_sr2)
       dummy2 = (0.104217d0 + 17.923d0 * TCa_sr2 + TCa_i2 *(1.75583d6 + 7.61673d6 * TCa_sr2) &
            + TCa_i2**2 *(6.08463d11 + 4.50544d11 * TCa_sr2))  ! D_cycle
       vcycle_pump = dummy1/dummy2  ![1/sec]
       J_serca = cpl_serca*vcycle_pump*Ap  ![uM/sec]


       !!TUANHOANG
       J_ryr_nj_T = +0d0
       !J_ryr_nj_T = SUM(pR_nj, MASK=(U_Ro_nj .EQ. 1 )) * v_ryr_nj_T * (Ca_nsr - Ca_myo) ![uM/s]
       
       !flux forming compound CaB_htrpn in myoplasm
       J_htrpn = kon_htrpn*Ca_myo*(B_htrpnT - CaB_htrpn) - koff_htrpn*CaB_htrpn  ![uM/s]
       J_ltrpn = kon_ltrpn*Ca_myo*(B_ltrpnT - CaB_ltrpn) - koff_ltrpn*CaB_ltrpn 

       !! Euler's==========================================
       m_na = m_na + dt*(a_m*(1-m_na)-b_m*m_na)
       h_na = h_na + dt*(a_h*(1-h_na)-b_h*h_na)
       j_na = j_na + dt*(a_j*(1-j_na)-b_j*j_na)
       a_ktof = a_ktof + dt*(1d3*((0.18064d0*EXP(0.03577d0*(Vm+30d0)))*(1-a_ktof)-(0.3956d0*EXP(-0.06237d0*(Vm+30d0)))*a_ktof))
       
       in_ktof = in_ktof + dt*(1d3*((((0.000152d0*EXP(-(Vm+13.5d0)/7d0))/(0.067083d0*EXP(-(Vm+33.5d0)/7d0)+1d0)))* &
            (1d0-in_ktof)-(((0.00095d0*EXP((Vm+33.5d0)/7d0))/(0.051335d0*EXP((Vm+33.5d0)/7d0)+1d0)))*in_ktof))
       
       a_ktos = a_ktos + dt*(1d3*((1d0/(1d0+EXP(-(Vm+22.5d0)/7.7d0)))-a_ktos)/(0.493d0*EXP(-0.0629d0*Vm)+2.058d0))
       in_ktos = in_ktos + dt*(1d3*((1d0/(1d0+EXP( (Vm+45.2d0)/5.7d0)))-in_ktos)/((270.0d0+1050d0/(1d0+EXP( (Vm+45.2d0)/5.7d0 )))))
       a_kss  = a_kss  + dt*(1d3*((1d0/(1d0+EXP(-(Vm+22.5d0)/7.7d0)))-a_kss)/(0.393d0*EXP(-0.0862d0*Vm)+0.13d0))
       in_kss  = in_kss  + dt*0
       
       !TUAN: 
       !Na_i = Na_i + dt*(-(I_Na + I_bNa + 3*I_ncx + 3*I_NaK)*Am_FVmyo)  ![uM]
       !K_i = K_i + dt*(-(I_K1 + I_Kss + I_Ktof + I_Ktos + I_bK - 2*I_NaK)*Am_FVmyo)  ![uM]
       
       I_dhpr_total(1) = sum_reduce_double_2host(I_dhpr_dev, NSFU, odata_d)
       
#if _RUNMODE_ == _MODE_I_CLAMP
       !I-clamp (dVm/dt)
       Vm = Vm + dt*(-1.0d3/Csc * (I_Na + (I_dhpr_total(1) + I_dhpr_nj_T)/Am + I_K1 + I_Kss + I_Ktof + &
            I_Ktos + I_ncx  + I_NaK + I_pmca + I_bNa + I_bCa + I_bk + I_app))
       
#endif
       !! synchronize GPU data
       cuErr = cudaDeviceSynchronize()
          
#if _SAFE_MODE_ == _YES
       IF (cuErr .NE. cudaSuccess) THEN
          WRITE(log_string,*) "Error: kernel"
          CALL write2log(log_string)
          WRITE(log_string,*) cudaGetErrorString(cudaGetLastError())
          CALL write2log(log_string)
          STOP
       END IF
#endif

       J_refill_total(1) = sum_reduce_double_2host(J_refill_dev, NSFU, odata_d)
       
       Ca_nsr = Ca_nsr + dt*((beta_nsr*inv_lambda_nsr)  * (J_serca - J_leak - J_refill_total(1) - J_ryr_nj_T - J_ryr_orphan))
          
       !!beta_myo = 1.0d0 /(1.0d0 + (B_myoT * Km_myo)/((Km_myo + Ca_myo)**2))
       dummy1 = (Km_myo + Ca_myo)**2
       dummy2 = dummy1 +  B_myoT * Km_myo
       beta_myo = dummy1 / dummy2
       
       J_efflux_total(1) = sum_reduce_double_2host(J_efflux_dev, NSFU, odata_d)
       
       Ca_myo = Ca_myo + dt*(beta_myo * (J_efflux_total(1) + J_leak + J_Ca_sarmem + J_ryr_nj_T - J_serca  - J_htrpn - J_ltrpn + J_ryr_orphan))

       !!TUANHOANG
       !PRINT *, J_leak, J_ryr_nj_T, time
       
       CaB_htrpn = CaB_htrpn + dt*J_htrpn
       CaB_ltrpn = CaB_ltrpn + dt*J_ltrpn
       
       !tmp_pL_nj = pL_nj
       !DO ii=1, irow_L_nj
          !pL_nj(ii) = tmp_pL_nj(ii) + dt*SUM(tmp_pL_nj(:)*(kva*aKv1_nj(:,ii) + kvi*aKv2_nj(:,ii) + aKm1_nj(:,ii) + &
          !     Ca_myo**eta_LCC * aKp1_nj(:,ii)))
       !END DO
       !tmp_pR_nj = pR_nj
       !DO ii=1, irow_R_nj
       !  pR_nj(ii) = tmp_pR_nj(ii) + dt*SUM(tmp_pR_nj(:)*(aKm2_nj(:,ii) + (k_jsr0*Ca_nsr+k_jsr1)*Ca_myo**eta_RyR * aKp2_nj(:,ii)))
       !ENDDO

       CaF  = CaF  + dt*(-koff_F*CaF + kon_F*Ca_myo*(F_T - CaF));
       
       !WRITE DATA OUT
       IF (time .GE. (t_end - duration_write_lastbeats)) THEN 
          IF ((time - time_prev1) .GE. write_interval) THEN 
             CALL data_io_write()
             time_prev1 = time
          END IF
       ENDIF
       !WRITE DATA OUT (First-Beat) every 100 time-steps
       IF (write_firstbeat .NE. 0) THEN
          IF (time .LE. (t_start+1.0d0)) THEN
             !! write out the first beat 
             IF ((time - time_prev_fb) .GE. write_interval) THEN 
                CALL data_io_write_fb()
                time_prev_fb = time
             END IF
          ELSEIF (is_done_firstbeat .EQ. .FALSE.) THEN
             CALL data_io_write_fb_last()
             CALL spark_summarize_fb()
             is_done_firstbeat = .TRUE.
          ENDIF
       ENDIF
    END DO INNERLOOP
!!!TUAN: switch back if want to output every 100 time steps
!!$       call data_io_write()
    IF (write_firstbeat .NE. 0) THEN
       IF (time .LE. (t_start+1.0d0)) THEN
          !! spark analysis for the first beat 
#if _RUNMODE_ == _MODE_RESTING || _RUNMODE_ == _MODE_CLOSEDCELL
          CALL analyze_Sparks <<<gridSize, blocksize>>> (REAL(time), REAL(Ca_nsr))
#endif
       ENDIF
    ENDIF
    
    ! print out data analysis for resting simulation
    IF (time .GE. time_start_spark_analysis) THEN
#if _RUNMODE_ == _MODE_RESTING || _RUNMODE_ == _MODE_CLOSEDCELL
       ! for resting or closed-cell, only do spark analysis
       CALL analyze_Sparks <<<gridSize, blocksize>>> (REAL(time), REAL(Ca_nsr))
#endif
       ! APD analysis + spark analysis in 1Hz stimulation of diastolic phase
       J_ryr_total(1) = sum_reduce_double_2host(J_ryr_dev, NSFU, odata_d)
       
       IF (NSFU_normal .LT. NSFU)  J_ryr_orphan_total(1) = &
            sum_reduce_double_2host(J_ryr_dev(NSFU_normal+1:NSFU), NSFU_orphan, odata_d)

       J_dhpr_total = -I_dhpr_total(1) * inv_2FVmyo

       accumulate_J_ryr = accumulate_J_ryr + (J_ryr_total(1)+J_leak+J_ryr_nj_T) * (time - time_prev2)
       accumulate_J_dhpr = accumulate_J_dhpr + J_dhpr_total * (time - time_prev2)
       accumulate_J_ryr_orphan = accumulate_J_ryr_orphan + J_ryr_orphan_total(1) * (time - time_prev2)
       time_prev2 = time
       
       DO ii = 1, NUM_APD_TYPES
          IF (Vm .GE. Vmthreshold(ii) .AND. is_inside_AP(ii) .EQ. .FALSE.) THEN
             is_inside_AP(ii) = .TRUE.
             APstart(ii) = time
          END IF
          IF (Vm .LE. Vmthreshold(ii) .AND. is_inside_AP(ii) .EQ. .TRUE.) THEN
             is_inside_AP(ii) = .FALSE.
             APend(ii) = time
             IF (APstart(ii) .LT. APend(ii)) THEN
                APD(ii) = APD(ii) + (APend(ii) - APstart(ii))
                APcount(ii) = APcount(ii) + 1
             END IF
             IF (ii .EQ. 1) THEN
                is_in_DiastolicPhase = .TRUE.
                time_start_DiastolicPhase = time
             END IF
          END IF
       END DO
       IF (Hz .LE. 1.0d0) THEN
          ! only do diastolic spark analysis when Hz .le. 1
          IF (is_in_DiastolicPhase) THEN
             !end an AP, can measure diastolic spark statistics
             !until the end of that second. This works when Hz < 1.
             CALL analyze_Sparks <<<gridSize, blocksize>>> (REAL(time), REAL(Ca_nsr))
             IF (time .GT. CEILING(time_start_DiastolicPhase)) is_in_DiastolicPhase = .FALSE.
          END IF
       ENDIF
    ENDIF
    
  END FUNCTION  simulate_NextSteps

  
  !!>****f* tscell_gpu_utility/get_compPdt_3
  !! SYNOPSIS  
  !!    [Device] Get the transition rate at which the SFU doesn't
  !!    change state. This will help determining the appropriate time
  !!    step so that 90% of the CRUs doesn't change the state. 
  !!   
  !! INPUTS
  !!   kva, kvi :: 2 parameters associated with transition rates formula
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  attributes(global) SUBROUTINE get_compPdt_2(kva, kvi, bigRos)
    USE cudafor
    IMPLICIT NONE
    
    REAL(KIND=DP), VALUE:: kva, kvi,bigRos
    ! local
    REAL(KIND=dp) :: bigC1, bigC2, chiC_mode1,chiC_mode2, chiO_mode1,chiO_mode2
    REAL(KIND=dp) :: phi, dtemp, dtemp2
    INTEGER:: tx, ii !,  ct   ! ct = current state
    REAL(KIND=dp), DIMENSION(BLOCKSIZE_1D), shared :: s_cads, s_ryrgate,s_Nclose_RyRs_mode1, &
                         s_Nclose_RyRs_mode2, s_Nopen_RyRs_mode1, s_Nopen_RyRs_mode2
    INTEGER,VOLATILE :: lcc_ct, ryr_ct
    
    ii =  (blockIdx%x-1) * blockDim%x + threadIdx%x
    
    IF (ii .LE. int_ca(1)) THEN
       tx = threadIdx%x
       s_cads(tx) = Ca_ds_dev(ii)
       s_ryrgate(tx) = RYRgate_dev(ii)
       s_Nclose_RyRs_mode1(tx) = Nclose_RyRs_mode1_dev(ii)
       s_Nclose_RyRs_mode2(tx) = Nclose_RyRs_mode2_dev(ii)
       s_Nopen_RyRs_mode1(tx) = Nopen_RyRs_mode1_dev(ii)
       s_Nopen_RyRs_mode2(tx) = Nopen_RyRs_mode2_dev(ii)
    ENDIF
    CALL syncthreads()
       
    IF (ii .LE. int_ca(1)) THEN
       CALL state_breakdown_ker(CRU_states_dev(ii), int_ca(10), lcc_ct, ryr_ct)
       !bigC2 = Ca_ds_dev(ii)**eta_RyR
       !bigC1 = Ca_ds_dev(ii)**eta_LCC
       dtemp = LOG(s_cads(tx))
       bigC2 = EXP(dtemp*dp_ca(14))
       bigC1 = EXP(dtemp*dp_ca(15))
       
       !chiC = EXP(-Ej*0.5d0*((N_R-RYRgate_dev(ii))*Ecc - (RYRgate_dev(ii)-1d0)*Eoo))
       dtemp = (s_Nclose_RyRs_mode1(tx)+s_Nclose_RyRs_mode2(tx))*dp_ca(2)
       dtemp2 = dtemp - (s_Nopen_RyRs_mode1(tx)+s_Nopen_RyRs_mode2(tx) - 1.d0)*dp_ca(3)
       dtemp = dp_ca(1)*dtemp2
       chiC_mode1 = EXP(dtemp)

       dtemp = (s_Nclose_RyRs_mode2(tx)+s_Nclose_RyRs_mode1(tx))*dp_ca(2)
       dtemp2 = dtemp - (s_Nopen_RyRs_mode2(tx) +s_Nopen_RyRs_mode1(tx)- 1.d0)*dp_ca(3)
       dtemp = dp_ca(1)*dtemp2
       chiC_mode2 = EXP(dtemp)
       
       !chiO = EXP(-Ej*0.5d0*(RYRgate_dev(ii)*Eoo - (N_R-RYRgate_dev(ii)- 1d0)*Ecc))
       dtemp = (s_Nopen_RyRs_mode1(tx)+s_Nopen_RyRs_mode2(tx))*dp_ca(3)
       dtemp2 = dtemp - (s_Nclose_RyRs_mode1(tx) + s_Nclose_RyRs_mode2(tx)- 1.d0)*dp_ca(2)
       dtemp = dp_ca(1)*dtemp2
       chiO_mode1 = EXP(dtemp)
       
       dtemp = (s_Nopen_RyRs_mode2(tx)+s_Nopen_RyRs_mode1(tx))*dp_ca(3)
       dtemp2 = dtemp - (s_Nclose_RyRs_mode2(tx) + s_Nclose_RyRs_mode1(tx)- 1.d0)*dp_ca(2)
       dtemp = dp_ca(1)*dtemp2
       chiO_mode2 = EXP(dtemp)

       !phi = k_jsr0*Ca_jsr_dev(ii) + k_jsr1
       phi = dp_ca(4)*Ca_jsr_dev(ii) + dp_ca(5)
       !Rostemp = RosDS**eta_Ros

       ! use 1st column for negative value
       dtemp = chiC_mode1 * ryr_compKfree_mode1_dev(ryr_ct,1)  
       dtemp2 = dtemp + bigC2 * chiO_mode1 * phi * ryr_compKcal_act_mode1_dev(ryr_ct,1)
       dtemp = dtemp2 + chiC_mode2*ryr_compKfree_mode2_dev(ryr_ct,1)
       dtemp2 = dtemp + bigC2 * chiO_mode2 * phi * ryr_compKcal_act_mode2_dev(ryr_ct,1)
       dtemp = dtemp2 + bigRos * ryr_compKRos_dev(ryr_ct,1)
       dtemp2 = dtemp + ryr_compKfree_dev(ryr_ct,1)
       dtemp = dtemp2 + kvi * lcc_compKv_inact_dev(lcc_ct,1) 
       dtemp2 = dtemp + kva * lcc_compKv_act_dev(lcc_ct,1)
       dtemp = dtemp2 + bigC1 * lcc_compKcal_inact_dev(lcc_ct,1) 
       compPdt_dev(ii) = dtemp + lcc_compKfree_dev(lcc_ct,1) 
    END IF
  END SUBROUTINE get_compPdt_2

  !=============================================================================
  ! Update CRU's state + calcium concentration (Ca-ds, Ca-jsr)
  !=============================================================================
  attributes(global) SUBROUTINE updateSystem_2(Xrnd, dt, invdt, &
       Ca_nsr, Ca_myo, Vm, kva, kvi, dp_arg1, dp_arg2, bigRos) 
    !new method
    USE cudafor
    IMPLICIT NONE
    REAL(KIND=dp), DIMENSION(:), device :: Xrnd
    REAL(KIND=dp), VALUE :: dt,  invdt, Ca_nsr, Ca_myo, Vm, kva, kvi, dp_arg1, dp_arg2, bigRos

    !local
    INTEGER :: ii, jj, kk, tx
    REAL(KIND=dp) :: dtemp, dtemp2, compP, J_ryr
    !INTEGER :: nt, ct, cnt
    !nt = new state, ct = current state, cnt= candidate next state
    REAL(KIND=dp) :: bigC1, bigC2, chiC_mode1, chiC_mode2, chiO_mode1, chiO_mode2, phi, J_dhpr
    REAL(KIND=dp) :: nryr ! new RYRgate
    REAL(KIND=dp) :: ndhpr !new dhprgate


    INTEGER :: itemp
    REAL(KIND=dp), DIMENSION(BLOCKSIZE_1D), shared :: s_cads, s_ryrgate, s_Nclose_RyRs_mode1, & 
                          s_Nclose_RyRs_mode2, s_Nopen_RyRs_mode1, s_Nopen_RyRs_mode2
    INTEGER, VOLATILE :: lcc_ct, ryr_ct
    INTEGER :: lcc_nt, ryr_nt
    INTEGER :: is_examining
    LOGICAL :: found_nt
    INTEGER :: candidate_new_state !can be any cluster types (LCC or RyR)
    INTEGER, PARAMETER :: EXAMINING_LCC = 1, EXAMINING_RYR = 2

    ! thread Index
    ii =  (blockIdx%x-1) * blockDim%x + threadIdx%x

    IF (ii .LE. int_ca(1)) THEN
       tx = threadIdx%x
       s_cads(tx) = Ca_ds_dev(ii)
       s_ryrgate(tx) = RYRgate_dev(ii)
       s_Nclose_RyRs_mode1(tx) = Nclose_RyRs_mode1_dev(ii)
       s_Nclose_RyRs_mode2(tx) = Nclose_RyRs_mode2_dev(ii)
       s_Nopen_RyRs_mode1(tx) = Nopen_RyRs_mode1_dev(ii)
       s_Nopen_RyRs_mode2(tx) = Nopen_RyRs_mode2_dev(ii)
    ENDIF
    CALL syncthreads()

    IF (ii .LE. int_ca(1)) THEN
       CALL state_breakdown_ker(CRU_states_dev(ii), int_ca(10), lcc_ct, ryr_ct)

       lcc_nt = lcc_ct
       ryr_nt = ryr_ct
       found_nt = .FALSE.
       compP = (compPdt_dev(ii) + invdt) * dt

       IF (compP .GE. Xrnd(ii)) THEN ! channel doesn't change state (easy)
          found_nt = .TRUE.
       ELSE ! channel change to some new state
          dtemp = LOG(s_cads(tx))
          dtemp2 = dtemp *dp_ca(14)
          bigC2 = EXP(dtemp2)
          dtemp2 = dtemp*dp_ca(15)
          bigC1 = EXP(dtemp2)

          !!dtemp = (int_ca(2)-s_ryrgate(tx))*dp_ca(2)
          !!dtemp2 = dtemp - (s_ryrgate(tx) - 1.d0)*dp_ca(3)
          !!dtemp = dp_ca(1)*dtemp2
          !!chiC = EXP(dtemp)
          
          !chiC = EXP(-Ej*0.5d0*((N_R-RYRgate_dev(ii))*Ecc - (RYRgate_dev(ii)-1d0)*Eoo))
          dtemp = (s_Nclose_RyRs_mode1(tx)+s_Nclose_RyRs_mode2(tx))*dp_ca(2)
          dtemp2 = dtemp - (s_Nopen_RyRs_mode1(tx) +s_Nopen_RyRs_mode2(tx)- 1.d0)*dp_ca(3)
          dtemp = dp_ca(1)*dtemp2
          chiC_mode1 = EXP(dtemp)

          dtemp = (s_Nclose_RyRs_mode2(tx)+ s_Nclose_RyRs_mode1(tx))*dp_ca(2)
          dtemp2 = dtemp - (s_Nopen_RyRs_mode2(tx) +s_Nopen_RyRs_mode1(tx)- 1.d0)*dp_ca(3)
          dtemp = dp_ca(1)*dtemp2
          chiC_mode2 = EXP(dtemp)

          !dtemp = (int_ca(2)-RYRgate_dev(ii))*dp_ca(2)
          !dtemp2 = dtemp - (RYRgate_dev(ii) - 1d0)*dp_ca(3)
          !dtemp = dp_ca(1)*dtemp2
          !chiC = EXP(dtemp)

          !!dtemp = s_ryrgate(tx)*dp_ca(3)
          !!dtemp2 = dtemp - (int_ca(2) - s_ryrgate(tx) - 1d0)*dp_ca(2)
          !!dtemp = dp_ca(1)*dtemp2
          !!chiO = EXP(dtemp)

          !chiO = EXP(-Ej*0.5d0*(RYRgate_dev(ii)*Eoo - (N_R-RYRgate_dev(ii)- 1d0)*Ecc))          
          dtemp = (s_Nopen_RyRs_mode1(tx)+s_Nopen_RyRs_mode2(tx))*dp_ca(3)
          dtemp2 = dtemp - (s_Nclose_RyRs_mode1(tx) + s_Nclose_RyRs_mode2(tx)- 1.d0)*dp_ca(2)
          dtemp = dp_ca(1)*dtemp2
          chiO_mode1 = EXP(dtemp)
       
          dtemp = (s_Nopen_RyRs_mode2(tx)+s_Nopen_RyRs_mode1(tx))*dp_ca(3)
          dtemp2 = dtemp - (s_Nclose_RyRs_mode2(tx) +s_Nclose_RyRs_mode1(tx) - 1.d0)*dp_ca(2)
          dtemp = dp_ca(1)*dtemp2
          chiO_mode2 = EXP(dtemp)
          !dtemp = RYRgate_dev(ii)*dp_ca(3)
          !dtemp2= dtemp - (int_ca(2)-RYRgate_dev(ii)- 1d0)*dp_ca(2)
          !dtemp = dp_ca(1)*dtemp2
          !chiO = EXP(dtemp)

          phi = dp_ca(4)*Ca_jsr_dev(ii) + dp_ca(5)
          !RosDS = 1.0d0
          DO jj=2, int_ca(15)
             IF (found_nt .EQ. .FALSE.) THEN
                IF (jj .LE. int_ca(3)) THEN
                   kk = jj
                   dtemp = bigC2 * chiO_mode1 * phi * ryr_compKcal_act_mode1_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKcal_act_mode1_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)-1) THEN
                   kk = jj-(int_ca(3)-1)
                   dtemp = chiC_mode1 * ryr_compKfree_mode1_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKfree_mode1_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)-2) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)-2)
                   dtemp = bigC2*chiO_mode2 *phi* ryr_compKcal_act_mode2_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKcal_act_mode2_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)-3) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)-3)
                   dtemp = chiC_mode2 * ryr_compKfree_mode2_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKfree_mode2_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)-4) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)-4)
                   dtemp = ryr_compKfree_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKfree_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)+int_ca(14)-5) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)+int_ca(13)-5)
                   dtemp = bigRos*ryr_compKRos_dev(ryr_ct,kk)
                   candidate_new_state = ryr_indxKRos_dev(ryr_ct,kk)
                   is_examining = EXAMINING_RYR
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)+int_ca(14)+int_ca(5)-6) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)+int_ca(13)+int_ca(14)-6)
                   dtemp = kva*lcc_compKv_act_dev(lcc_ct,kk)
                   candidate_new_state = lcc_indxKv_act_dev(lcc_ct,kk)
                   is_examining = EXAMINING_LCC

                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)+int_ca(14)+int_ca(5)+int_ca(6)-7) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)+int_ca(13)+int_ca(14)+int_ca(5)-7)
                   dtemp = kvi*lcc_compKv_inact_dev(lcc_ct,kk)
                   candidate_new_state = lcc_indxKv_inact_dev(lcc_ct,kk)
                   is_examining = EXAMINING_LCC
                   
                ELSE IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)+int_ca(14)+int_ca(5)+int_ca(6)+int_ca(7)-8) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)+int_ca(13)+int_ca(14)+int_ca(5)+int_ca(6)-8)
                   dtemp = bigC1 * lcc_compKcal_inact_dev(lcc_ct,kk)
                   candidate_new_state = lcc_indxKcal_inact_dev(lcc_ct,kk)
                   is_examining = EXAMINING_LCC                   
                   
                ELSE !IF (jj .LE. int_ca(3)+int_ca(4)+ int_ca(11)+int_ca(12)+ int_ca(13)+int_ca(14)+int_ca(5)+int_ca(6)+int_ca(7)+int_ca(8)-9) THEN
                   kk = jj-(int_ca(3)+ int_ca(4)+ int_ca(11)+int_ca(12)+int_ca(13)+int_ca(14)+int_ca(5)+int_ca(6)+int_ca(7)-9)
                   dtemp = lcc_compKfree_dev(lcc_ct,kk)
                   candidate_new_state = lcc_indxKfree_dev(lcc_ct,kk)
                   is_examining = EXAMINING_LCC                   
                   
                   !if we need more else, need to retain IF in the previous one
                END IF
                compP = compP + dtemp * dt 
                IF (compP .GE. Xrnd(ii)) THEN
                   found_nt = .TRUE.
                   SELECT CASE (is_examining)
                   CASE(EXAMINING_LCC)
                      !to avoid LCC gating at negative potential
                      !IF (Vm .GE. -70.0d0)  
                      lcc_nt = candidate_new_state
                   CASE(EXAMINING_RYR)
                      ryr_nt = candidate_new_state
                   END SELECT
                END IF
             END IF
          END DO
       END IF

       nryr = U_Ro_dev(ryr_nt)
       ndhpr = U_Lo_dev(lcc_nt)

       J_efflux_dev(ii) = dp_ca(11)*(s_cads(tx)-Ca_myo)
       J_refill_dev(ii) = dp_ca(10)*(Ca_nsr-Ca_jsr_dev(ii))

#if (_RUNMODE_ == _MODE_I_CLAMP || _RUNMODE_ == _MODE_V_CLAMP)
       !I-clamp or V-clamp
       I_dhpr_dev(ii)  = ndhpr * dp_ca(13) * dp_ca(16) * Vm * &
            (s_cads(tx) * dp_arg1 - dp_ca(18)) * dp_arg2  ! [uA]
       !J_dhpr_dev(ii)  = -I_dhpr_dev(ii) * dp_ca(17)  ! [uA/s]
       J_dhpr  = -I_dhpr_dev(ii) * dp_ca(17)  ! [uA/s]
#else
       !J_dhpr_dev(ii) = 0d0
       J_dhpr = 0d0
#endif

       !!J_ryr = nryr * dp_ca(9) * (Ca_jsr_dev(ii)-Ca_ds_dev(ii)) !J_ryr
       J_ryr = nryr * dp_ca(9) * (Ca_jsr_dev(ii)-s_cads(tx))
       J_ryr_dev(ii) = J_ryr
#if _USE_CaDS_APPROX_ == _NO
       !exact formula
       !s_cads(tx) = s_cads(tx) + dt * dp_ca(6) * dp_ca(7) * (J_ryr - J_efflux_dev(ii) + J_dhpr_dev(ii))
       s_cads(tx) = s_cads(tx) + dt * dp_ca(6) * dp_ca(7) * (J_ryr - J_efflux_dev(ii) + J_dhpr)
#else
       dtemp = ndhpr * dp_ca(13) * dp_ca(16) * Vm * dp_arg2 * dp_ca(17)
       dtemp2 = nryr * dp_ca(9)
       s_cads(tx) = (dtemp2 * Ca_jsr_dev(ii) + dp_ca(11) * Ca_myo + dtemp * dp_ca(18))/&
            (dtemp2 + dp_ca(11) + dtemp)
       !dtemp = (dp_ca(13)*dp_ca(16) * Vm * dp_ca(18) * dp_arg2) * &
       !     dp_ca(17) * ndhpr + dp_ca(11) * Ca_myo + nryr * dp_ca(9) * Ca_jsr_dev(ii)
       !dtemp2 = nryr * dp_ca(9) + dp_ca(11) + (dp_ca(13) * dp_ca(16) * Vm * dp_arg1 * dp_arg2)  * dp_ca(17) * ndhpr
       !s_cads(tx) = dtemp / dtemp2       
#endif
       dtemp = dp_ca(19) + Ca_jsr_dev(ii)
       dtemp2 = dtemp * dtemp
       dtemp = dtemp2 / (dp_ca(19)*dp_ca(20) + dtemp2)   !beta_jsr
       Ca_jsr_dev(ii) = Ca_jsr_dev(ii) + dt*(dtemp * dp_ca(8)) * (J_refill_dev(ii) - J_ryr)

       CRU_states_dev(ii) =  state_combine_single_ker(lcc_nt, ryr_nt, int_ca(10))
       !CRU_states_dev(ii) = (lcc_nt-1) * int_ca(10) + ryr_nt
       RYRgate_dev(ii) = nryr
       DHPRgate_dev(ii) = ndhpr
       Ca_ds_dev(ii) = s_cads(tx)
       Nclose_RyRs_mode1_dev(ii) = U_RyR_State_dev(ryr_nt,1)
       Nclose_RyRs_mode2_dev(ii) = U_RyR_State_dev(ryr_nt,4)
       Nopen_RyRs_mode1_dev(ii) = U_RyR_State_dev(ryr_nt,2)
       Nopen_RyRs_mode2_dev(ii) = U_RyR_State_dev(ryr_nt,3)
       

    END IF
  END SUBROUTINE updateSystem_2

  attributes(global) SUBROUTINE analyze_Sparks(time, Ca_nsr)
    USE cudafor
    IMPLICIT NONE
    REAL, VALUE :: time
    REAL, VALUE :: Ca_nsr

    INTEGER :: ii
    REAL :: dtemp1
    
    ! thread Index
    ii =  (blockIdx%x-1) * blockDim%x + threadIdx%x

    IF (ii .LE. int_ca(1)) THEN
       IF (Ryrgate_dev(ii) .GE. 1 .AND. flag1(ii) .EQ. .FALSE.) THEN
          ! spark/blink/quark start
          flag1(ii) = .TRUE.; flag2(ii) = .TRUE.; flag3(ii) = .TRUE.;
          ts_spark(ii) = time; ! time start spark/quark
          ts_blink(ii) = time; ! time start blink
       ELSEIF (RyRgate_dev(ii) .EQ. 0 .AND. flag1(ii) .EQ. .TRUE.  & 
            .AND. flag4(ii) .EQ. .FALSE.) THEN
          ! detect quarks
          ! count # of quarks with 1,2,3... opening RyRs
          Nquarks_dev(ii)= Nquarks_dev(ii) + 1 
          ! reset flags
          flag1(ii) = .FALSE.; flag2(ii) = .FALSE.; flag3(ii) = .FALSE.;
       !ELSEIF  (RyRgate_dev(ii) .GE. float(int_ca(2))*0.909d0 .AND.  &
       ELSEIF  (RyRgate_dev(ii) .GE. float(int_ca(2)) .AND.  &
            flag1(ii) .EQ. .TRUE. .AND. flag3(ii) .EQ. .TRUE.) THEN
          ! mature spark
!!$          IF (RyRgate_dev(ii) .LT. prev_RyRgate(ii)) THEN
!!$             ! spark peak (increase spark count)
!!$             Nsparks_dev(ii) = Nsparks_dev(ii) + 1
!!$             SparkPeaks(ii, Nsparks_dev(ii)) = prev_Ca_ds(ii)
!!$             flag3(ii) = .FALSE.; flag4(ii) = .TRUE.;
!!$          ELSE
!!$             prev_Ca_ds(ii) = Ca_ds_dev(ii)
!!$          END IF
          flag3(ii) = .FALSE.
          flag4(ii) = .TRUE.
          Nsparks_dev(ii) = Nsparks_dev(ii) + 1
          tot_SparkPeaks(ii) = tot_SparkPeaks(ii) + Ca_ds_dev(ii)  
       ELSE
          IF (RyRgate_dev(ii) .EQ. 0 .AND. flag1(ii) .EQ. .TRUE. .AND. &
            flag4(ii) .EQ. .TRUE.) THEN
             ! spark termination 
             flag1(ii) = .FALSE.; 
             tot_SparkDurations(ii) = tot_SparkDurations(ii) + (time-ts_spark(ii))
          ENDIF
          IF  (Ca_jsr_dev(ii) .GT. 0.99 * Ca_nsr .AND. flag2(ii) .EQ. .TRUE. .AND. &
               flag4(ii) .EQ. .TRUE.)  THEN
             ! blink detection
             flag4(ii) = .FALSE.; flag2(ii) = .FALSE.
             Nblinks_dev(ii) = Nblinks_dev(ii) + 1;
             tot_BlinkDurations(ii) = tot_BlinkDurations(ii) + (time - ts_blink(ii))
          ENDIF
       ENDIF
       !prev_RyRgate(ii) = RyRgate_dev(ii)
    END IF
  END SUBROUTINE analyze_Sparks


  SUBROUTINE free_device()
    CALL cleanup_gpudata_ecc()
    DEALLOCATE(Ca_ds_dev, Ca_jsr_dev)
  END SUBROUTINE free_device

  SUBROUTINE cleanup_gpudata_ecc()
    IF (ALLOCATED(lcc_compKv_act_dev)) DEALLOCATE(lcc_compKv_act_dev)
    IF (ALLOCATED(lcc_compKv_inact_dev)) DEALLOCATE(lcc_compKv_inact_dev)
    IF (ALLOCATED(lcc_compKfree_dev)) DEALLOCATE(lcc_compKfree_dev)
    IF (ALLOCATED(ryr_compKfree_dev)) DEALLOCATE(ryr_compKfree_dev)
    IF (ALLOCATED(ryr_compKRos_dev)) DEALLOCATE(ryr_compKRos_dev)
    IF (ALLOCATED(ryr_compKfree_mode1_dev)) DEALLOCATE(ryr_compKfree_mode1_dev)
    IF (ALLOCATED(ryr_compKfree_mode2_dev)) DEALLOCATE(ryr_compKfree_mode2_dev)
    IF (ALLOCATED(lcc_compKcal_inact_dev)) DEALLOCATE(lcc_compKcal_inact_dev)
    IF (ALLOCATED(ryr_compKcal_act_mode1_dev)) DEALLOCATE(ryr_compKcal_act_mode1_dev)
    IF (ALLOCATED(ryr_compKcal_act_mode2_dev)) DEALLOCATE(ryr_compKcal_act_mode2_dev)
    
    IF (ALLOCATED(lcc_indxKv_act_dev)) DEALLOCATE(lcc_indxKv_act_dev)
    IF (ALLOCATED(lcc_indxKv_inact_dev)) DEALLOCATE(lcc_indxKv_inact_dev)
    IF (ALLOCATED(lcc_indxKfree_dev)) DEALLOCATE(lcc_indxKfree_dev)
    IF (ALLOCATED(ryr_indxKfree_dev)) DEALLOCATE(ryr_indxKfree_dev)
    IF (ALLOCATED(ryr_indxKRos_dev)) DEALLOCATE(ryr_indxKRos_dev)
    IF (ALLOCATED(ryr_indxKfree_mode1_dev)) DEALLOCATE(ryr_indxKfree_mode1_dev)
    IF (ALLOCATED(ryr_indxKfree_mode2_dev)) DEALLOCATE(ryr_indxKfree_mode2_dev)    
    IF (ALLOCATED(lcc_indxKcal_inact_dev)) DEALLOCATE(lcc_indxKcal_inact_dev)
    IF (ALLOCATED(ryr_indxKcal_act_mode1_dev)) DEALLOCATE(ryr_indxKcal_act_mode1_dev)
    IF (ALLOCATED(ryr_indxKcal_act_mode2_dev)) DEALLOCATE(ryr_indxKcal_act_mode2_dev)

    IF (ALLOCATED(compPdt_dev)) DEALLOCATE(compPdt_dev)
    IF (ALLOCATED(U_Ro_dev)) DEALLOCATE(U_Ro_dev)
    IF (ALLOCATED(U_Lo_dev)) DEALLOCATE(U_Lo_dev)
    IF (ALLOCATED(U_RyR_State_dev)) DEALLOCATE(U_RyR_State_dev)
    IF (ALLOCATED(CRU_states_dev)) DEALLOCATE(CRU_states_dev)
    IF (ALLOCATED(ryrgate_dev)) DEALLOCATE(ryrgate_dev)
    IF (ALLOCATED(dhprgate_dev)) DEALLOCATE(dhprgate_dev)
    IF (ALLOCATED(J_efflux_dev)) DEALLOCATE(J_efflux_dev)
    IF (ALLOCATED(J_refill_dev)) DEALLOCATE(J_refill_dev)
    IF (ALLOCATED(J_efflux_total)) DEALLOCATE(J_efflux_total)
    IF (ALLOCATED(J_refill_total)) DEALLOCATE(J_refill_total)
    IF (ALLOCATED(I_dhpr_total)) DEALLOCATE(I_dhpr_total)
    IF (ALLOCATED(I_dhpr_dev)) DEALLOCATE(I_dhpr_dev)
    IF (ALLOCATED(Cads_total)) DEALLOCATE(Cads_total)
  END SUBROUTINE cleanup_gpudata_ecc

  SUBROUTINE data_io_init_fb()
    LOGICAL :: is_exist
    INTEGER :: ii

#if _DATAMODE_ == _IO_TEST
    CALL datamode1_init_fb()
#endif
  END SUBROUTINE

  SUBROUTINE data_io_init()
    LOGICAL :: is_exist
    INTEGER :: ii

#if _DATAMODE_ == _IO_TEST
    CALL datamode1_init()
#elif _DATAMODE_ == _IO_2_HDF
    CALL datamode2_init()
#elif _DATAMODE_ == _IO_2_ASCII
    CALL datamode3_init()
#elif _DATAMODE_ == _IO_2_SILO_3D
    CALL datamode4_init()
#elif _DATAMODE_ == _IO_2_SILO_SLICE
    CALL datamode5_init()
#elif _DATAMODE_ == _IO_NONE
    WRITE(log_string,*) "No I/O"
    CALL write2log(log_string)
#elif _DATAMODE_ == _IO_2_HDF2
    CALL datamode7_init()
#endif  !DATAMODE
    
  END SUBROUTINE data_io_init

#if _DATAMODE_ == _IO_TEST
  SUBROUTINE datamode1_init_fb()

    sys_command = "mkdir " // TRIM(where2save) // "/firstbeat"
    CALL system(sys_command)
    fstring = TRIM(where2save) // "/firstbeat/ca_ds.dat"
    OPEN(unit=11, file=TRIM(fstring))
    OPEN(unit=12, file=TRIM(where2save) //"/firstbeat/ca_jsr.dat")
    !OPEN(unit=18, file=TRIM(where2save) //"/firstbeat/state_Indx.dat")
  END SUBROUTINE datamode1_init_fb
  SUBROUTINE datamode1_print_fb(done_fb)
    integer, optional :: done_fb
    INTEGER :: jj
    CHARACTER(100) :: formatstring

    WRITE(formatstring, "(i0)") NSFU+1
    formatstring = "(" // TRIM(formatstring) //  "(2x,1pe12.5))"

       Ca_ds = Ca_ds_dev
       Ca_jsr = Ca_jsr_dev
       RYRgate = Ryrgate_dev
       DHPRgate = DHPRgate_dev
       CRU_states = CRU_states_dev
       CALL state_breakdown_all(CRU_states, irow_R, LCCcluster_states, RYRcluster_states)
       WRITE(11,TRIM(formatstring)) time,(Ca_ds(jj), jj=1,NSFU)
       WRITE(12,TRIM(formatstring)) time,(Ca_jsr(jj), jj=1,NSFU)
       !WRITE(18,TRIM(formatstring)) time,(real(RYRcluster_states(jj)), jj=1,NSFU)

   if(present(done_fb)) then 
     if (done_fb .eq. 1) then 
         CLOSE(11) 
         CLOSE(12)
         !CLOSE(18)
         call datamode1_init()
     endif
    endif

  END SUBROUTINE datamode1_print_fb

  SUBROUTINE datamode1_init()
      
    WRITE(log_string,*) "Mode 1"
    CALL write2log(log_string)
    
    fstring = TRIM(where2save) // "/ca_ds.dat"
    !fstring = TRIM(where2save) // "/CaB_htrpn.dat"
    OPEN(unit=11, file=TRIM(fstring))
    !WRITE(log_string,*) "CaB_ltrpn.dat: calcium bound to high-affinity troponin "
    WRITE(log_string,"(A,A)") "calcium diadic subspace: ", TRIM(fstring)
    CALL write2log(log_string)

    OPEN(unit=12, file=TRIM(where2save) // "/ca_jsr.dat")
    !WRITE(log_string,*) "CaB_ltrpn.dat: calcium bound to low-affinity troponin"
    WRITE(log_string,*) "ca_jsr.dat: calcium in junctional SR"
    CALL write2log(log_string)
    
    OPEN(unit=13, file=TRIM(where2save) // "/RosDS.dat")
    !WRITE(log_string,*) "nRo.dat: number of RyR open"
    !CALL write2log(log_string)

    OPEN(unit=14, file=TRIM(where2save) // "/RosI.dat")!nLo.dat")      
    WRITE(log_string,*) "nLo.dat: number of DHPR open"
    CALL write2log(log_string)
    
    OPEN(unit=15, file=TRIM(where2save) // "/Vm.dat")    
    OPEN(unit=16, file=TRIM(where2save) // "/Cmyo.dat")
    OPEN(unit=17, file=TRIM(where2save) // "/Cnsr.dat")
    OPEN(unit=18, file=TRIM(where2save) // "/state_Indx.dat")
    OPEN(unit=19, file=TRIM(where2save) // "/CaB_htrpn.dat")
    OPEN(unit=20, file=TRIM(where2save) // "/CaB_ltrpn.dat")
    OPEN(unit=21, file=TRIM(where2save) // "/CaB_trpn.dat")
    !OPEN(unit=18, file=TRIM(where2save) // "/Idhprtotal.dat")
  END SUBROUTINE datamode1_init
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE datamode1_print()
    INTEGER :: jj
    CHARACTER(100) :: formatstring

    WRITE(formatstring, "(i0)") NSFU+1
    formatstring = "(" // TRIM(formatstring) //  "(2x,1pe12.5))"
    
    ! write out the last 2 second only
    !IF (time .GE. t_end - 2.0d0) THEN          
       !! print ca_ds
       CaB_trpn = CaB_htrpn + CaB_ltrpn
       Ca_ds = Ca_ds_dev
       Ca_jsr = Ca_jsr_dev
       RYRgate = Ryrgate_dev
       DHPRgate = DHPRgate_dev
       CRU_states = CRU_states_dev
       CALL state_breakdown_all(CRU_states, irow_R,LCCcluster_states, RYRcluster_states)

       IF ((time - time_prev_lessfreq) .GE. 5.0d-3) THEN 
         !WRITE(11,TRIM(formatstring)) time,(Ca_ds(jj), jj=1,NSFU)
         !WRITE(12,TRIM(formatstring)) time,(Ca_jsr(jj), jj=1,NSFU)
         !WRITE(18,TRIM(formatstring)) time, (real(RYRcluster_states(jj)), jj =1,NSFU)            
         time_prev_lessfreq = time
       END IF


       WRITE(13,TRIM(formatstring)) time, RosDS!(RYRgate(jj), jj=1,NSFU)
       WRITE(14,TRIM(formatstring)) time, RosI !(DHPRgate(jj), jj=1,NSFU)
       WRITE(15,TRIM(formatstring)) time, Vm
       WRITE(16,TRIM(formatstring)) time, Ca_myo
       WRITE(17,TRIM(formatstring)) time, Ca_nsr
       WRITE(19,TRIM(formatstring)) time, CaB_htrpn !(RYRgate(jj), jj=1,NSFU)
       WRITE(20,TRIM(formatstring)) time, CaB_ltrpn !(RYRgate(jj), jj=1,NSFU)
       WRITE(21,TRIM(formatstring)) time, CaB_trpn !(RYRgate(jj), jj=1,NSFU)
       
       !PRINT *,"RosDS: ",RosDS,"RosProd: ",RosProduced
       !WRITE(18, *) I_dhpr_total(1)             
    !END IF
  END SUBROUTINE datamode1_print
  SUBROUTINE datamode1_close()
    ! close output files
    CLOSE(11)
    CLOSE(12)
    CLOSE(13)
    CLOSE(14)
    CLOSE(15)
    CLOSE(16)
    CLOSE(17)
    CLOSE(18)
    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
  END SUBROUTINE datamode1_close
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode2_init()
    IMPLICIT NONE
    LOGICAL :: is_exist
    INTEGER :: hdf_error
    
    WRITE(log_string,*) "Mode 2 (HDF)"
    CALL write2log(log_string)

    hdf_error =  hdf_init()
    CHUNK_SIZE = TRACK_LENGTH
    CHUNKS_PER_FILE = 1 ! set it to one for now
    NUM_STATES_NJ = mR * mL
    NUM_FIELDS4CURRENTS = 11
    FIELDS_PER_RECORD = 4 + 3*NSFU + NUM_STATES_NJ + NUM_FIELDS4CURRENTS


    ! + 2 = CaF, J_ryr_orphan
    FIELDS_PER_RECORD = FIELDS_PER_RECORD + 2
    
    num_record = 0
    file_index = 0
    ALLOCATE(data_buffer(FIELDS_PER_RECORD, CHUNK_SIZE))
    !remove existing files
    sys_command = "rm " // TRIM(where2save) // "/*.h5"
    CALL system(sys_command)
  END SUBROUTINE datamode2_init
  SUBROUTINE datamode2_close()
    INTEGER :: hdf_error
    
    hdf_error = hdf_close()
  END SUBROUTINE datamode2_close
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode2_print()
    INTEGER :: offset
    INTEGER :: ii
    ! track and output data
    IF (num_record .EQ. CHUNK_SIZE) THEN
       file_index = file_index + 1
       WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save), file_index         
       CALL hdf_open_simple(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
       CALL hdf_write_simple(file_struct, data_buffer, CHUNK_SIZE, FIELDS_PER_RECORD)
       CALL hdf_spark_close(file_struct)
       num_record = 0
    ENDIF
    num_record = num_record + 1
    data_buffer(1, num_record) = time
    data_buffer(2, num_record) = Vm
    data_buffer(3, num_record) = Ca_myo
    data_buffer(4, num_record) = Ca_nsr
    offset = 4
    data_buffer(offset+1:offset+NSFU, num_record) = CRU_states_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) =  Ca_ds_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) = Ca_jsr_dev
    offset = offset + NSFU
    !DO ii = 1, mL
    !   p_nj((ii-1)*mR+1:ii*mR) = pL_nj(ii) * pR_nj(1:mR)
    !END DO
    !data_buffer(offset+1:offset+NUM_STATES_NJ, num_record) = p_nj
    !offset = offset + NUM_STATES_NJ
    data_buffer(offset+1:offset+NUM_FIELDS4CURRENTS, num_record) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)
    
    offset = offset + NUM_FIELDS4CURRENTS
    data_buffer(offset+1,num_record) = CaF
    data_buffer(offset+2,num_record) = J_ryr_orphan
  END SUBROUTINE datamode2_print
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode2_print_last()
    file_index = file_index + 1
    WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save), file_index
    CALL hdf_open_simple(file_out, NSFU, num_record, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
    CALL hdf_write_simple(file_struct, data_buffer, num_record, FIELDS_PER_RECORD)
    CALL hdf_spark_close(file_struct)
  END SUBROUTINE datamode2_print_last
  
  SUBROUTINE datamode2_print_fb()
    !print out the first-beat
    INTEGER :: offset
    INTEGER :: ii
    ! track and output data
    IF (num_record .EQ. CHUNK_SIZE) THEN
       file_index = file_index + 1
       WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save)//'/'//TRIM(first_beat), file_index         
       CALL hdf_open_simple(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
       CALL hdf_write_simple(file_struct, data_buffer, CHUNK_SIZE, FIELDS_PER_RECORD)
       CALL hdf_spark_close(file_struct)
       num_record = 0
    ENDIF
    num_record = num_record + 1
    data_buffer(1, num_record) = time
    data_buffer(2, num_record) = Vm
    data_buffer(3, num_record) = Ca_myo
    data_buffer(4, num_record) = Ca_nsr
    offset = 4
    data_buffer(offset+1:offset+NSFU, num_record) = CRU_states_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) =  Ca_ds_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) = Ca_jsr_dev
    offset = offset + NSFU
    !DO ii = 1, mL
    !   p_nj((ii-1)*mR+1:ii*mR) = pL_nj(ii) * pR_nj(1:mR)
    !END DO
    !data_buffer(offset+1:offset+NUM_STATES_NJ, num_record) = p_nj
    !offset = offset + NUM_STATES_NJ
    data_buffer(offset+1:offset+NUM_FIELDS4CURRENTS, num_record) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)

    offset = offset + NUM_FIELDS4CURRENTS
    data_buffer(offset+1,num_record) = CaF
    data_buffer(offset+2,num_record) = J_ryr_orphan
  END SUBROUTINE datamode2_print_fb
  SUBROUTINE datamode2_print_fb_last()
    file_index = file_index + 1
    WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save)//'/'//TRIM(first_beat), file_index
    CALL hdf_open_simple(file_out, NSFU, num_record, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
    CALL hdf_write_simple(file_struct, data_buffer, num_record, FIELDS_PER_RECORD)
    CALL hdf_spark_close(file_struct)
    file_index =0
    num_record = 0
  END SUBROUTINE datamode2_print_fb_last

    
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!
#if _DATAMODE_ == _IO_2_ASCII
  SUBROUTINE datamode3_init()
    IMPLICIT NONE
    LOGICAL :: is_exist
    INTEGER :: offset
      
    WRITE(log_string,*) "Mode 3"
    CALL write2log(log_string)
      
    INQUIRE(FILE = TRIM(where2save), EXIST = is_exist)
    IF (.NOT. is_exist) THEN 
       sys_command = 'mkdir ' // TRIM(where2save) 
       CALL system(sys_command)
    END IF
    WRITE(log_string,*) "location: ", TRIM(where2save)
    CALL write2log(log_string)

    num_record = 0
    file_index = 0
    offset = 5 ! should be the value the same as the one being used in datamode3_print
    NUM_STATES_NJ = mR*mL
    NUM_FIELDS4CURRENTS = 11
    FIELDS_PER_RECORD = offset + 3*NSFU + NUM_STATES_NJ + NUM_FIELDS4CURRENTS
    ALLOCATE(data_buffer(FIELDS_PER_RECORD, TRACK_LENGTH) )

    WRITE(fstring, "(i0a)") offset, PRECISION_OUTPUT
    WRITE(fstring, "(a,',',i0a)") TRIM(fstring), 3*NSFU, PRECISION_OUTPUT
    WRITE(fstring, "(a,',',i0a)") TRIM(fstring), NUM_STATES_NJ, PRECISION_OUTPUT
    WRITE(fstring, "(a,',',i0a)") TRIM(fstring), NUM_FIELDS4CURRENTS, PRECISION_OUTPUT
      
    fstring = "(" // TRIM(fstring) // ")"
  END SUBROUTINE datamode3_init
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode3_close()
    CLOSE(unit=11)
  END SUBROUTINE datamode3_close
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode3_print()
    INTEGER :: offset
    INTEGER :: fileID
    
    fileID = get_file_ID()
    IF (num_record .EQ. TRACK_LENGTH) THEN
       file_index = file_index + 1
       WRITE(file_out, "(a,'/data-', i3.3,'.txt')") TRIM(where2save), file_index
       OPEN(unit=fileID, file=TRIM(file_out))
       DO ii = 1, TRACK_LENGTH
          WRITE(fileID,TRIM(fstring)) data_buffer(:, ii)
       END DO
       CLOSE(unit = fileID)
       num_record = 0
    END IF
    num_record = num_record + 1
    data_buffer(1, num_record) = time
    data_buffer(2, num_record) = Vm
    data_buffer(3, num_record) = I_dhpr_total(1)
    data_buffer(4, num_record) = Ca_myo
    data_buffer(5, num_record) = Ca_nsr
    offset = 5
    data_buffer(offset+1:offset+NSFU, num_record) = CRU_states_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) =  Ca_ds_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) = Ca_jsr_dev
    offset = offset + NSFU
    !DO ii = 1, mL
    !   p_nj((ii-1)*mR+1:ii*mR) = pL_nj(ii) * pR_nj(1:mR)
    !END DO
    !data_buffer(offset+1:offset+NUM_STATES_NJ, num_record) = p_nj
    !offset = offset + NUM_STATES_NJ
    data_buffer(offset+1:offset+NUM_FIELDS4CURRENTS, num_record) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)
    offset = offset + NUM_FIELDS4CURRENTS
  END SUBROUTINE datamode3_print
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode3_print_last()
    INTEGER :: fileID
    
    fileID = get_file_ID()
    
    file_index = file_index + 1
    WRITE(file_out, "(a,'/data-', i3.3,'.txt')") TRIM(where2save), file_index
    OPEN(unit=fileID, file=TRIM(file_out))
    
    DO ii = 1, num_record
       WRITE(fileID,TRIM(fstring)) data_buffer(:, ii)
    END DO

    CLOSE(unit = fileID)
    
  END SUBROUTINE datamode3_print_last

#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE data_io_write()
    !WRITE INITIAL DATA OUT
#if _DATAMODE_ == _IO_TEST
    CALL datamode1_print()
#elif _DATAMODE_ == _IO_2_HDF
    CALL datamode2_print()
#elif _DATAMODE_ == _IO_2_ASCII
    CALL datamode3_print()
#elif _DATAMODE_ == _IO_2_HDF2
    CALL datamode7_print()
#endif
  END SUBROUTINE data_io_write

  SUBROUTINE data_io_finalize()
    INCLUDE "lib3f.h"
    INTEGER :: ierr, ii, fileID
      
    ! writing the last time
#if _DATAMODE_ == _IO_TEST
    CALL datamode1_close()
#elif _DATAMODE_ == _IO_2_HDF
    CALL datamode2_print_last()
    CALL datamode2_close()
#elif  _DATAMODE_ == _IO_2_ASCII
    CALL datamode3_print_last()  ! print the last time
    CALL datamode3_close()
      !   CALL datamode3_finalize()
!!$#elif _DATAMODE_ == 3
!!$    CLOSE(11)
!!$    sys_command = 'gzip -c ' // trim(filename) // ' > ' // trim(filename) // '.gz'
!!$    CALL system(sys_command)
!!$    sys_command = 'rm ' // trim(filename)
!!$    CALL system(sys_command)
#elif _DATAMODE_ == _IO_2_HDF2
    CALL datamode7_print_last()
    CALL datamode7_close()      
#endif
    IF (write_restartfile .NE. 0) THEN
       CALL write2restartfile()
    ENDIF
            
    !!send email message
    ierr = hostnm(hostname)
    sys_command = 'echo "' // TRIM(full_command) // ' \n run on machine ' //  TRIM(hostname) // ' complete" | mail -s "PROGRAM COMPLETED" hoangtrongminhtuan@gmail.com'
    ierr = system(TRIM(sys_command))

  END SUBROUTINE data_io_finalize

  SUBROUTINE loaddata_from_conf()
    IMPLICIT NONE
      
    ! Concentration parameters + initial condition
    Ca_ds = Ca_myo
    Ca_jsr = Ca_nsr
    CaF = CaF_infty
    
    ! init data use for EC-coupling model
    CALL init_data_ECC()

    CaB_htrpn = (Ca_myo * B_htrpnT) / ( koff_htrpn/kon_htrpn + Ca_myo)
    CaB_ltrpn = (Ca_myo * B_ltrpnT) / ( koff_ltrpn/kon_ltrpn + Ca_myo)    
    
#if _TEST_ORPHAN_CRU_ == _YES
    ALLOCATE(pR_orphan(irow_R_orphan), tmp_pR_orphan(irow_R_orphan))
    pR_orphan = 0d0
    ALLOCATE(chiC_orphan(irow_R_orphan), chiO_orphan(irow_R_orphan))
#endif
  END SUBROUTINE loaddata_from_conf

  !! we need
  !!  - gating variables
  !!  - Camyo, Cansr, Cads, Cajsr, Na_i, K_i...
  !!  - CRU_states
  SUBROUTINE loaddata_from_restartfile()
    IMPLICIT NONE
    LOGICAL :: is_exist
    INTEGER :: offset, ii, fileID
    CHARACTER(100) :: formatstring
    
    WRITE(log_string, *) "Running simulation using input from restart file"
    CALL write2log(log_string)
    ! check if restart file existed
    INQUIRE(FILE = TRIM(restart_file), EXIST = is_exist)
    IF (.NOT. is_exist) THEN
       WRITE(log_string, *) "Restart " // TRIM(restart_file) // " file not found"
       CALL write2log(log_string)
       STOP
    END IF

    ! check if restart file data fit wit conf. data file
    NUM_FIELDS4CURRENTS = 11
    NUM_INDIVIDUAL_FIELDS = 7+NUM_APD_TYPES
    
#if _TEST_ORPHAN_CRU_ == _YES
    RF_FIELDS_PER_RECORD = 3*FULL_NSFU + mL + mR + NUM_FIELDS4CURRENTS + NUM_INDIVIDUAL_FIELDS
#else
    RF_FIELDS_PER_RECORD = 3*NSFU + mL + mR + NUM_FIELDS4CURRENTS + NUM_INDIVIDUAL_FIELDS
#endif

    ALLOCATE(data_buffer(RF_FIELDS_PER_RECORD,1))
    
    WRITE(formatstring, "(i0)") RF_FIELDS_PER_RECORD  
    formatstring = "(" // TRIM(formatstring) //  "(2x,1pe21.13))"

    fileID = get_file_ID()
    OPEN (unit = fileID, file = TRIM(restart_file))
    READ (fileID, TRIM(formatstring)) data_buffer
    CLOSE (fileID)

#if _TEST_ORPHAN_CRU_ == _YES
    ! a fraction of CRU becomes orphan, i.e. LCC becomes non-junctional
    ! and RyR clusters see Ca_myo, not Ca_ds
    IF (TRIM(restart_file) .EQ. "") THEN
       WRITE (6, *) "Orphan mode only work with restart file, for now"
       WRITE (6, *) "and requires 20,000 CaRUs in control simul"
       STOP
    END IF

    offset = 0
    CRU_states = data_buffer(offset+1:offset+NSFU,1)
    DO ii = 1, irow_R_orphan
       pR_orphan(ii) = COUNT(data_buffer(4+NSFU+1:4+FULL_NSFU,1) .EQ. ii) 
    END DO

    offset = offset + FULL_NSFU
    Ca_ds = data_buffer(offset+1:offset+NSFU,1)
    offset = offset+ FULL_NSFU
    Ca_jsr = data_buffer(offset+1:offset+NSFU,1)
    offset = offset + FULL_NSFU
    !pL_nj = data_buffer(offset+1:offset+mL,1)
    !offset = offset + mL
    !pR_nj = data_buffer(offset+1:offset+mR,1)
    !offset = offset + mR
               
#else
    ! set up data
    offset = 0
    CRU_states = data_buffer(offset+1:offset+NSFU,1)
    offset = offset + NSFU
    Ca_ds = data_buffer(offset+1:offset+NSFU,1)
    offset = offset+NSFU
    Ca_jsr = data_buffer(offset+1:offset+NSFU,1)
    offset = offset + NSFU
    !pL_nj = data_buffer(offset+1:offset+mL,1)
    !offset = offset + mL
    !pR_nj = data_buffer(offset+1:offset+mR,1)
    !offset = offset + mR
#endif
      
    Na_i = data_buffer(offset + 1, 1)
    K_i = data_buffer(offset + 2,1)
    m_na = data_buffer(offset + 3,1)
    h_na = data_buffer(offset + 4,1)
    j_na = data_buffer(offset + 5,1)
    a_ktof = data_buffer(offset + 6,1)
    in_ktof = data_buffer(offset + 7,1)
    a_ktos = data_buffer(offset + 8,1)
    in_ktos = data_buffer(offset + 9,1)
    a_kss = data_buffer(offset + 10,1)
    in_kss = data_buffer(offset + 11,1)
    offset = offset + NUM_FIELDS4CURRENTS

    Vm = data_buffer(offset+1,1)
    Ca_myo = data_buffer(offset+2,1)
    Ca_nsr = data_buffer(offset+3,1)
    is_being_clamped = data_buffer(offset + 4,1)
    CaB_htrpn = data_buffer(offset + 5, 1)
    CaB_ltrpn = data_buffer(offset + 6, 1)
    CaF = data_buffer(offset + 7, 1)
    is_inside_AP = data_buffer(offset + 8:offset+7+NUM_APD_TYPES, 1)
    offset = offset + 7 + NUM_APD_TYPES
    
    real_timestart = data_buffer(offset + 1, 1)
    real_timestart = 0.0d0 !I want it to start from 0 and not from where the restart_file ended
    offset =  offset + 1
  END SUBROUTINE loaddata_from_restartfile

  !!>****f* 
  !! SYNOPSIS  
  !!   Setup kernel configuration 
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  FUNCTION setup_KernelConfiguration() RESULT(ierr)
    IMPLICIT NONE
    INTEGER :: ierr
    !local
    INTEGER :: err
    
    ierr = SUCCESS
    ! create cuda stream
    err = cudaStreamCreate(stream0)
!!$    IF (err .NE. cudaSuccess) THEN
!!$       WRITE (error_message,*) "ERROR: Create cuda stream"
!!$       ierr = ERR_CUDA_CREATE_STREAM
!!$       RETURN
!!$    ENDIF
    err = cudaStreamCreate(stream1)
!!$    IF (err .NE. cudaSuccess) THEN
!!$       WRITE (error_message,*) "ERROR: Create cuda stream"
!!$       ierr = ERR_CUDA_CREATE_STREAM
!!$       RETURN
!!$    ENDIF

    !for 1D thread BLOCK
    blocksize = BLOCKSIZE_1D
    !gridSize = INT((NSFU+blocksize-1)/blocksize)  ! ceiling value
    gridSize = CEILING(REAL(NSFU)/blocksize)  ! ceiling value
    RETURN
  END FUNCTION  setup_KernelConfiguration

  
  SUBROUTINE init_PRNG()
    IMPLICIT NONE
    INTEGER :: global_seed
    REAL(KIND=dp) :: tempsec, tempsec2

    !choose PRNG
#if _USE_PRNG_ == _PRNG_SARU
    WRITE(6,*) "Use SARU"
    ! initialize the pseudo-RNG
    IF (fixed_seed .EQ. 1) THEN
       global_seed = rand_seed
       WRITE(log_string,*) "use fixed-seed"
       CALL write2log(log_string)
    ELSE
       tempsec = 372.d0*SECNDS(0.0)-1.d0
       tempsec2 = 2.d0*SECNDS(0.0)-1.d0
       global_seed = INT(tempsec) + MOD(INT(tempsec2),127)*975121
    END IF
    saruObj = init_saru(NSFU, global_seed)
#elif _USE_PRNG_ == _PRNG_MT
    WRITE(log_string,*) "Use MT GPU"
    CALL write2log(log_string)
    d_status = init_mt_status()
    d_texture = init_mt_texture(d_status)
!!$    num_data = 1000
!!$    res = cudaMalloc(d_data, num_data)
#elif _USE_PRNG_ == _PRNG_MTCPU
    WRITE(log_string,*) "Use MT CPU"
    CALL write2log(log_string)
    IF (fixed_seed .EQ. 1) THEN
       CALL init_rndgen_testonly(mt_STATE, INFO, rand_seed)
    ELSE
       CALL init_rndgen(mt_STATE, INFO)
    END IF
    ALLOCATE(X_r(NSFU*copy_interval))
    mt_LB = 0.0
    mt_UB = 1.0
#endif

  END SUBROUTINE init_PRNG


  SUBROUTINE spark_summarize()
    IMPLICIT NONE
    INTEGER :: ii, fileID
    REAL(KIND=dp), DIMENSION(NSFU) :: tmp
    
#if _RUNMODE_ == _MODE_V_CLAMP || _RUNMODE_ == _MODE_I_CLAMP
    IF (accumulate_J_dhpr .LE. EPS) THEN
       GAIN = 0.0d0
    ELSE
       GAIN = accumulate_J_ryr / accumulate_J_dhpr
    ENDIF
    DO ii = 1, NUM_APD_TYPES
       IF (APcount(ii) .EQ. 0) THEN
          APD(ii) = 0
       ELSE
          APD(ii) = APD(ii)/ APcount(ii)
       ENDIF
    END DO
#else
    GAIN = 0.0d0
    APD = 0
#endif
    summary_data(1) = GAIN
    description_summarydata(1) = 'GAIN'
    count_summarydata = 1
    
    summary_data(count_summarydata+1:count_summarydata+NUM_APD_TYPES) = APD
    description_summarydata(count_summarydata+1:count_summarydata+NUM_APD_TYPES) = (/ "APD90   ", "APD50", 'APD20', "APD0mV" /)
    count_summarydata = NUM_APD_TYPES + count_summarydata
    
    summary_data(count_summarydata+1) = t_start 
    summary_data(count_summarydata+2) = t_end
    description_summarydata(count_summarydata+1) = "time_start"
    description_summarydata(count_summarydata+2) = "time_end"
    count_summarydata = count_summarydata + 2
    
    summary_data(count_summarydata+1) = (accumulate_J_ryr-accumulate_J_ryr_orphan)/(t_end-time_start_spark_analysis) * Vol_ref/V_myo_T
    summary_data(count_summarydata+2) = accumulate_J_ryr_orphan/(t_end-time_start_spark_analysis) * Vol_ref/V_myo_T  
    description_summarydata(count_summarydata+1) = "mass_Ca_release_normal_per_sec"
    description_summarydata(count_summarydata+2) = "mass_Ca_release_orphan_per_sec"
    count_summarydata = count_summarydata + 2

    Nsparks = Nsparks_dev
    Nblinks = Nblinks_dev
    Nquarks = Nquarks_dev
#if _TEST_NEW_ORPHAN_CRU_ == _YES    
    num_Sparks = SUM(Nsparks(1:NSFU_normal))
    num_Blinks = SUM(Nblinks(1:NSFU_normal))
    num_Quarks = SUM(Nquarks(1:NSFU_normal))
    IF (num_Sparks .EQ. 0)  num_Sparks = 1
    IF (num_Blinks .EQ. 0) num_Blinks = 1
    IF (num_Quarks .EQ. 0) num_Quarks = 1

    tmp = tot_SparkDurations
    avg_SparkDuration = SUM(tmp(1:NSFU_normal))/num_Sparks
    tmp = tot_BlinkDurations
    avg_BlinkDuration = SUM(tmp(1:NSFU_normal))/num_Blinks
    tmp = tot_SparkPeaks
    avg_SparkPeak = SUM(tmp(1:NSFU_normal))/num_Sparks

#else    
    num_Sparks = SUM(Nsparks)
    num_Blinks = SUM(Nblinks)
    num_Quarks = SUM(Nquarks)
    IF (num_Sparks .EQ. 0)  num_Sparks = 1
    IF (num_Blinks .EQ. 0) num_Blinks = 1
    IF (num_Quarks .EQ. 0) num_Quarks = 1

    tmp = tot_SparkDurations
    avg_SparkDuration = SUM(tmp)/num_Sparks
    tmp = tot_BlinkDurations
    avg_BlinkDuration = SUM(tmp)/num_Blinks
    tmp = tot_SparkPeaks
    avg_SparkPeak = SUM(tmp)/num_Sparks
#endif

    IF (num_Quarks+num_Sparks .GT. 0) spark_Freq = float(num_Sparks)/(num_Quarks+num_Sparks)

    summary_data(count_summarydata+1) = num_Sparks/(t_end - time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+2) = num_Blinks/(t_end - time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+3) = num_Quarks/(t_end - time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+4) = avg_SparkDuration
    summary_data(count_summarydata+5) = avg_BlinkDuration
    summary_data(count_summarydata+6) = avg_SparkPeak
    description_summarydata(count_summarydata+1) = "num_DiastolicSparks_per_sec"
    description_summarydata(count_summarydata+2) = "num_DiastolicBlinks_per_sec"
    description_summarydata(count_summarydata+3) = "num_DiastolicQuarks_per_sec"
    description_summarydata(count_summarydata+4) = "avg_SparksDuration"
    description_summarydata(count_summarydata+5) = "avg_BlinkDuration"
    description_summarydata(count_summarydata+6) = "avg_SparkPeak"
    count_summarydata = count_summarydata + 6

    summary_data(count_summarydata+1) = spark_Freq
!!$    summary_data(count_summarydata+2) = SFWHM
!!$    summary_data(count_summarydata+3) = BFWHM
    description_summarydata(count_summarydata+1) = "Spark_Freq"
!!$    description_summarydata(count_summarydata+2) = "Spark_FWHM"
!!$    description_summarydata(count_summarydata+3) = "Blink_FWHM"
!!$    count_summarydata = count_summarydata + 3
    count_summarydata = count_summarydata + 1

    summary_data(count_summarydata+1) = time_start_DiastolicPhase
    description_summarydata(count_summarydata+1) = "time_start_DiastolicPhase"
    count_summarydata = count_summarydata + 1
    
#if _TEST_NEW_ORPHAN_CRU_ == _YES
    num_Sparks = SUM(Nsparks(NSFU_normal+1:))
    num_Blinks = SUM(Nblinks(NSFU_normal+1:))
    num_Quarks = SUM(Nquarks(NSFU_normal:))
    avg_SparkDuration = SUM(tot_SparkDurations(NSFU_normal+1:))/num_Sparks
    avg_BlinkDuration = SUM(tot_BlinkDurations(NSFU_normal+1:))/num_Blinks
    avg_SparkPeak = SUM(tot_SparkPeaks(NSFU_normal+1:))/num_Sparks

    
    IF (num_Quarks+num_Sparks .GT. 0) spark_Freq = float(num_Sparks)/(num_Quarks+num_Sparks)
    
    summary_data(count_summarydata+1) = num_Sparks/(t_end-time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+2) = num_Blinks/(t_end-time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+3) = num_Quarks/(t_end-time_start_spark_analysis) * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+4) = avg_SparkDuration
    summary_data(count_summarydata+5) = avg_BlinkDuration
    summary_data(count_summarydata+6) = avg_SparkPeak
    description_summarydata(count_summarydata+1) = "num_Sparks_per_sec_Orphan"
    description_summarydata(count_summarydata+2) = "num_Blinks_per_sec_Orphan"
    description_summarydata(count_summarydata+3) = "num_Quarks_per_sec_Orphan"
    description_summarydata(count_summarydata+4) = "avg_SparksDuration_Orphan"
    description_summarydata(count_summarydata+5) = "avg_BlinkDuration_Orphan"
    description_summarydata(count_summarydata+6) = "avg_SparkPeak_Orphan"
    count_summarydata = count_summarydata + 6
#endif
    !write Camyo, Cansr at the end
    summary_data(count_summarydata+1) = Ca_myo
    summary_data(count_summarydata+2) = Ca_nsr
    description_summarydata(count_summarydata+1) = "Camyo"
    description_summarydata(count_summarydata+2) = "Cansr"
    count_summarydata = count_summarydata + 2
    
    fileID = get_file_ID()
    !! write summary data to file
    fstring = TRIM(where2save) // "/" // summary_filename
    OPEN(unit=fileID, file=TRIM(fstring))
    !! write config file name to the first line
    WRITE(fileID, *) TRIM(conf_file)
    DO ii = 1, count_summarydata
       WRITE(fileID,"(1pe12.5,'   ', A)") summary_data(ii), TRIM(description_summarydata(ii))
    END DO
    CLOSE(fileID)
  END SUBROUTINE spark_summarize

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode7_init()
    IMPLICIT NONE
    LOGICAL :: is_exist
    INTEGER :: hdf_error
    
    WRITE(log_string,*) "Mode 2 (HDF)"
    CALL write2log(log_string)
    
    hdf_error =  hdf_init()
    CHUNK_SIZE = TRACK_LENGTH
    CHUNKS_PER_FILE = 1 ! set it to one for now
    NUM_STATES_NJ = mR * mL
    NUM_FIELDS4CURRENTS = 11
    FIELDS_PER_RECORD = 4 + 3*NSFU + NUM_STATES_NJ + NUM_FIELDS4CURRENTS


    ! + 2 = CaF, J_ryr_orphan
    FIELDS_PER_RECORD = FIELDS_PER_RECORD + 2
    
    num_record = 0
    file_index = 0
    ALLOCATE(data_buffer(FIELDS_PER_RECORD, CHUNK_SIZE))
    !remove existing files
    sys_command = "rm " // TRIM(where2save) // "/*.h5"
    CALL system(sys_command)
  END SUBROUTINE datamode7_init
  SUBROUTINE datamode7_close()
    INTEGER :: hdf_error
    
    hdf_error = hdf_close()
  END SUBROUTINE datamode7_close
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode7_print()
    INTEGER :: offset
    INTEGER :: ii
    ! track and output data
    IF (num_record .EQ. CHUNK_SIZE) THEN
       file_index = file_index + 1
       WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save), file_index         
       CALL hdf_ecc_open(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
       CALL hdf_ecc_write(file_struct, data_buffer, CHUNK_SIZE, FIELDS_PER_RECORD, NSFU)
       CALL hdf_ecc_close(file_struct)
       num_record = 0
    ENDIF
    num_record = num_record + 1
    data_buffer(1, num_record) = time
    data_buffer(2, num_record) = Vm
    data_buffer(3, num_record) = Ca_myo
    data_buffer(4, num_record) = Ca_nsr
    offset = 4
    !DO ii = 1, mL
    !   p_nj((ii-1)*mR+1:ii*mR) = pL_nj(ii) * pR_nj(1:mR)
    !END DO
    !data_buffer(offset+1:offset+NUM_STATES_NJ, num_record) = p_nj
    !offset = offset + NUM_STATES_NJ
    data_buffer(offset+1:offset+NUM_FIELDS4CURRENTS, num_record) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)
    offset = offset + NUM_FIELDS4CURRENTS
    data_buffer(offset+1,num_record) = CaF
    data_buffer(offset+2,num_record) = J_ryr_orphan
    offset = offset + 2
    
    data_buffer(offset+1:offset+NSFU, num_record) = CRU_states_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) =  Ca_ds_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) = Ca_jsr_dev
    offset = offset + NSFU
  END SUBROUTINE datamode7_print
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE datamode7_print_last()
    file_index = file_index + 1
    WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save), file_index
    CALL hdf_ecc_open(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
    CALL hdf_ecc_write(file_struct, data_buffer, num_record, FIELDS_PER_RECORD, NSFU)
    CALL hdf_ecc_close(file_struct)
  END SUBROUTINE datamode7_print_last

  SUBROUTINE datamode7_print_fb()
    !print out the first-beat
    INTEGER :: offset
    INTEGER :: ii
    ! track and output data
    IF (num_record .EQ. CHUNK_SIZE) THEN
       file_index = file_index + 1
       WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save)//'/'//TRIM(first_beat), file_index         
       CALL hdf_ecc_open(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
       CALL hdf_ecc_write(file_struct, data_buffer, CHUNK_SIZE, FIELDS_PER_RECORD, NSFU)
       CALL hdf_ecc_close(file_struct)
       num_record = 0
    ENDIF
    num_record = num_record + 1
    data_buffer(1, num_record) = time
    data_buffer(2, num_record) = Vm
    data_buffer(3, num_record) = Ca_myo
    data_buffer(4, num_record) = Ca_nsr
    offset = 4
    !DO ii = 1, mL
    !   p_nj((ii-1)*mR+1:ii*mR) = pL_nj(ii) * pR_nj(1:mR)
    !END DO
    !data_buffer(offset+1:offset+NUM_STATES_NJ, num_record) = p_nj
    !offset = offset + NUM_STATES_NJ
    data_buffer(offset+1:offset+NUM_FIELDS4CURRENTS, num_record) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)
    
    offset = offset + NUM_FIELDS4CURRENTS
    data_buffer(offset+1,num_record) = CaF
    data_buffer(offset+2,num_record) = J_ryr_orphan
    offset = offset + 2
    data_buffer(offset+1:offset+NSFU, num_record) = CRU_states_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) =  Ca_ds_dev
    offset = offset + NSFU
    data_buffer(offset+1:offset+NSFU, num_record) = Ca_jsr_dev
    offset = offset + NSFU
    
  END SUBROUTINE datamode7_print_fb
  SUBROUTINE datamode7_print_fb_last()
    file_index = file_index + 1
    WRITE(file_out, "(a,'/data-', i3.3,'.h5')") TRIM(where2save)//'/'//TRIM(first_beat), file_index
    CALL hdf_ecc_open(file_out, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, FIELDS_PER_RECORD)
    CALL hdf_ecc_write(file_struct, data_buffer, CHUNK_SIZE, FIELDS_PER_RECORD, NSFU)
    CALL hdf_ecc_close(file_struct)
    file_index =0
    num_record = 0
  END SUBROUTINE datamode7_print_fb_last


  SUBROUTINE data_io_write_fb()
    !WRITE INITIAL DATA OUT
#if _DATAMODE_ == _IO_TEST
    CALL datamode1_print_fb()
#elif _DATAMODE_ == _IO_2_HDF
    CALL datamode2_print_fb()
#elif _DATAMODE_ == _IO_2_ASCII
    !CALL datamode3_print()
#elif _DATAMODE_ == _IO_2_HDF2
    CALL datamode7_print_fb()
#endif
  END SUBROUTINE data_io_write_fb
  
  SUBROUTINE data_io_write_fb_last()
    integer :: done_fb = 1
    !WRITE INITIAL DATA OUT
#if _DATAMODE_ == _IO_TEST
    CALL datamode1_print_fb(done_fb)
#elif _DATAMODE_ == _IO_2_HDF
    CALL datamode2_print_fb_last()
#elif _DATAMODE_ == _IO_2_ASCII
    !CALL datamode3_print()
#elif _DATAMODE_ == _IO_2_HDF2
    CALL datamode7_print_fb_last()
#endif
  END SUBROUTINE data_io_write_fb_last

  SUBROUTINE write2restartfile()
    IMPLICIT NONE
    INTEGER :: fileID
    CHARACTER(100) :: formatstring, fstring
    INTEGER :: offset
    
    WRITE(log_string,*) "Print out restart file"
    CALL write2log(log_string)
      
    NUM_FIELDS4CURRENTS = 11
    NUM_INDIVIDUAL_FIELDS = 7+NUM_APD_TYPES
    RF_FIELDS_PER_RECORD = 3*NSFU + mL + mR + NUM_FIELDS4CURRENTS + NUM_INDIVIDUAL_FIELDS

    ALLOCATE(data_2restartfile(RF_FIELDS_PER_RECORD,1))
    ! write out data to restart file
    offset = 0
    data_2restartfile(offset+1:offset+NSFU,1) = CRU_states_dev
    offset = offset + NSFU
    data_2restartfile(offset+1:offset+NSFU,1) =  Ca_ds_dev
    offset = offset + NSFU
    data_2restartfile(offset+1:offset+NSFU,1) = Ca_jsr_dev
    offset = offset + NSFU
    !data_2restartfile(offset+1:offset+mL, 1) = pL_nj
    !offset = offset + mL
    !data_2restartfile(offset+1:offset+mR, 1) = pR_nj
    !offset = offset + mR
    data_2restartfile(offset+1:offset+NUM_FIELDS4CURRENTS, 1) = (/ Na_i, K_i, m_na, h_na, j_na, a_ktof, in_ktof, a_ktos, in_ktos, a_kss, in_kss /)
    offset = offset + NUM_FIELDS4CURRENTS
    data_2restartfile(offset+1,1) = Vm
    data_2restartfile(offset+2,1) = Ca_myo
    data_2restartfile(offset+3,1) = Ca_nsr
    data_2restartfile(offset+4,1) = is_being_clamped
    data_2restartfile(offset+5,1) = CaB_htrpn
    data_2restartfile(offset+6,1) = CaB_ltrpn
    data_2restartfile(offset+7,1) = CaF
    data_2restartfile(offset+8:offset+7+NUM_APD_TYPES,1) = is_inside_AP
    offset = offset + 7 + NUM_APD_TYPES
    data_2restartfile(offset+1, 1) = t_end
    
    WRITE(formatstring, "(i0)") RF_FIELDS_PER_RECORD
    formatstring = "(" // TRIM(formatstring) //  "(2x,1pe21.13))"
    
    fstring = TRIM(where2save) // "/" // TRIM(restart_filename) 
    fileID = get_file_ID()
    OPEN (unit = fileID, file = TRIM(fstring))
    WRITE (fileID, TRIM(formatstring)) data_2restartfile
    CLOSE (fileID)
  END SUBROUTINE write2restartfile
    
  SUBROUTINE spark_summarize_fb()
    IMPLICIT NONE
    INTEGER :: ii
    INTEGER :: fileID
    REAL(KIND=dp), DIMENSION(NSFU) :: tmp

#if _RUNMODE_ == _MODE_V_CLAMP || _RUNMODE_ == _MODE_I_CLAMP
    GAIN = accumulate_J_ryr / accumulate_J_dhpr
    APD = APD/ APcount
#else
    GAIN = 0.0d0
    APD = 0
#endif
    summary_data(1) = GAIN
    description_summarydata(1) = 'GAIN'
    count_summarydata = 1
    
    summary_data(count_summarydata+1:count_summarydata+NUM_APD_TYPES) = APD
    description_summarydata(count_summarydata+1:count_summarydata+NUM_APD_TYPES) = (/ "APD90   ", "APD50", 'APD20', "APD0mV" /)
    count_summarydata = NUM_APD_TYPES + count_summarydata
    
    summary_data(count_summarydata+1) = t_start 
    summary_data(count_summarydata+2) = t_end
    description_summarydata(count_summarydata+1) = "time_start"
    description_summarydata(count_summarydata+2) = "time_end"
    count_summarydata = count_summarydata + 2
    
    summary_data(count_summarydata+1) = (accumulate_J_ryr-accumulate_J_ryr_orphan) * Vol_ref/V_myo_T 
    summary_data(count_summarydata+2) = accumulate_J_ryr_orphan * Vol_ref/V_myo_T 
    description_summarydata(count_summarydata+1) = "mass_Ca_release_normal_per_sec"
    description_summarydata(count_summarydata+2) = "mass_Ca_release_orphan_per_sec"
    count_summarydata = count_summarydata + 2

    Nsparks = Nsparks_dev
    Nblinks = Nblinks_dev
    Nquarks = Nquarks_dev
#if _TEST_NEW_ORPHAN_CRU_ == _YES    
    num_Sparks = SUM(Nsparks(1:NSFU_normal))
    num_Blinks = SUM(Nblinks(1:NSFU_normal))
    num_Quarks = SUM(Nquarks(1:NSFU_normal))
    IF (num_Sparks .EQ. 0)  num_Sparks = 1
    IF (num_Blinks .EQ. 0) num_Blinks = 1
    IF (num_Quarks .EQ. 0) num_Quarks = 1

    tmp = tot_SparkDurations
    avg_SparkDuration = SUM(tmp(1:NSFU_normal))/num_Sparks
    tmp = tot_BlinkDurations
    avg_BlinkDuration = SUM(tmp(1:NSFU_normal))/num_Blinks
    tmp = tot_SparkPeaks
    avg_SparkPeak = SUM(tmp(1:NSFU_normal))/num_Sparks
!!$    avg_SparkDuration = SUM(tot_SparkDurations(1:NSFU_normal))/num_Sparks 
!!$    avg_BlinkDuration = SUM(tot_BlinkDurations(1:NSFU_normal))/num_Blinks
!!$    avg_SparkPeak = SUM(tot_SparkPeaks(1:NSFU_normal))/num_Sparks
#else
    num_Sparks = SUM(Nsparks)
    num_Blinks = SUM(Nblinks)
    num_Quarks = SUM(Nquarks)
    IF (num_Sparks .EQ. 0)  num_Sparks = 1
    IF (num_Blinks .EQ. 0) num_Blinks = 1
    IF (num_Quarks .EQ. 0) num_Quarks = 1

!!$    avg_SparkDuration = SUM(tot_SparkDurations)/num_Sparks
!!$    avg_BlinkDuration = SUM(tot_BlinkDurations)/num_Blinks
!!$    avg_SparkPeak = SUM(tot_SparkPeaks)/num_Sparks
    tmp = tot_SparkDurations
    avg_SparkDuration = SUM(tmp)/num_Sparks
    tmp = tot_BlinkDurations
    avg_BlinkDuration = SUM(tmp)/num_Blinks
    tmp = tot_SparkPeaks
    avg_SparkPeak = SUM(tmp)/num_Sparks
#endif

    IF (num_Quarks+num_Sparks .GT. 0) spark_Freq = float(num_Sparks)/(num_Quarks+num_Sparks)
    
    summary_data(count_summarydata+1) = num_Sparks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+2) = num_Blinks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+3) = num_Quarks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+4) = avg_SparkDuration
    summary_data(count_summarydata+5) = avg_BlinkDuration
    summary_data(count_summarydata+6) = avg_SparkPeak
    description_summarydata(count_summarydata+1) = "num_DiastolicSparks_per_sec"
    description_summarydata(count_summarydata+2) = "num_DiastolicBlinks_per_sec"
    description_summarydata(count_summarydata+3) = "num_DiastolicQuarks_per_sec"
    description_summarydata(count_summarydata+4) = "avg_SparksDuration"
    description_summarydata(count_summarydata+5) = "avg_BlinkDuration"
    description_summarydata(count_summarydata+6) = "avg_SparkPeak"
    count_summarydata = count_summarydata + 6

    summary_data(count_summarydata+1) = spark_Freq
!!$    summary_data(count_summarydata+2) = SFWHM
!!$    summary_data(count_summarydata+3) = BFWHM
    description_summarydata(count_summarydata+1) = "Spark_Freq"
!!$    description_summarydata(count_summarydata+2) = "Spark_FWHM"
!!$    description_summarydata(count_summarydata+3) = "Blink_FWHM"
!!$    count_summarydata = count_summarydata + 3
    count_summarydata = count_summarydata + 1
    
    summary_data(count_summarydata+1) = time_start_DiastolicPhase
    description_summarydata(count_summarydata+1) = "time_start_DiastolicPhase"
    count_summarydata = count_summarydata + 1
    
#if _TEST_NEW_ORPHAN_CRU_ == _YES
    num_Sparks = SUM(Nsparks(NSFU_normal+1:))
    num_Blinks = SUM(Nblinks(NSFU_normal+1:))
    num_Quarks = SUM(Nquarks(NSFU_normal:))
    avg_SparkDuration = SUM(tot_SparkDurations(NSFU_normal+1:))/num_Sparks
    avg_BlinkDuration = SUM(tot_BlinkDurations(NSFU_normal+1:))/num_Blinks
    avg_SparkPeak = SUM(tot_SparkPeaks(NSFU_normal+1:))/num_Sparks
    
    IF (num_Quarks+num_Sparks .GT. 0) spark_Freq = float(num_Sparks)/(num_Quarks+num_Sparks)
    
    summary_data(count_summarydata+1) = num_Sparks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+2) = num_Blinks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+3) = num_Quarks * (float(MAX_CaRUs)/NSFU)
    summary_data(count_summarydata+4) = avg_SparkDuration
    summary_data(count_summarydata+5) = avg_BlinkDuration
    summary_data(count_summarydata+6) = avg_SparkPeak
    description_summarydata(count_summarydata+1) = "num_Sparks_per_sec_Orphan"
    description_summarydata(count_summarydata+2) = "num_Blinks_per_sec_Orphan"
    description_summarydata(count_summarydata+3) = "num_Quarks_per_sec_Orphan"
    description_summarydata(count_summarydata+4) = "avg_SparksDuration_Orphan"
    description_summarydata(count_summarydata+5) = "avg_BlinkDuration_Orphan"
    description_summarydata(count_summarydata+6) = "avg_SparkPeak_Orphan"
    count_summarydata = count_summarydata + 6
#endif
    !write Camyo, Cansr at the end
    summary_data(count_summarydata+1) = Ca_myo
    summary_data(count_summarydata+2) = Ca_nsr
    description_summarydata(count_summarydata+1) = "Camyo"
    description_summarydata(count_summarydata+2) = "Cansr"
    count_summarydata = count_summarydata + 2
      
    fileID = get_file_ID()
    !! write summary data to file
    fstring = TRIM(where2save) // "/" // TRIM(first_beat) // "/summary_data.txt"
    OPEN(unit=fileID, file=TRIM(fstring))
    !! write config file name to the first line
    WRITE(fileID, *) TRIM(conf_file), count_summarydata
    DO ii = 1, count_summarydata
       WRITE(11,"(1pe12.5,'   ', A)") summary_data(ii), TRIM(description_summarydata(ii))
    END DO
    CLOSE(fileID)

    GAIN = 0d0
    accumulate_J_dhpr = 0.0d0
    accumulate_J_ryr = 0.0d0; accumulate_J_ryr_orphan = 0.0d0
    APD = 0d0; APstart = 0d0; APend = 0d0;
    Vmthreshold = (/ V_rest * 0.9d0, V_rest * 0.5d0, V_rest * 0.2d0, +0.0d0 /)
    APcount = +0
    is_in_DiastolicPhase = .FALSE.
    is_done_firstbeat = .FALSE.
    
  END SUBROUTINE spark_summarize_fb

  attributes(device) SUBROUTINE state_breakdown_ker(cru_state, irow_R, lcc_state, ryr_state)
    USE cudafor
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: cru_state, irow_R
    !INTEGER, INTENT(INOUT) :: lcc_state, ryr_state
    INTEGER :: lcc_state, ryr_state
    
    lcc_state = (cru_state-1)/irow_R + 1
    ryr_state = cru_state - (lcc_state-1)*irow_R

  END SUBROUTINE state_breakdown_ker

  attributes(device) FUNCTION state_combine_single_ker(LCCcluster_state, RYRcluster_state, irow_R) RESULT(cru_state)
    USE cudafor
    INTEGER :: cru_state
    INTEGER, INTENT(IN) :: LCCcluster_state, RYRcluster_state
    INTEGER, INTENT(IN) :: irow_R
    
    cru_state = (LCCcluster_state-1) * irow_R + RYRcluster_state

    RETURN
  END FUNCTION state_combine_single_ker

  INTEGER FUNCTION print_mode_info() RESULT(res)
    !! print out information which mode is being used
    IMPLICIT NONE

    res = SUCCESS

#if _RUNMODE_ == _MODE_RESTING
    WRITE(log_string,*) "Resting"
    CALL write2log(log_string)
#elif _RUNMODE_ == _MODE_V_CLAMP
    WRITE(log_string,*) "V-clamped"
    CALL write2log(log_string)
    WRITE(log_string, "(A, F6.2, A, F6.2, A, A, F3.0)") "V_stim = ", V_stim, "mV; duration = ", (t_V_stim_end - t_V_stim_start) * 1d3, "ms", "; Hz=", Hz
    CALL write2log(log_string)

#elif _RUNMODE_ == _MODE_I_CLAMP
    WRITE(log_string,*) "I-clamped"
    CALL write2log(log_string)
    WRITE(log_string, "(A, F6.2, A, F6.2, A, A, F3.0)") "I_app = ", I_stim, "uA/cm^2; duration = ", (t_I_stim_end - t_I_stim_start)*1d3, "ms", "; Hz=", Hz
    CALL write2log(log_string)

#elif _RUNMODE_ == _MODE_CLOSEDCELL
    WRITE(log_string,*) "Closed cell (no ion in/out)"
    CALL write2log(log_string)

#else
    WRITE(log_string,*) "Error: unknown mode"
    CALL write2log(log_string)
    STOP
#endif    

#if _USE_CaDS_APPROX_ == _NO
    IF (dt_UB .GT. 1.0d-7) THEN
       WRITE(log_string,*) "WARNING: Need to reduce dt for exact calculation"
       CALL write2log(log_string)
    END IF
#endif
        
  END FUNCTION print_mode_info

  SUBROUTINE gpu_cleanup()
    IMPLICIT NONE
    INTEGER :: err
    
    err = cudaStreamDestroy(stream0)
    err = cudaStreamDestroy(stream1)
    
    err = cudaDeviceReset()
  END SUBROUTINE gpu_cleanup

  SUBROUTINE calc_new_RosDS(dt, RosProduced, RosEffluxed)
    IMPLICIT NONE
    REAL(KIND=dp),INTENT(IN) :: dt
    REAL(KIND=dp),INTENT(OUT) :: RosProduced, RosEffluxed
    REAL(KIND=dp) :: RosDS_new,RosRemoved
    RosDS_new  = 0.0d0
    CALL calcRosProd(RosProduced)
    
    CALL calcRosEfflux(RosEffluxed)
    CALL calcRosRemovedDS(RosDS,RosRemoved,dt)
    RosDS_new = RosDS+dt*inv_lambda_ds_Ros*(RosProduced-RosEffluxed-RosRemoved)
    RosDS = RosDS_new
    
  END SUBROUTINE calc_new_RosDS

  SUBROUTINE calc_new_RosI(dt,RosEffluxed)
    IMPLICIT NONE
    REAL(KIND=dp),INTENT(IN) :: dt, RosEffluxed
    REAL(KIND=dp) :: RosI_new, RosRemoved
    RosI_new = 0.0d0

    CALL calcRosRemovedI(RosI,RosRemoved,dt)
    RosI_new = RosI + dt*(effective_NSFU*RosEffluxed - RosRemoved)    
    RosI = RosI_new
  END SUBROUTINE calc_new_RosI

  SUBROUTINE calcRosProd(RosProd)
    IMPLICIT NONE
    REAL(KIND=dp), INTENT(OUT) :: RosProd
    REAL(KIND=dp) :: t_after_stretch, t_after_unstretch
    t_this_cycle = MOD((time-letgo-t_stretch_start),cycle_size)

    IF(time .LT. (letgo+t_stretch_start)) THEN
      RosProd = baseRos_prod
      kon_htrpn = temp_kon_htrpn
      kon_ltrpn = temp_kon_ltrpn
    !ELSEIF((time .LT. 4.0d1) .AND. (t_this_cycle .LT. stretch_duration) ) THEN
    ELSEIF(t_this_cycle .LT. stretch_duration) THEN
      t_after_stretch = t_this_cycle 
      !kon_htrpn = (2.0d1)*(1.2d0)  !(2.0d1) !*(1.1d0)!(2.37d0)*(1.1d0)  !(1.0d2)*(1.1d0)
      !koff_htrpn = (6.60d-2)*(0.7d0)
      !kon_ltrpn = (4.0d1)*(1.2d0) !(4.0d1) !*(1.1d0)!(4.00d0)*(1.1d0)  !(3.9d1)*(1.1d0)
      !koff_ltrpn = (4.0d1)*(0.7d0)
      kon_htrpn = temp_kon_htrpn*(1.2d0)
      kon_ltrpn = temp_kon_ltrpn*(1.2d0)
      IF(t_after_stretch .LT. prod_rise_time) THEN
        !RosProd = 3.25d-2*t_after_stretch !!m35,m35a
         !RosProd = 8.333d-3*t_after_stretch !!m41b
         !RosProd = 1.666d-3*t_after_stretch !!m41e
         !RosProd = RosProd_atstart_str+((2.d-4-RosProd_atstart_str)/1.2d-1)*t_after_stretch !!m41e
         RosProd = RosProd_atstart_str+((2.0d-4-RosProd_atstart_str)/1.2d-1)*t_after_stretch !!m41g,i,j

      ELSEIF(t_after_stretch .LT. 1.5d0) THEN
        !RosProd = 3.8d-4*t_after_stretch**(2.0d0)-2.0d-3*t_after_stretch+4.134528d-3 !m35,m35a
        !RosProd = 1.8d-4*t_after_stretch**(2.0d0)-6.1d-4*t_after_stretch+1.070608d-3 !m41b
        !RosProd = 3.6d-5*t_after_stretch**(2.0d0)-1.22d-4*t_after_stretch+2.14d-4 !m41e
        RosProd = 4.5d-5*t_after_stretch**(2.0d0)-1.525d-4*t_after_stretch+2.676d-4 !m41g,i,j
      !ELSEIF(t_after_stretch .LT. 5.0d0) THEN
      ELSEIF(t_after_stretch .LT. 4.0d0) THEN
        !RosProd = 2.7d-5*t_after_stretch**(2.0d0)-4.022d-4*t_after_stretch+2.350928d-3 !m35a
        !RosProd = 4.4d-6*t_after_stretch**(2.0d0)-8.0d-5*t_after_stretch+6.70708d-4 !m41b
        !RosProd = 8.8d-7*t_after_stretch**(2.0d0)-1.6d-5*t_after_stretch+1.34d-4 !m41e
        !RosProd = 1.1d-6*t_after_stretch**(2.0d0)-2.0d-5*t_after_stretch+1.67677d-4 !m41g,i
        RosProd = 1.25d-6*t_after_stretch**(2.0d0)-1.9d-5*t_after_stretch+1.658395d-4 !m41j
      ELSE
        !RosProd = 1.384571d-3*exp(-.07*t_after_stretch) !m35a
        !RosProd = 5.97069d-4*exp(-.09*t_after_stretch) !m41b
        !RosProd = 7.29262d-4*exp(-.13*t_after_stretch) !m41d
        !RosProd = 1.4585d-4*exp(-.13*t_after_stretch) !m41e
        !RosProd = 1.82315d-4*exp(-.13*t_after_stretch) !m41g
        !RosProd = 1.49267d-4*exp(-.09*t_after_stretch) !m41i
        RosProd = 1.57436d-4*exp(-.09*t_after_stretch) !m41j
      ENDIF
      RosProd_atstart_unstr = RosProd
    ELSE
      !IF(time .GT. 4.0d1) THEN
      !  t_after_unstretch = time - 4.0d1
      !ELSE
        t_after_unstretch = t_this_cycle - stretch_duration
      !ENDIF

      !kon_htrpn = 2.0d1*(7.5d-1) !1.0d2
      !koff_htrpn = 6.60d-2
      !kon_ltrpn = 4.0d1*(7.5d-1) !4.00d0 !3.9d1
     ! koff_ltrpn = 4.0d1
      kon_htrpn = temp_kon_htrpn
      kon_ltrpn = temp_kon_ltrpn

      IF(t_after_unstretch .LE. Ros_relaxTime) THEN
         RosProd =(RosProd_atstart_unstr/Ros_relaxTime)*(Ros_relaxTime-t_after_unstretch)
      ELSE
         RosProd = baseRos_prod   
         !RosProd = 0.0d0
      ENDIF 
      RosProd_atstart_str = RosProd
    ENDIF
  END SUBROUTINE calcRosProd

  SUBROUTINE calcRosEfflux(RosEffluxed)
    IMPLICIT NONE
    REAL(KIND=dp), INTENT(OUT) :: RosEffluxed

    RosEffluxed = Vefflux_Ros*(RosDS-RosI)
  END SUBROUTINE calcRosEfflux

  SUBROUTINE calcRosRemovedI(Ros,RosRemoved,dt)
    IMPLICIT NONE
    REAL(KIND=dp), INTENT(IN) :: Ros, dt
    REAL(KIND=dp), INTENT(OUT) :: RosRemoved

      !RosRemoved = v_RosRemoval_i*(Ros-baseRos_I)
      RosRemoved = v_RosRemoval_i*Ros-v_RosRemoval_base

  END SUBROUTINE calcRosRemovedI
  SUBROUTINE calcRosRemovedDS(Ros,RosRemoved,dt)
    IMPLICIT NONE
    REAL(KIND=dp), INTENT(IN) :: Ros, dt
    REAL(KIND=dp), INTENT(OUT) :: RosRemoved

      RosRemoved = v_RosRemoval_ds*Ros

  END SUBROUTINE calcRosRemovedDS
END MODULE ecc_gpu_utility


