 !!>****h*  NON-SPATIAL-CELL/data_global
!! SYNOPSIS
!!   This file contains data to be read from configuration file
!!   
!! CREATOR 
!!      Hoang-Trong Minh Tuan 
!! AFFILIATION
!!      George Mason University 
!! 
!! HISTORY
!!    * Rev: 0.1 (Feb, 2, 2010)
!!    * 0.2: add new data for spatial model, I-clamp & V-clamp (Hoang-Trong Minh Tuan, Oct, 18, 2010)
!!    * 0.3: change from using v_ryr_nj and P_dhpr_nj to prcent_ryr_nj and prcent_dhpr_nj (Tuan, Jul, 14, 2011)
!!    * 0.4: add 3 different testing conditions (orphan CaRUs, change capacitance, change I_LCC) (Tuan, Aug, 23, 2011)
!!    * 1.0: major revision - use 3 new parameter files (remove unused rate transition)
!!           - first 2 lines tell the Markov files to use for LCC and RyR (this will help switching from one model to another easier)
!!_****
MODULE data_global
  USE error_codes
  USE precisioninfo
  USE cuda_interop
  USE makeKs_utility
  IMPLICIT NONE
#include "macro.h"
  
#define _DEBUG_  0

  !! plan to use
  REAL(KIND=dp) :: i_1dhpr, i_1ryr
  REAL(KIND=dp) :: delta_Cs = 1.0d0, & ! the effect of Cesium
       K_Cs, Vhalf

  !! constant
  CHARACTER(12) :: PRECISION_OUTPUT = "(2x,1pe15.7)"
  INTEGER, PARAMETER :: MAX_CaRUs = 20000  ! number of CaRU in normal heart cells
  CHARACTER(20), PARAMETER :: first_beat = "first_beat"  !subfolder where to store data for first-beat
  CHARACTER(20), PARAMETER ::  FormatIn = "(D14.7, A35)", &
       FormatOut = "(A35, D14.7)"      !! matlab/Fortran parameter wrapper
  INTEGER, PARAMETER :: GROUP_LENGTH = 25
  CHARACTER(FILENAME_LENGTH), PARAMETER :: testdata_filename = "testdata_file.txt", &
       summary_filename = "summary_data.txt", &
       LCC_Markov_filename = "LCC_Markov.conf", &
       RyR_Markov_filename = "RYR_Markov.conf", &
       restart_filename = "restart_file.dat"
  INTEGER, PARAMETER :: MIN_SIZE_DIM = 30   ! minimum grid dimension in a CUDA kernel

  REAL(KIND=dp), PARAMETER :: zCa = 2.d0, zNa = 1.d0, zK = 1.d0   !! covalent numbers

  INTEGER, PARAMETER :: eta_pmca = 2  !Hill coefficient
  REAL(KIND=dp), PARAMETER :: Km_pmca = 0.5d0
  REAL(KIND=dp), PARAMETER :: eta_ncx = 0.35d0, &
       Ksat_ncx = 0.1d0, &  !saturation
       Km_ncxNa = 87.5d3, & !uM
       Km_ncxCa = 1.38d3 ! uM
  REAL(KIND=dp), PARAMETER :: cpl_serca = 2.0

  
#if _GPU_MODEL_ == _TESLA
  INTEGER, PARAMETER :: BLOCKSIZE_1D = 196
  INTEGER, PARAMETER :: BLOCK_DIMX = 32, BLOCK_DIMY = 8 ! 256 threads/block (Tesla)
#elif _GPU_MODEL_ == _FERMI  
  INTEGER, PARAMETER :: BLOCKSIZE_1D = 256
  !INTEGER, PARAMETER :: BLOCK_DIMX = 32, BLOCK_DIMY = 16 ! 512 threads/block (Fermi)
  !INTEGER, PARAMETER :: BLOCK_DIMX = 64, BLOCK_DIMY = 8 ! 512 threads/block (Fermi)
  INTEGER, PARAMETER :: BLOCK_DIMX = 32, BLOCK_DIMY = 8 ! 256 threads/block (Fermi)
  !INTEGER, PARAMETER :: BLOCK_DIMX = 16, BLOCK_DIMY = 16 ! 256 threads/block (Fermi)
#endif
  !! GPU-related info
  INTEGER :: gpuID  ! ID of GPU (Fermi, Tesla)
  INTEGER :: blocksize  ! ThreadsPerBlock (assuming 1D)
  INTEGER :: gridSize   ! BlocksPerGrid (assuming 1D)
      
  !!! unconventional testing
  REAL(KIND=dp) :: Ej_small, Ej_orphan
  INTEGER :: NSFU_bigRyR, NSFU_smallRyR, NSFU_normal, NSFU_orphan, FULL_NSFU, N_R_small
  REAL(KIND=dp) :: V_ds_orphan, inv_lambda_ds_orphan
  INTEGER :: is_dualRyRcluster
  REAL(KIND=dp) :: per_largeRyR

  !! global data
  CHARACTER(FILENAME_LENGTH) :: conf_file, prog_name
  INTEGER :: ryr_initialstate, lcc_initialstate
  !sys_command   = contain the command to be executed from the terminal
  CHARACTER(MAX_STR_LEN) :: sys_command = ""
  REAL(KIND=dp) :: time ! keep the current time in simulation
  LOGICAL :: is_EndSim
  REAL(KIND=dp) :: real_timestart = 0D0  !being used only when restart-file is read
  real(KIND=dp) :: V_myo_T_true
  
  !+ computed variables
  REAL(KIND=dp) :: RT_F, V_ds, V_jsr, &
       lambda_ds, lambda_jsr, lambda_nsr &
       , v_refill, v_efflux, v_ryr  &
       , Acap, Am_2FVmyo, Am_FVmyo, &
       Am_2FV, Am_FV, inv_2FV
  REAL(KIND=dp) :: RT_FzCa, F_RT, invKp_myo, invKp_nsr
  REAL(KIND=dp) :: P_dhpr, F2zCa2_RT, FzCa_RT, inv_2FVmyo
  REAL(KIND=dp) :: t_V_stim_end, t_I_stim_end
  REAL(KIND=dp) :: Vol_ref   ! volume reference on which concentrations, transfer rates are defined
  REAL(KIND=dp) :: v_ryr_nj_T, P_dhpr_nj
  REAL(KIND=dp) :: inv_lambda_ds, inv_lambda_jsr, inv_lambda_nsr
  INTEGER :: TOTAL_GRID_ELEL  !keep the "effector" total number of grid element in 3D cells
  REAL(KIND=dp) :: CaF_infty
  REAL(KIND=dp) :: dx2, dy2, dz2
  REAL(KIND=dp) :: invzCa
  INTEGER :: MAXX, MAXY, MAXZ ! maximum number of grid cells

  !!Markov-model data
  !REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: ryr_KCal_act, ryr_Kfree
  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: ryr_KCal_act_mode1, ryr_Kfree_mode1, ryr_KCal_act_mode2, ryr_Kfree_mode2, &
                                              ryr_Kfree, ryr_KRos

  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: lcc_Kv_act, lcc_Kv_inact, lcc_KCal_inact, lcc_Kfree
  INTEGER, DIMENSION(:), ALLOCATABLE :: iU_R, iU_L



  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!Data from configuration file
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  INTEGER :: NSFU, N_L, N_R, mL, mR   !channels number + state number
  REAL(KIND=dp) :: V_rest, V_stim, I_stim, t_start, t_end, dt_LB, dt_UB  !Voltage + Current + time
  INTEGER :: write_restartfile
  INTEGER :: write_firstbeat  
  REAL(KIND=dp) :: duration_write_lastbeats
  REAL(KIND=dp) :: t_V_stim_start, t_I_stim_start,  t_V_stim_duration, t_I_stim_duration, Hz !when to stimulate + frequency

  REAL(KIND=dp) :: write_interval !(sec)
  REAL(KIND=dp) :: Na_o, Na_i, Ca_o, K_o, K_i, Csc
  !+ calcium ionic concentrations
  REAL(KIND=dp) :: Ca_myo, Ca_nsr
  REAL(KIND=dp) :: zF, zR, zT
  REAL(KIND=dp) :: eta_RyR, eta_LCC, Ej, Ecc, Eoo
  !+ volume
  REAL(KIND=dp) :: V_myo_T, V_nsr_T, V_jsr_T, V_ds_T

  !+ buffering
  REAL(KIND=dp) :: beta_ds, beta_nsr, B_myoT, B_jsrT, B_htrpnT, B_ltrpnT
  REAL(KIND=dp) :: Km_myo, Km_jsr, kon_htrpn, koff_htrpn, kon_ltrpn, koff_ltrpn
  !+ rate
  REAL(KIND=dp) :: v_refill_T, v_efflux_T, v_leak, &
       v_ryr_T, P_dhpr_T,  prcent_ryr_nj, prcent_dhpr_nj
  REAL(KIND=dp) :: Ap, Kp_myo, Kp_nsr, Am
  REAL(KIND=dp) :: g_bCa, g_Na, g_K1, g_Kss, g_Ktof, &
       g_Ktos, g_bNa, g_bK
  !+ maximum current
  REAL(KIND=dp) :: I_pmca_bar, I_ncx_bar, I_NaK_bar
  !+ 
  REAL(KIND=dp) :: k_jsr0, k_jsr1
  REAL(KIND=dp) :: V_theta2, V_sigma2, vtheta_dhpr, sigma_v_dhpr
  !+ for random number generator
  INTEGER :: rand_seed, fixed_seed
  INTEGER :: TRACK_LENGTH 
  !+ cell geometry
  REAL(KIND=dp) :: cell_x_len, cell_y_len, cell_z_len
  !+ size of grid element
  REAL(KIND=dp) :: dX_, dY_, dZ_   
  !+ diffusion coefficient
  REAL(KIND=dp) :: Dmyo_x, Dmyo_y, Dmyo_z   
  REAL(KIND=dp) :: Dnsr_x, Dnsr_y, Dnsr_z
  !+ dye data
  REAL(KIND=dp) :: Ddye
  REAL(KIND=dp) :: F_T, kon_F, koff_F
  !+ inter CRU distance
  REAL(KIND=dp) :: lspacing, tspacing, dist2membraneX, dist2membraneYZ

  !!!!!!!!!!!!!!!!!!
  PRIVATE :: read_param
  SAVE  
CONTAINS
  !!>****f* data_global/read_param
  !! SYNOPSIS  
  !!   Read in parameters from shared fortran-Matlab configuration
  !!   file:
  !!   * first line is omitted
  !!   * skip blank lines
  !!   * comments start WITH # at the first column
  !!   Return 0 (success); otherwise (error)
  !! INPUTS
  !!   conf_file : name of the configuration file (optional)
  !!               default is parameters.conf
  !! OUTPUTS
  !! RETURN
  !!   SUCCESS or an error code
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  FUNCTION read_param(conf_file) RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    CHARACTER(*), OPTIONAL :: conf_file
    !!local
    REAL(KIND=dp), DIMENSION(150) :: var_vect
    CHARACTER(40) :: var_name
    INTEGER :: ii
    CHARACTER :: dummy_char
    REAL(KIND=dp), DIMENSION(GROUP_LENGTH) :: temp_data
    INTEGER :: fileID, end_idx, error_code
    CHARACTER(FILENAME_LENGTH) :: filename
    LOGICAL :: fexist
    REAL(KIND=dp) :: Vm, prcent_current_change, prcent_cap_change, prcent_orphan_sfu = 0d0
    REAL(KIND=dp) :: dummy, factor_Vds, K1_reduction, Kf_reduction, Ks_reduction, &
         factor_NCX, SERCA_reduction, Cm_reduction, &
         Ejorphan_reduction, Ejnormal_reduction, factor_kR12, factor_P_dhpr
    
    res = SUCCESS

    IF (PRESENT(conf_file)) THEN
       filename = conf_file
    ELSE
       filename = "./parameters.conf"
    END IF
    WRITE(log_string, "(A,A)") "parameter file: ", TRIM(filename)
    CALL write2log(log_string)

    INQUIRE(FILE=filename, EXIST=fexist)
    IF (.NOT. fexist) THEN
       res  = ERR_FILE_NOT_EXIST
       WRITE(error_message, *), "File " // TRIM(filename) // " not exist"
       RETURN 
    END IF
    IF (TRIM(where2save) .NE. ".") THEN
       sys_command = 'cp -rL ' // TRIM(conf_file) // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
       sys_command = 'cp -rL ' // TRIM(prog_name) // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
       sys_command = 'cp -rL ' // "src/macro.h" // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
    END IF

    
    fileID = get_file_ID()
    OPEN(unit=fileID, file=filename, ACTION='READ')
    !! skip the first line
    READ(fileID, '(A1)', IOSTAT=error_code) dummy_char 

    ii = 1
    DO WHILE (error_code .eq. 0)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (error_code .EQ. 0 .AND. dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,FormatIn, IOSTAT=error_code) var_vect(ii), var_name
          ii = ii + 1
          var_name = ""
       END IF
    END DO
    CLOSE(unit = fileID)
    
    end_idx = 1

    !! data-related info
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    fixed_seed = int(temp_data(1))
    rand_seed = int(temp_data(2))
    TRACK_LENGTH = int(temp_data(3))
    end_idx = end_idx + 3
#if _DEBUG_ == 1    
    WRITE(*,*) "fixed_seed", fixed_seed
    WRITE(*,*) "rand_seed", rand_seed
    WRITE(*,*) "TRACK_LENGTH", TRACK_LENGTH
#endif
    
    !! release site configuration
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    NSFU = int(temp_data(1))
    N_R  = int(temp_data(2))
    N_L  = int(temp_data(3))
    end_idx = end_idx + 3
#if _DEBUG_ == 1
    WRITE(*,*) "NSFU, N_R, N_L"
    WRITE(*,*) NSFU, N_R, N_L
#endif
    
    !! simulation settings
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    V_rest  = temp_data(1)
    V_stim  = temp_data(2)
    I_stim  = temp_data(3)  ! [uA/cm^2]
    t_start = temp_data(4)
    t_end   = temp_data(5)
    write_firstbeat = INT(temp_data(6))
    duration_write_lastbeats = temp_data(7)
    write_restartfile = INT(temp_data(8))
    dt_LB  = temp_data(9)
    dt_UB  = temp_data(10)
    t_V_stim_start = temp_data(11)
    t_V_stim_duration = temp_data(12)
    t_I_stim_start = temp_data(13)
    t_I_stim_duration = temp_data(14)
    Hz  = temp_data(15)
    write_interval = temp_data(16)  ! (sec)
    end_idx = end_idx + 16

    t_V_stim_end = t_V_stim_duration + t_V_stim_start
    t_I_stim_end = t_I_stim_duration + t_I_stim_start

    IF (t_V_stim_end .GE. 1d0/Hz+t_V_stim_start) THEN
       PRINT *, "V_stim start-end:", t_V_stim_start, t_V_stim_end
       PRINT *, "WARNING: Voltage stimulus is ON all the time"
    END IF
    IF (t_I_stim_end .GE. 1d0/Hz+t_I_stim_start) THEN
       PRINT *, "I_stim start-end:", t_I_stim_start, t_I_stim_end
       PRINT *, "WARNING: Current stimulus is ON all the time"
    END IF
#if _DEBUG_ == 1
    WRITE (*,*) "V_rest, V_stim, I_stim, t_start, t_end, dt_LB, dt_UB, t_V_stim_start" // &
         "t_V_stim_end, t_I_stim_start, t_I_stim_end, Hz, write_interval"
    WRITE(*,"(20(D14.5))") V_rest, V_stim, I_stim, t_start, t_end, dt_LB, dt_UB, t_V_stim_start, &
         t_V_stim_end, t_I_stim_start, t_I_stim_end, Hz
    WRITE(*, *) write_interval 
#endif
    
    !! ion concentrations
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    Na_i = temp_data(1)      !uM
    K_i  = temp_data(2)
    Ca_myo = temp_data(3)
    Ca_nsr = temp_data(4)
    end_idx = end_idx+4
#if _DEBUG_ == 1
    WRITE(*,*) "Na_i, K_i, Ca_myo, Ca_nsr"
    WRITE(*,"(20(D14.5))") Na_i, K_i, Ca_myo, Ca_nsr
#endif
    
    !! constants
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    zF = temp_data(1)         !Faraday's constant, C/mol
    zR = temp_data(2)         !universal gas constant, mJ mol-1 K-1
    zT = temp_data(3)         !temperature, K
    Na_o = temp_data(4)      !uM
    Ca_o = temp_data(5)      !uM
    K_o  = temp_data(6)
    end_idx = end_idx + 6
#if _DEBUG_ == 1
    WRITE(*,*) "zF, zR, zT, Na_o, Ca_o, K_o"
    WRITE(*,"(10(D14.5))") zF, zR, zT, Na_o, Ca_o, K_o
#endif
    
    !! cell geometry variables (volumes) !pL
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    V_myo_T = temp_data(1)
    V_nsr_T = temp_data(2)
    V_jsr_T = temp_data(3)
    V_ds_T = temp_data(4)
    end_idx = end_idx + 4

#if _DEBUG_ == 1
    WRITE(*,*) "V_myo_T, V_nsr_T, V_jsr_T, V_ds_T"
    WRITE(*,"(D14.5)") V_myo_T, V_nsr_T, V_jsr_T, V_ds_T
#endif    
    !! buffering parameters
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    beta_ds = temp_data(1)
    beta_nsr = temp_data(2)
    B_myoT = temp_data(3)
    B_jsrT = temp_data(4)
    Km_myo = temp_data(5)
    Km_jsr = temp_data(6)
    B_htrpnT = temp_data(7)
    kon_htrpn = temp_data(8)
    koff_htrpn = temp_data(9)
    B_ltrpnT = temp_data(10)
    kon_ltrpn = temp_data(11)
    koff_ltrpn = temp_data(12)
    end_idx = end_idx + 12
#if _DEBUG_ == 1
    WRITE(*,*) "beta_ds, beta_nsr, B_myoT, B_jsrT, Km_myo, Km_jsr, B_htrpnT, kon_htrpn, koff_htrpn"
    WRITE(*,"(20(D14.5))") beta_ds, beta_nsr, B_myoT, B_jsrT, Km_myo, Km_jsr, B_htrpnT, kon_htrpn, koff_htrpn
#endif
    
    !! diffusion parameters
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    v_refill_T = temp_data(1)
    v_efflux_T = temp_data(2)
    v_leak     = temp_data(3)
    v_ryr_T = temp_data(4)
    P_dhpr_T   = temp_data(5)
    prcent_ryr_nj = temp_data(6)
    prcent_dhpr_nj = temp_data(7)
    end_idx = end_idx + 7
#if _DEBUG_ == 1
    WRITE(*, *) "v_refill_T, v_efflux_T, v_leak, v_ryr_T, P_dhpr_T, prcent_ryr_nj, prcent_dhpr_nj"
    WRITE(*,"(20(D14.5))") v_refill_T, v_efflux_T, v_leak, v_ryr_T, P_dhpr_T, prcent_ryr_nj, prcent_dhpr_nj
#endif    

    !! SERCA variables
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    Ap = temp_data(1)           !Ap = SERCA pumps concentration (uM)
    Kp_myo = temp_data(2)       
    Kp_nsr = temp_data(3)       
    end_idx = end_idx + 3
#if _DEBUG_ == 1
    WRITE(*,*) "Ap, Kp_myo, Kp_nsr"
    WRITE(*, "(20(D14.5))") Ap, Kp_myo, Kp_nsr
#endif
    
    !! membrane currents variables
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    g_bCa = temp_data(1)
    g_Na = temp_data(2)
    g_K1 = temp_data(3)
    g_Kss = temp_data(4)
    g_Ktof = temp_data(5)
    g_Ktos = temp_data(6)
    g_bNa = temp_data(7)
    g_bK = temp_data(8)
    Am = temp_data(9)   !cm^2
    Csc   = temp_data(10)      !specific capacitance uF/cm^2
    end_idx = end_idx + 10
#if _DEBUG_ == 1
    WRITE(*,*) "g_bca, g_Na, g_K1, g_Kss, g_Ktof, g_Ktos, g_bNa, g_bK, Am, Csc"
    WRITE(*, "(20(D14.5))") g_bca, g_Na, g_K1, g_Kss, g_Ktof, g_Ktos, g_bNa, g_bK, Am, Csc
#endif
    
    !! PMCA and NCX rates
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    I_pmca_bar = temp_data(1)   !uA/uF
    I_ncx_bar = temp_data(2)    !uA/uF
    I_NaK_bar = temp_data(3)
    end_idx = end_idx + 3
#if _DEBUG_ == 1
    WRITE(*,*) "I_pmca_bar, I_ncx_bar, I_NaK_bar"
    WRITE(*, "(20(D14.5))") I_pmca_bar, I_ncx_bar, I_NaK_bar
#endif
    
    !! channels parameters: rate constants, constants
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    eta_RyR = temp_data(1)
    eta_LCC = temp_data(2)
    Ej = temp_data(3)
    Ecc = temp_data(4)
    Eoo = temp_data(5)
    k_jsr0 = temp_data(6)
    k_jsr1 = temp_data(7)
    end_idx = end_idx + 7
#if _DEBUG_ == 1
    WRITE(*,*) "eta_RyR, eta_LCC, Ej, Ecc, Eoo, k_jsr0, k_jsr1"
    WRITE(*,"(20(D14.5))") eta_RyR, eta_LCC, Ej, Ecc, Eoo, k_jsr0, k_jsr1
    STOP
#endif
    
    !! LCC kinetics parameters
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    Vtheta_dhpr = temp_data(1)  ! mV
    sigma_v_dhpr = temp_data(2) ! mV
    V_theta2 = temp_data(3)
    V_sigma2 = temp_data(4)
    end_idx = end_idx + 4
#if _DEBUG_ == 1
    WRITE(*,*) "Vtheta_dhpr, sigma_v_dhpr, V_theta2, V_sigma2"
    WRITE(*,"(20(D14.5))") Vtheta_dhpr, sigma_v_dhpr, V_theta2, V_sigma2
#endif

    !! Spatial info (grid elements)
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    cell_x_len = temp_data(1)
    cell_y_len = temp_data(2)
    cell_z_len = temp_data(3)
    dX_ = temp_data(4)
    dY_ = temp_data(5)
    dZ_ = temp_data(6)
    end_idx = end_idx + 6
#if _DEBUG_ == 1
    WRITE(*, *) "Cell dimension (true): "
    WRITE(*,"(20(D14.5))") cell_x_len, cell_y_len, cell_z_len 
    WRITE(*, *) "Grid size"
    WRITE(*,"(20(D14.5))") dX_, dY_, dZ_    
#endif
    
    !! Spatial info (diffusion constant)
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    Dmyo_x = temp_data(1)
    Dmyo_y = temp_data(2)
    Dmyo_z = temp_data(3)
    Dnsr_x = temp_data(4)
    Dnsr_y = temp_data(5)
    Dnsr_z = temp_data(6)
    Ddye   = temp_data(7)  ! assume isotropic diffusion of dye
    F_T    = temp_data(8)
    kon_F  = temp_data(9)
    koff_F = temp_data(10)    
    end_idx = end_idx + 10
#if _DEBUG_ == 1
    WRITE(*,*) "Dmyo(x,y,z), Dnsr(x,y,z)"
    WRITE(*,"(20(D14.5))") Dmyo_x, Dmyo_y, Dmyo_z, Dnsr_x, Dnsr_y, Dnsr_z
#endif
    
    !! Spatial info (CaRU distances)
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    lspacing = temp_data(1)
    tspacing = temp_data(2)
    end_idx = end_idx + 2

    !! SINGLE CURRENTS
    temp_data = var_vect(end_idx:end_idx+GROUP_LENGTH)
    i_1ryr = temp_data(1)
    end_idx = end_idx + 1

    !! computed constants
    RT_F = zR*zT/zF
    F_RT = zF/(zR*zT)
    RT_FzCa = RT_F/zCa
    !!++permeability at each LCC channel
    !convert from 10^⁻3 cm^3/s to L/s    
    P_dhpr_T = P_dhpr_T * 1d-6  ! [L/s]
    
    P_dhpr = P_dhpr_T/(dfloat(NSFU*N_L))  ! [L/s]
#if _USE_SCALING_MODE_ == _NO
    P_dhpr = P_dhpr_T/(MAX_CaRUs*N_L)  ! [L/s]
#endif    
    
    !!TUAN PLAN: replace P_dhpr_T by i_1dhpr from parameter file
    !P_dhpr_T = i_1dhpr / ((zCa**2*zF*Vm)/RTF * ((0.1d0*EXP(2*Vm/RTF)-0.341 * &
    !     Ca_o)/(EXP(2*Vm/RTF)-1))) * NSFU * N_L
    i_1dhpr = -0.15d0 !pA

    NSFU_orphan = 0
    NSFU_normal = NSFU
    inv_lambda_ds_orphan = inv_lambda_ds
    
    !!++ volume [pL]
    V_ds = V_ds_T/dfloat(NSFU)
    V_jsr = V_jsr_T/dfloat(NSFU)
#if _USE_SCALING_MODE_ == _NO
    V_ds = V_ds_T/dfloat(MAX_CaRUs)
    V_jsr = V_jsr_T/dfloat(MAX_CaRUs)
#endif    

    Vol_ref = V_myo_T   ! use total myoplasmic volume as volume reference

    !!read in Markov-model
    !!res = read_RyR_Markov(RyR_Markov_filename, ryr_KCal_act, ryr_Kfree, mR, iU_R, ryr_initialstate)
    res = read_RyR_Markov(RyR_Markov_filename, ryr_KCal_act_mode1, ryr_Kfree_mode1, &
	ryr_KCal_act_mode2, ryr_Kfree_mode2, ryr_Kfree, ryr_KRos, mR, iU_R, ryr_initialstate)

    IF (res .NE. SUCCESS) CALL cleanup()  !ERROR OCCUR
    res = read_LCC_Markov(LCC_Markov_filename, lcc_Kv_act, lcc_Kv_inact, &
         lcc_KCal_inact, lcc_Kfree, mL, iU_L, lcc_initialstate)
    IF (res .NE. SUCCESS) CALL cleanup()  !ERROR OCCUR    
    !! copy Markov model
    IF (TRIM(where2save) .NE. '.') THEN
       sys_command = "cp -fL " // TRIM(RyR_Markov_filename) // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
       sys_command = "cp -fL " // TRIM(LCC_Markov_filename) // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
    ENDIF

#if _USE_TEST_FILE_ == _YES
    !! copy testdata file
    IF (TRIM(where2save) .NE. '.') THEN
       sys_command = "cp -fL " // TRIM(testdata_filename) // '  ' // &
            TRIM(where2save) // '/'
       CALL system(sys_command)
    ENDIF
    res = read_testdata(var_vect, testdata_filename)
    IF (res .NE. SUCCESS) RETURN  !ERROR OCCUR
    factor_Vds = var_vect(1)
    prcent_orphan_sfu = var_vect(2)
    K1_reduction = var_vect(3)
    Ks_reduction = var_vect(4)
    Kf_reduction = var_vect(5)
    factor_NCX = var_vect(6)
    SERCA_reduction = var_vect(7)
    Cm_reduction = var_vect(8)
    Ejorphan_reduction = var_vect(9)
    Ejnormal_reduction = var_vect(10)
    factor_kR12 = var_vect(11)
    factor_P_dhpr = var_vect(12)

    P_dhpr = P_dhpr * factor_P_dhpr ! LCC current change
    WRITE(log_string, *) "LCC current change: ", factor_P_dhpr "-fold"
    CALL write2log(log_string)
    ! reduce membrane area
    Am = Am *(1d0 - Cm_reduction) 
    WRITE(log_string, *) "percent capacitance reduce: ", Cm_reduction
    CALL write2log(log_string)
    
    NSFU_orphan = INT(NSFU * (prcent_orphan_sfu))
    NSFU_normal = NSFU - NSFU_orphan
    V_ds_orphan = V_ds * factor_Vds
    P_dhpr_T = P_dhpr * DBLE(NSFU * N_L)
 #if _USE_SCALING_MODE_ == _NO
    P_dhpr_T = P_dhpr * DBLE(MAX_CaRUs * N_L)
 #endif
    
    inv_lambda_ds_orphan = Vol_ref/V_ds_orphan
    Ej_orphan = Ej * (1D0 - Ejorphan_reduction)
    
    WRITE(log_string, *) "percent orphan CaRUs (see large Vds): ", prcent_orphan_sfu
    call write2log(log_string)
    WRITE(log_string, *) "scaling factor Vds: ", factor_Vds
    call write2log(log_string)
    WRITE(log_string, *) "Ej reduction: ", Ejorphan_reduction
    call write2log(log_string)
    WRITE(log_string, *) "NSFU normal = ", NSFU_normal
    call write2log(log_string)
    WRITE(log_string, *) "NSFU orphan = ", NSFU_orphan
    call write2log(log_string)

    g_k1 = g_k1 * ( 1d0 - K1_reduction)
    g_ktof = g_ktof  * ( 1d0 - Kf_reduction)
    g_ktos = g_ktos * ( 1d0 - Ks_reduction)

    Ap = Ap * (1d0 - SERCA_reduction)
    I_ncx_bar = I_ncx_bar * factor_NCX
    Ej = Ej * (1d0 - Ejnormal_reduction)
    kR12 = kR12 * factor_kR12
    
    WRITE(log_string, *) "HF condition:"
    CALL write2log(log_string)

#endif
        
    !!++volume ratio
    lambda_ds = V_ds/Vol_ref
    lambda_jsr = V_jsr/Vol_ref
    lambda_nsr = V_nsr_T/Vol_ref
    ! for computing efficient, use inv_lambda values
    inv_lambda_ds = Vol_ref/V_ds
    inv_lambda_jsr = Vol_ref/V_jsr
    inv_lambda_nsr = Vol_ref/V_nsr_T    

    !! for small size RyR cluster
    is_dualRyRcluster = 0  !! 0 or 1
    IF (is_dualRyRcluster .EQ. 1) THEN
       N_R_small  = 16
       per_largeRyR = 50  !percent small RyR clusters  %
       Ej_small = 0.2  !!(a* for 4x4 RyR cluster)

       NSFU_bigRyR = INT(REAL(NSFU * per_largeRyR)/100)
       NSFU_smallRyR = NSFU - NSFU_bigRyR
       !v_ryr_smallRyR = v_ryr_T/dfloat(NSFU*N_R_small)
    END IF
    
    ! work for non-spatial only
    v_refill = v_refill_T/dfloat(NSFU)
    v_efflux = v_efflux_T/dfloat(NSFU)
    v_ryr = v_ryr_T/dfloat(NSFU*N_R)
 #if _USE_SCALING_MODE_ == _NO
    v_refill = v_refill_T/dfloat(MAX_CaRUs)
    v_efflux = v_efflux_T/dfloat(MAX_CaRUs)
    v_ryr = v_ryr_T/dfloat(MAX_CaRUs*N_R)    
 #endif
    
    v_ryr_nj_T = v_ryr_T*prcent_ryr_nj
    P_dhpr_nj = P_dhpr_T*prcent_dhpr_nj

    !!++may be used by EC-coupling model
    F2zCa2_RT=(zF**2)*(zCa**2)/(zR*zT) !  C/(mol.mV)
    FzCa_RT=(zF*zCa)/(zR*zT)     ! 1/mV
    invKp_myo = 1.0d0/Kp_myo
    invKp_nsr = 1.0d0/Kp_nsr

    ! TUAN PLAN: plan to remove Ej from parameter file and use
    !Ej = (2*4 + 3 * 4 * (SQRT(REAL(N_R)) - 2) + 4 * (SQRT(REAL(N_R))-2)**2)/ (N_R*(N_R-1))
    !PRINT *, "Ej = ", Ej
    ! this requires reference volume must be in pL
    inv_2FVmyo = 1.0d12/(zCa*zF*Vol_ref)  !mol/(C.L) !use when volume in pL

    
    IF (ABS(V_stim) .LT. EPS)  V_stim = 0.0001d0
    
    !!++may be used by spatial whole-cell model
    invzCa = 1.0d0/zCa
    
    dx2 = dX_**2
    dy2 = dY_**2
    dz2 = dZ_**2

    MAXX = INT(cell_x_len/dX_)
    MAXY = INT(cell_y_len/dY_)
    MAXZ = INT(cell_z_len/dZ_)
    TOTAL_GRID_ELEL = MAXX*MAXY*MAXZ
    
#if _TEST_ORPHAN_CRU_ == _YES && _TEST_NEW_ORPHAN_CRU_ == _NO
    !!TUAN TEST
    ! orphan CaRUs x%
    ! by moving L-type from subspace to non-junctional
    ! and RyRs even though in groups of 49 but see the Ca_myo, not Ca_ds
    !prcent_orphan_sfu = 10d0/100
    prcent_orphan_sfu = 25d0/100

    P_dhpr_nj = P_dhpr_nj + P_dhpr * DBLE(N_L) * NSFU * prcent_orphan_sfu
 #if _USE_SCALING_MODE_ == _NO
    P_dhpr_nj = P_dhpr_nj + P_dhpr * DBLE(N_L) * MAX_CaRUs * prcent_orphan_sfu
 #endif
    
    FULL_NSFU = NSFU
    NSFU = INT(NSFU * (1d0-prcent_orphan_sfu))
    
    P_dhpr_T = P_dhpr * DBLE(NSFU * N_L)
    WRITE(log_string, *) "percent orphan CaRUs (see Cmyo): ", prcent_orphan_sfu
    CALL write2log(log_string)
    WRITE(log_string, *) "NSFU new = ", NSFU
    CALL write2log(log_string)
#endif

    Acap = Am*Csc  ! [uF]
    !! use in calcium current calculation (use when volume reference in pL)
    Am_2FVmyo = (Am*1.0d12)/(zCa*zF*Vol_ref)  ! [mol.cm^2/(C.L)]
    !! use in Na, K current calculation
    Am_FVmyo = (Am*1.0d12)/(zF*Vol_ref)  ! [mol.cm^2/(C.L)]

        
#if _DEBUG_ == 1
    STOP
#endif
  END FUNCTION  read_param

  !!>****f* data_global/check_spatial_info
  !! SYNOPSIS  
  !!   Check if data for spatial cell model is legal
  !!   Return 0 (success); otherwise (error)
  !! INPUTS 
  !! OUTPUTS
  !! RETURN
  !!  SUCCESS = okay
  !!  otherwise = error
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  FUNCTION check_spatial_info() RESULT(res)
    INTEGER :: res
    
    res = SUCCESS

    !! check input
    IF (ABS(dX_ - dY_) .GT. EPS .OR. ABS(dX_ - dZ_) .GT. EPS .OR. ABS(dY_ - dZ_) .GT. EPS) THEN
       WRITE(error_message, *) "Current limit: grid size must be of equal in all dimension"
       res = ERR_SPATIAL_GRID_SIZE_NOT_UNIFORM
    END IF
    IF (ABS(dX_ - 0.2d0) .GT. EPS .AND. ABS(dX_ - 0.1d0) .GT. EPS) THEN
       WRITE(error_message,*) "Current limit: grid size must be either 0.2um or 0.1 um"
       res = ERR_SPATIAL_GRID_SIZE_NOT_VALID
    END IF
    !!check requirement
    IF (MAXX .LE. MIN_SIZE_DIM .OR. MAXY .LE. MIN_SIZE_DIM .OR. MAXZ .LE. MIN_SIZE_DIM) THEN
       WRITE(error_message, *) "Each dimension should be larger than ", MIN_SIZE_DIM
       res = ERR_SPATIAL_MINIMUM_DIM_SIZE
    END IF
  END FUNCTION check_spatial_info

  !!>****f* data_global/prog_help
  !! SYNOPSIS  
  !!   Print out help information of how to use the program
  !! INPUTS 
  !! OUTPUTS 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  SUBROUTINE print_prog_help(prog)
    CHARACTER(*) :: prog
    WRITE(*,*) prog, " <arg1> <arg2> <arg3> <arg4>"
    WRITE(*,*) "<arg1> is the name of the configuration file"
    WRITE(*,*) "  use ... to skip the argument (default: parameters.conf)"
    WRITE(*,*) "<arg2> is the gpu ID"
    WRITE(*,*) "  use ... to skip the argument (default: device with maximum GFLOPS)"
    WRITE(*,*) "<arg3> is the location where to save the output files"
    WRITE(*,*) "  use ... to skip the argument (default: . )"
    WRITE(*,*) "<arg4> tells this is a continuing simulation, and the value is the file"
    WRITE(*,*) "   containing the necessary data for the continuous simulation."
    WRITE(*,*) " NOTE: <arg4> has not been implemented"
    WRITE(*,*) " "
    WRITE(*,*) "The name of single-output file is automatically generated based on DATETIME()"
    WRITE(*,*) "The name of multiple-output file are based on \"data-\" + increasing number"
    STOP
  END SUBROUTINE print_prog_help

  !!>****f* data_global/prog_get_info
  !! SYNOPSIS  
  !!   Read command-line arguments
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN
  !!    0 = success
  !!    otherwise = error
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  FUNCTION  prog_get_info() RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    !local
    INTEGER :: num_args
    CHARACTER(FILENAME_LENGTH) :: dest, strgpuID
    LOGICAL :: is_exist
    
    res = SUCCESS

    where2save = "."

    num_args = COMMAND_ARGUMENT_COUNT()
    CALL GET_COMMAND_ARGUMENT(0, prog_name)
    full_command = TRIM(full_command) // ' ' // trim(prog_name)
        
    IF (num_args .GE. 3) THEN
       CALL GET_COMMAND_ARGUMENT(3, dest)
       where2save = TRIM(dest)
    END IF
    INQUIRE(FILE = TRIM(where2save), EXIST = is_exist)
    IF (.NOT. is_exist) THEN 
       sys_command = 'mkdir ' // TRIM(where2save)
       CALL system(sys_command)
       INQUIRE(FILE = TRIM(where2save) // "/" //TRIM(first_beat), EXIST = is_exist)
       IF (.NOT. is_exist) THEN
          sys_command = 'mkdir ' // TRIM(where2save) // "/" //TRIM(first_beat)
          CALL system(sys_command)
       ENDIF
    END IF
    IF (num_args .EQ. 0) THEN
       res = read_param()
    ELSE IF (num_args .EQ. 1) THEN
       CALL GET_COMMAND_ARGUMENT(1, conf_file)
       
       IF ("-h" .EQ. TRIM(conf_file)) THEN
          CALL print_prog_help(trim(prog_name))
       ELSE
          res = read_param(TRIM(conf_file))
       END IF
    ELSEIF (num_args .GE. 2) THEN
       CALL GET_COMMAND_ARGUMENT(1, conf_file)
       CALL GET_COMMAND_ARGUMENT(2, strgpuID)
       full_command = TRIM(full_command) // ' ' // trim(conf_file)
       full_command = TRIM(full_command) // ' ' // trim(strgpuID)

       IF ("..." .EQ. TRIM(conf_file)) THEN
          res = read_param()
       ELSE
          res = read_param(TRIM(conf_file))
       END IF

       IF ("..." .EQ. TRIM(strgpuID)) THEN
          gpuID = -1
       ELSE
          READ(strgpuID, '(I10)') gpuID
       END IF
    END IF

    IF (num_args .GE. 3) THEN
       full_command = TRIM(full_command) // ' ' // TRIM(dest)       
       IF (num_args .EQ. 4) THEN
          CALL GET_COMMAND_ARGUMENT(4, dest)
          restart_file = TRIM(dest)
          full_command = TRIM(full_command) // ' ' // TRIM(dest)
       END IF
    ENDIF

    WRITE(log_string, *) TRIM(full_command) 
    CALL write2log(log_string)
    
  END FUNCTION prog_get_info
  
  !!>****f* data_global/check_ecc_data
  !! SYNOPSIS  
  !!   Validate data from configuration file
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN
  !!    0 = success
  !!    otherwise = error
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  FUNCTION check_ecc_data() RESULT(res)
    INTEGER :: res

    res = SUCCESS
!!$    IF (dt_UB - 1.d-6 .GE. EPS) THEN
!!$       PRINT *, "dt_UB should be small enough for Euler method to work"
!!$       PRINT *, "expected value: .LE. 1.0d-6"
!!$       PRINT *, "current value:", dt_UB
!!$       res = 1
!!$    END IF
    IF (res .NE. SUCCESS) WRITE(error_message, *) "Error with data validation"
    RETURN 
  END FUNCTION check_ecc_data
  
  !!>****f* data_global/read_testdata
  !! SYNOPSIS  
  !!   Read testing data
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN
  !!   SUCCESS if ok
  !!   an error code if not
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  FUNCTION read_testdata(var_vect, filename) RESULT(res) 
    IMPLICIT NONE
    REAL(KIND=dp), DIMENSION(*) :: var_vect
    CHARACTER(FILENAME_LENGTH) :: filename
    INTEGER :: res
    !!local
    CHARACTER(40) :: var_name
    INTEGER :: ii
    CHARACTER :: dummy_char
    REAL(KIND=dp), DIMENSION(GROUP_LENGTH) :: temp_data
    INTEGER :: fileID, end_idx, error_code
    LOGICAL :: fexist
    REAL(KIND=dp) :: Vm, prcent_current_change, prcent_cap_change, prcent_orphan_sfu = 0d0
    REAL(KIND=dp) :: dummy
    
    res = SUCCESS

    INQUIRE(FILE=TRIM(filename), EXIST=fexist)
    IF (.NOT. fexist) THEN
       WRITE(error_message, *), "File " // TRIM(filename) // " not exist"
       res  = ERR_FILE_NOT_EXIST
       RETURN
    END IF

    IF (TRIM(where2save) .NE. '.') THEN
       sys_command = 'cp ' // TRIM(filename) // ' ' // TRIM(where2save)
       CALL system(sys_command)
    ENDIF
    
    fileID = get_file_ID()
    OPEN(unit=fileID, file=TRIM(filename), ACTION='READ')
    ii = 1
    error_code = 0
    DO WHILE (error_code == 0)
       READ(fileID, '(A1)', IOSTAT=error_code) dummy_char
       IF (dummy_char .NE. '#' .AND. dummy_char .NE. ' ') THEN
          BACKSPACE(fileID)
          READ(fileID,FormatIn, IOSTAT=error_code) var_vect(ii), var_name
          !WRITE(*,FormatOut) var_name, var_vect(ii)
          ii = ii + 1
          var_name = ""
       END IF
    END DO
    CLOSE(unit = fileID)
  END FUNCTION read_testdata
  
  !!>****f* data_global/cleanup
  !! SYNOPSIS  
  !!   Write out the last information
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN
  !!   
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  SUBROUTINE cleanup()
    IMPLICIT NONE
    
    !write out the last error message

    CALL write2log(error_message)
    PRINT *, error_message

    IF (TRIM(error_message) .EQ. "") THEN 
       IF (TRIM(where2save) .NE. ".") THEN
          !sys_command = "if [ `hostname` != leak ]; then scp -r " // TRIM(where2save) // " leak:/leak1/minhtuan; fi"
          !CALL system(sys_command)
          !PRINT *, "Data copied to leak"
       ENDIF
    ENDIF
    STOP
  END SUBROUTINE cleanup

END MODULE data_global

