C	This is a code written in FORTRAN that can be compiled on any machine with
C	an appropriate compiler. Compilation of the program requires IMSL mathematics 
C	library. We have tested the code using Digital fortran on PC.
C	The objective of the code is to calculate the cell conc. in the parralel-plate
C	flow chamber in regions near the substrate. The  number of cell-substrate collisions
C	is calcuated based on the trajectory.

C	Essential differences between the program code and the paper published in JIM are:
C	i) cgs units are used in this program code for selected parameters, like shear stress
C	ii) An additional parameter (theta_ar) is used here to account for the reversible release of 
C		adherent cells from the flow chamber substrate. This has been excluded from the 
C		manuscript text for simplicity reasons. This parameter can be set to 0 by the user
C		since this is typically the case in most experiments.
C	ii) Rolling Cell Conc (yc3) is expressed in per volume basis in the code, while in the 
C		paper it is expressed with units of cell density (i.e. per unit area basis). Essentially,
C		Cell density reported in the manuscript equals yc3*a 

C	Note: 
C	i)	This source code is accompanied by the file DATA.TXT, which contains some of the
C		model parameters values in text format. The executable (.exe) file should thus be placed
C		in the same	folder as the DATA.TXT file.
C	ii)	The output data from this program is saved in a file in text format. It can be 
C		imported into spread-sheet softwares, including Excel
C	iii) The default time step for calculation is 0.02 (this time is the dimensionless time 
C	    defined in the paper, t*=t*vset0/(2*b). This number is a good value for flow chamber runs
C		that last 5-10min. If you wish to only analyze data for short runs (like a few seconds),
C		 you must reduce the dimensionless time step to a smaller value to get accurate results 


	Program Main

	USE Numerical_Libraries
	external gcn1, gcn2, gcn3

	real*8 fycrit		! function to calculate the ycrit point for trajectory
	real*8 fNcoltot		! function to calculate the Ncoltot


C	geometry of the flow chamber
	real*8 b			! gasket half thickness 
	real*8 w			! gasket width 
	real*8 L			! gasket length 

C	flow properties
	real*8 visc			! viscosity 
	real*8 cd,md		! cell and media density 
	real*8 a			! cell radius 
	real*8 ml			! microville length
	real*8 vset0		! setting velocity 
	real*8 Q			! flow rate

	real*8 Cb			! inlet cell concentration 
	
	real*8 sstress		! shear stress in dyn/cm^2
	real*8 ur			! rolling velocity in m/s
	real*8 rw			! shear rate in /s
	
C	free stream velocity
	real*8 uf2			! cell free velocity in region 2

C	the maximum substrate occupation
	real*8 fmax			! the maximum percentage of substrate occupation area
						! in the flow chamber bottom		
	
C	make small grid in the flow chamber
	integer m			! m grid in x direction for the flow chamber
	parameter (m=11)
	real*8 x(m)			! the point along the flow chamber, dimensionless
	real*8 xx			! equal to each x(mm)
	real*8 delx			! x increment to calculate the trajectory

C	variables to calculate the Zcs and Zcc.	
	real*8 dely			! y increment to calculate the trajectory
	real*8 ycrit(m)		! corresponding point for the trajection to x(m)

C	primary tethering
	real*8 Ncoltot(m)	! total Ncol in cells/s below ycrit(m)
	real*8 Ncol(m)		! Ncol in cells/s for a specific region
						! between x(m-1) and x(m)

C	cell conc. (dimentionless form) for different region
	real*8 yc1(m), yc2(m), yc3(m), yprc1(m), yprc2(m), yprc3(m)
							! cell conc. for region 1, 2 and 3
	real*8 C1(m)			! cell conc. for region 1 =yc1(m)

C	time
	real*8 t1, t2, t3, delt		! time point for calculation, the time step
	real*8 tottime			! total running time for each run
	integer nstep			! total time steps

C	capture/release, transmigration flux
	real*8 Rp(m), R1(m), Ad(m), Ad1(m), TM(m)
							! tethering flux/reverse (cells/area/s)
C	secondary tethering
	real*8 Zcc1(m), Zcc2(m)	! cell-cell collision rate in R1 and R2
	real*8 Rs1(m), Rs2(m)	! secondary tethering flux in R1 and R2
	real*8 Rstot(m)			! total secondary tethering flux

C	total density
	real*8 CA(m)			! total adhesion cells at time t
	real*8 CTM(m)			! total transmigration cells at time t

C	all the frequency numbers
	real*8 Pfr, Pra, Par, Prf, Pcc, Pat		! frequency and possibility
	real*8 timeat			! average time used to transmigrate(s)

C	parameter for simulations where frequency values vary with iteration
	real*8 Pfr0, Pra0, Par0, Prf0, Pcc0 	! initial value
	real*8 delPfr, delPra, delPrf, delPar, delPcc	! interval
	integer istepPfr,istepPra,istepPcc,istepPrf,istepPar
						! the number of possible changes in frequency
	integer	i1, i2, i3, i4, i5
						! number of pfr, pra, pcc, prf, par counted in the code
	integer i,j,k,kk		! integer used for steps

C	for the output of the file
	integer optpara			! variable to control the output format
							! based on the variable which may be changed
	real*8  opt				! check how many variables have been changed
	real*8  option			! tell if the variable will be changed

C	save the data into file
	integer rn				! number of your run
	character out*2			! name of your run
	character outx1*15, outx2*15		! character for output axis
	real*8 output1(10,10)	! matrix used to store rolling density data
	real*8 output2(10,10)	! matrix used to store total adhesion data
	real*8 output3(10,10)	! matrix used to store total transmigration date
	real*8 x1(10), x2(10)	! array used to store the variables

	common/var1/a, ml
	common/var2/b,w,L
	common/var3/vset0
	common/var4/rw
	common/var5/Cb
	common/var6/Ncol
	common/var7/Rp, R1
	common/var8/Rs1
	common/var8a/Rs2
	common/var8b/Rstot
	common/var9/Ad, Ad1
	common/var10/ur
	common/var11/C1
	common/var12/outx2
	common/var13/x2
	common/var14/output1, output2, output3

CC	start the calculation	

C	define the run number
	write(*,*) 'If you make a MISTAKE when inputting parameters,'
	write(*,*) 'Please type Ctrl+C to abort the program'
	write(*,*)
	write(*,*)  '******START******'

C	get some default parameters for flow chamber and cells
C	default parameters are saved in file DATA.txt
C	they are able to change directly from that file
	
	open(unit=10, file='DATA.TXT', status='unknown')
	read(10,*) b		! gasket half thickness in inch
	read(10,*) w		! gasket width in cm
	read(10,*) L		! gasket length in m

	read(10,*) visc		! viscosity in cp
	read(10,*) cd		! cell density in g/ml
	read(10,*) md		! media density in g/ml
	read(10,*) a		! cell radius in micro-meter
	read(10,*) ml		! microvilli length in micro-meter
	read(10,*) fmax		! max. fraction of substrate area occupied by cells

	close(10)

	write(*,*) 'Please provide a run-number (any #less than 100)'
	read(*,*) rn

	if (rn.LT.10) then
	out='0'//char(rn+48)
	else
	out=char(INT(rn/10)+48)//
     &         char(rn-INT(10*INT(rn/10))+48)
	end if

	open(unit=11, file='RUN'//out//'.TXT', status='unknown')
	write(11, *) 'Run number   '//out//':', 'Simulation Result'
	write(*,*) 'Output from this run is stored in file: 
	1 RUN'//out//'.TXT'	
C	end of defining the output file	

C	Output some default parameters for flow chamber and cells

	write(11,*) 
	write(11,*) '*Default Parameters for Simulation'
	write(11,*)  'Flow Chamber Dimension'
	write(*,*)
	write(*,*)  'Flow Chamber Dimension'
	call outword('Half Flow Chamber Height b:',b,' inch')
	call outword('Flow Chamber Width w:', w, ' cm')
	call outword('Flow Chamber Length L:', L, ' cm')

	write(11,*)
	write(*,*)
	write(11,*) 'Cell and Medium properties'
	write(*,*) 'Cell and Medium properties'

	call outword('Medium Viscosity:', visc, ' cp')
	call outword('Cell Density:', cd,' g/ml')
	call outword('Medium Density:',md,' g/ml')
	call outword('Cell Radius a:', a, ' micro-meter')
	call outword('Cell Microvilli length:', ml, ' micro-meter')
	call outword('Fraction max. occupied area:', fmax, ' ')

C	change unit into SI
	b=b*2.54d-2			! unit: m
	w=w*0.01			! unit: m
	L=L*0.01			! unit: m
	cd=cd*1.0d3			! unit: kg/m^3
	md=md*1.0d3			! unit: kg/m^3
	a=a*1.0d-6			! unit: m
	ml=ml*1.0d-6		! unit: m
	
	write(*,*)
	write(*,*) 'If you wish to change the default parameters listed above'
	write(*,*) 'Please change them from file DATA.TXT'
	
	write(*,*)
	write(*,*) 'Please input the simulation variables requested below'
	write(11,*)
	write(11,*) '*Simulation Variables Input'
	
C	the shear stress
	write(*,*) 'Shear stress in dyn/cm^2?'
	read(*,*) sstress
	call outword('Shear stress:', sstress, ' dyn/cm^2')
		
C	the rolling velocity
	write(*,*) 'Rolling velocity in micro-meter/s?'
	read(*,*) ur
	call outword('Rolling veloicity ur:', ur,' micro-meter/s')
	ur=ur*1.0d-6		!m/s
	
C	shear rate
	rw=(sstress*0.1)/(visc*0.001)
	call outword('Shear rate:', rw, ' /s')

C	initial settling velocity
	vset0=(2*a*a*(cd-md)*9.81)/(9*visc*0.001)
C	free stream velocity
	uf2=rw*(a+ml/2)*(1-5/16*(a/(a+ml/2))**3)	!not dimensionlessed here
									!free stream velocity in R2	from Goldman
C	the inlet cell conc.
	write(*,*) 'Inlet cell conc. in million cells/ml?'
	write(*,*) 'Example, for 2x10^5 cells/ml, please input 0.2'
	read(*,*) Cb
	call outword('Inlet cell conc. Cb:', Cb, ' million cells/ml')
		
	Cb=Cb*1.0d6*1.0d6	! in /m^3
	write(*,*)

	write(*,*) 'To format the output file as discussed in the
	1manuscript you must choose two of the following parameters/
     2variables:'
	write(*,*) 't, x, Pfr, Pra, Pcc, Par, Prf'
	write(*,*)	
	write(*,*) 'These two parameters/variables can be varied over 
	1a resonable range. Once you have selected the parameters 
     2to vary, the program will ask you specific questions in order 
     3to assign the minimum and maximum values of the parameters. 
     4In order to limit the output file to a reasonable size, the 
     5results file will only provide cell rolling, adhesion and 
     6transmigration data for 10 different values for each variable 
     7at equal intervals.'
	write(*,*)
	call outword('Inlet cell conc. Cb:', Cb*1.0d-12, ' million cells/ml')
	write(*,*) 'The other five parameters/variables must be assigned 
	1fixed values, i.e. they cannot be varied during the simulation'
	write(*,*)
	write(*,*) 'See manuscript for details or just try a run!'
	
C	the variable which control the output format with different varialbes
	optpara=0
	opt=0

C	set the following optvalue as output control value for 7 variable
C	for t: opt=50
C	for x: opt=20
C	for theta_fr: opt=1
C	for theta-ra: opt=2
C	for theta_cc: opt=3
C	for theta_ar: opt=4
C	for theta_rf: opt=5

C	for any two variables, the sum of the optpara
C	<10, --->p and p change
C	between (20, 30) ---> x, and any one of p
C	between (50, 60) ---> t, and any one of p
C	>60, --->t and x change

C	start to choose the variables
C	time, t
	write(*,*)
	write(*,*) 'Is t one of the variables you wish to change?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option
	write(*,*)

CC	set the time for each run	
	write(*,*) 'How long is the flow chamber run, in second?'
	read(*,*) tottime
	call outword('Total time:', tottime, ' sec.')
	write(11,*)
	write(*,*)

	delt=0.02	! make time step to be 0.02, dimensionless
	nstep=int(tottime*1.0/(2*b/vset0*delt))+1

C	set optpara for later output use
	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have decided to output data as a function of
	1 variable t'
	write(*,*)
	option=50
	optpara=optpara+option
	else
	opt=opt+0
	write(*,*) 'You have decided not to obtain output as a function of
	1 t'
	write(*,*) 'Only the result at the end of the flow chamber run 
	1 will be output'
	write(*,*)
	optpara=optpara+0
	end if
	
CC	distance along the flow chamber, dimensionless x
	write(*,*) 'Is x one of the variables you wish to change?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option
	write(*,*)

C	set optpara for later output use
	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have decided to output data as a function of
	1 variable x'
	write(*,*)
	option=20
	optpara=optpara+option
	else
	opt=opt+0
	write(*,*) 'You have decided not to obtain output as a function of
	1 x'
	write(*,*) 'Only the result at the center of the flow chamber run 
	1 will be output (i.e. x*=0.5)'
	write(*,*)
	optpara=optpara+0
	end if
	
CC	set the frequency numbers

CCC	theta_fr
C	check if you have changed two variables
	if (opt.ge.2) then
	write(*,*) 'You have already chosen two variables'
	write(*,*) 'You now need to fix all the frequency numbers'
	write(*,*)
CC	theta_fr
	write(*,*) 'Please provide Theta_fr value (/m)'
	read(*,*) Pfr0
	call outword('Theta_fr=', Pfr0, ' /m')
	delPfr=0
	istepPfr=1
C	theta_ra
	write(*,*) 'Please provide Theta_ra value (/m)'
	read(*,*) Pra0
	call outword('Theta_ra=', Pra0,' /m')
	delPra=0
	istepPra=1
C	theta_cc
	write(*,*) 'Please provide Theta_cc value'
	read(*,*) Pcc0
	call outword('Theta_cc=', Pcc0, ' ')
	delPcc=0
	istepPcc=1
C	theta_ar
	write(*,*) 'Please provide Theta_ar value (/s)'
	read(*,*) Par0
	call outword('Theta_ar=', Par0, ' /s')
	delPar=0
	istepPar=1
C	theta_rf	
	write(*,*) 'Please provide Theta_rf value (/m)'
	read(*,*) Prf0
	call outword('Theta_rf=', Prf0, ' /m')
	delPrf=0
	istepPrf=1

C	*1: else
	else

	if (opt.eq.0) write(*,*) 'You still have to choose two parameters/
	1variables to change'
	if (opt.eq.1) write(*,*) 'You still have to choose one parameter/
	1variable to change'
	write(*,*)
	write(*,*) 'Do you wish to choose Theta_fr as an output variable?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option
	write(*,*)

	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have chosen Theta_fr as a variable'
	write(*,*) 'You need to set the interval over which Theta_fr
	1 will be varied. For this, please choose the initial/lowest Theta_
     2fr value and the increment below. Remember that the maximum Theta_
	3fr value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_fr value (/m)'
	read(*,*) Pfr0
	write(*,*) 'Increment for Theta_fr value (/m)'
	read(*,*) delPfr
	istepPfr=10
	call outword('Theta_fr0=', Pfr0, ' /m') 
	call outword('increment Theta_fr=', delPfr,' /m')
	call outword1('Total number of Theta_fr values=', istepPfr)
	option=1
	optpara=optpara+option

	else
	opt=opt+0
	write(*,*) 'You have decided not to choose Theta_fr as a variable'
	write(*,*)
	optpara=optpara+0
C	set the theta_fr value
	write(*,*) 'Please provide Theta_fr value (/m)'
	read(*,*) Pfr0
	call outword('Theta_fr=', Pfr0, ' /m')
	delPfr=0
	istepPfr=1
	end if

CCC	theta_ra
	if (opt.ge.2) then
	write(*,*) 'You have already chosen two variables'
	write(*,*) 'You now need to fix all remaining frequency numbers'
	write(*,*)
	
C	theta_ra
	write(*,*) 'Please provide Theta_ra value (/m)'
	read(*,*) Pra0
	call outword('Theta_ra=', Pra0, ' /m')
	delPra=0
	istepPra=1
C	theta_cc
	write(*,*) 'Please provide Theta_cc value'
	read(*,*) pcc0
	call outword('Theta_cc=', Pcc0, ' ')
	delPcc=0
	istepPcc=1
C	theta_ar
	write(*,*) 'Please provide Theta_ar value (/s)'
	read(*,*) Par0
	call outword('Theta_ar=', Par0,' /s')
	delPar=0
	istepPar=1
C	theta_rf	
	write(*,*) 'Please provide Theta_rf value (/m)?'
	read(*,*) Prf0
	call outword('Theta_rf=', Prf0, ' /m')
	delPrf=0
	istepPrf=1

C	*2: else
	else

	if (opt.eq.0) write(*,*) 'You still have to choose two parameters/
	1variables to change'
	if (opt.eq.1) write(*,*) 'You still have to choose one parameter/
	1variable to change'
	write(*,*)
	
	write(*,*) 'Do you wish to choose Theta_ra as an output variable?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option
	write(*,*)

	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have chosen Theta_ra as a variable'
	write(*,*) 'You need to set the interval over which Theta_ra
	1 will be varied. For this, please choose the initial/lowest Theta_
     2ra value and the increment below. Remember that the maximum Theta_
	3ra value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_ra value (/m)'
	read(*,*) Pra0
	write(*,*) 'Increment for Theta_ra value (/m)'
	read(*,*) delPra
	istepPra=10
	call outword('Theta_ra0=', Pra0, ' /m') 
	call outword('increment Theta_ra=', delPra,' /m')
	call outword1('Total number of Theta_ra values=', istepPra)
	option=2
	optpara=optpara+option

	else
	opt=opt+0
	write(*,*) 'You have decided not to choose Theta_ra as a variable'
	write(*,*)
	optpara=optpara+0
C	set the theta_ra value
	write(*,*) 'Please provide Theta_ra value (/m)'
	read(*,*) Pra0
	call outword('Theta_ra=', Pra0, ' /m')
	delPra=0
	istepPra=1
	end if


CCC	theta_cc
	write(*,*)

	if (opt.ge.2) then
	write(*,*) 'You have already chosen two variables'
	write(*,*) 'You now need to fix all remaining frequency numbers'
	write(*,*)

C	theta_cc
	write(*,*) 'Please provide Theta_cc value'
	read(*,*) pcc0
	call outword('Theta_cc=', Pcc0, ' ')
	delPcc=0
	istepPcc=1
C	theta_ar
	write(*,*) 'Please provide Theta_ar value (/s)'
	read(*,*) Par0
	call outword('Theta_ar=', Par0,' /s')
	delPar=0
	istepPar=1
C	theta_rf	
	write(*,*) 'Please provide Theta_rf value (/m)'
	read(*,*) Prf0
	call outword('Theta_rf=', Prf0, ' /m')
	delPrf=0
	istepPrf=1


C	*3: else
	else

	if (opt.eq.0) write(*,*) 'You still have to choose two parameters/
	1variables to change'
	if (opt.eq.1) write(*,*) 'You still have to choose one parameter/
	1variable to change'
	write(*,*)
	
	write(*,*) 'Do you wish to choose Theta_cc as an output variable?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option
	write(*,*)

	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have chosen Theta_cc as a variable'
	write(*,*) 'You need to set the interval over which Theta_cc
	1 will be varied. For this, please choose the initial/lowest Theta_
     2cc value and the increment below. Remember that the maximum Theta_
	3cc value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_cc value'
	read(*,*) Pcc0
	write(*,*) 'Increment for Theta_cc value'
	read(*,*) delPcc
	istepPcc=10
	call outword('Theta_cc0=', Pcc0, ' ') 
	call outword('Increment Theta_cc=', delPcc,' ')
	call outword1('Total number of Theta_cc values=', istepPcc)
	option=3
	optpara=optpara+option
	
	else
	opt=opt+0
	write(*,*) 'You have decided not to choose Theta_cc as a variable'
	write(*,*)
	optpara=optpara+0
C	set the theta_cc value
	write(*,*) 'Please provide Theta_cc value'
	read(*,*) Pcc0
	call outword('Theta_cc=', Pcc0, '  ')
	delPcc=0
	istepPcc=1
	end if
	
CCC	theta_ar

	if (opt.ge.2) then
	write(*,*) 'You have already chosen two variables'
	write(*,*) 'You now need to fix all remaining frequency numbers'
	write(*,*)

C	theta_ar
	write(*,*) 'Please provide Theta_ar value (/s)'
	read(*,*) Par0
	call outword('Theta_ar=', Par0,' /s')
	delPar=0
	istepPar=1
C	theta_rf	
	write(*,*) 'Please provide Theta_rf value (/m)'
	read(*,*) Prf0
	call outword('Theta_rf=', Prf0, ' /m')
	delPrf=0
	istepPrf=1

C	*4: else
	else

	if (opt.eq.0) then
	write(*,*) 'You need to change the last two frequencies'
	write(*,*)
C	theta_ar
	write(*,*) 'You need to set the interval over which Theta_ar
	1 will be varied. For this, please choose the initial/lowest Theta_
     2ar value and the increment below. Remember that the maximum Theta_
	3ar value will be =initial value+9*increment'
	write(*,*) 'Set the initial theta_ar value (/s)'
	read(*,*) Par0
	write(*,*) 'Interval between theta_ar value (/s)'
	read(*,*) delPar
	write(*,*)
	istepPar=10
	call outword('Theta_ar0=', Par0, ' /s') 
	call outword('Increment Theta_ar=', delPar,'  /s')
	call outword1('Total number of Theta_ar values=', istepPar)
	option=4
	optpara=optpara+option

C	theta_rf	
	write(*,*) 'You need to set the interval over which Theta_rf
	1 will be varied. For this, please choose the initial/lowest Theta_
     2rf value and the increment below. Remember that the maximum Theta_
	3rf value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_rf value, (/m)'
	read(*,*) Prf0
	write(*,*) 'Increment for Theta_rf value, (/m)'
	read(*,*) delPrf
	istepPrf=10
	call outword('Theta_rf0=', Prf0, ' /m ') 
	call outword('Increment Theta_rf=', delPrf,' /m')
	call outword1('Total number of Theta_rf values=', istepPrf)
	option=5
	optpara=optpara+option
	
	end if
	
	if (opt.eq.1) then
	if (opt.eq.1) write(*,*) 'You still have to choose one parameter/
	1variable to change'
	write(*,*)
	write(*,*) 'Do you wish to choose Theta_ar as an output variable?'
	write(*,*) 'If YES, input 1'
	write(*,*) 'If NO, input 0'
	read(*,*) option
	write(*,*)

	if (option.eq.1) then 
	opt=opt+1		! one variable has been changed here!
	write(*,*) 'You have chosen Theta_ar as a variable'
	write(*,*) 'You need to set the interval over which Theta_ar
	1 will be varied. For this, please choose the initial/lowest Theta_
     2ar value and the increment below. Remember that the maximum Theta_
	3ar value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_ar value (/s)'
	read(*,*) Par0
	write(*,*) 'Inrement for Theta_ar value (/s)'
	read(*,*) delPar
	istepPar=10
	call outword('Theta_ar0=', Par0, ' /s ') 
	call outword('increment Theta_ar=', delPar,' /s')
	call outword1('Total number of Theta_ar values=', istepPar)
	option=4
	optpara=optpara+option
	
C	theta_rf	
	write(*,*) 'Please provide Theta_rf value (/m)'
	read(*,*) Prf0
	call outword('Theta_rf=', Prf0, ' /m')
	delPrf=0
	istepPrf=1
	
	else

	write(*,*) 'Please provide Theta_ar value (/s)'
	read(*,*) Par0
	call outword('Theta_ar=', Par0, ' /s')
	delPar=0
	istepPar=1

	write(*,*) 'You have chosen Theta_rf as a variable'
	write(*,*) 'You need to set the interval over which Theta_rf
	1 will be varied. For this, please choose the initial/lowest Theta_
     2rf value and the increment below. Remember that the maximum Theta_
	3rf value will be =initial value+9*increment'
 	write(*,*) 'Set the initial Theta_rf value (/m)'
	read(*,*) Prf0
	write(*,*) 'Increment for Theta_rf value (/m)'
	read(*,*) delPrf
	istepPrf=10
	call outword('Theta_rf0=', Prf0, ' /m') 
	call outword('Increment Theta_rf=', delPrf,' /m')
	call outword1('Total number of Theta_rf values=', istepPrf)
	option=5
	optpara=optpara+option

	end if
	end if		
	end if		!	*4
	end if		!	*3
	end if		!	*2
	end if		!	*1

CCCCCC

C	set parameter for the transmigration
	write(*,*) 'Does the system have transmigration?'
	write(*,*) 'if YES, input 1'
	write(*,*) 'if NO, input 0'
	read(*,*) option				
						! final option: for the transmigration
	if (option.eq.1) then
	write(*,*) 'Input the time taken for adherent cells to 
	1transmigrate, (s)'
	read(*,*) timeat
	call outword('Time for transmigration:', timeat,'s')
	Pat=0.693/timeat
	else
	Pat=0				! 1/Pat: infinite time to transmigrate
	write(11,*) 'No transmigration'
	write(*,*) 'No transmigration'
	end if



CC	calculate the trajectory and flux
C	define the dely, the dimensionless form

	dely=1.0d-3
	dely=dely*a

C	define the x along the flow chamber
	delx=L/(m-1)

C	calculate each trajectory from (x(k),a+ml) to (0, ycrit(k))

	do k=1,m
		x(k)=delx*(k-1)
		ycrit(k)=fycrit(x(k), ml+3*a, dely, 0.001*a)
		Ncoltot(k)=fNcoltot(ycrit(k))
	end do

C	to calculate the collision in each part of x
	do k=1,m
		if (k.eq.1) then
		Ncol(k)=0
		else
		Ncol(k)=Ncoltot(k)-Ncoltot(k-1)
		end if
	end do


	i=0							! i: total number of calculation

	i1=0
9	continue
	i1=i1+1
	Pfr=Pfr0+(i1-1)*delPfr

	i2=0
99	continue
	i2=i2+1
	Pra=Pra0+(i2-1)*delPra

	i3=0
999	continue
	i3=i3+1
	Pcc=Pcc0+(i3-1)*delPcc

	i4=0
9999	continue
	i4=i4+1
	Par=Par0+(i4-1)*delPar

	i5=0
99999 continue
	i5=i5+1
	Prf=Prf0+(i5-1)*delPrf


C	setup the initial condition
C	calculate R1, R2 and R3 at the same time
C	time is dimensionless
	
	t1=0.0
	t2=t1
	t3=t2
	
	do k=1, m
		yc1(k)=0.0
		yc2(k)=0.0
		yc3(k)=0.0
		yprc1(k)=0.0
		yprc2(k)=0.0
		yprc3(k)=0.0
	end do

C	define Rp, R1, Ad, CA, Ad1, Rs
	do k=1, m
		Rp(k)=0
		Rs1(k)=0
		Rs2(k)=0
		Rstot(k)=0
		R1(k)=0
		Ad(k)=0
		CA(k)=0
		CTM(k)=0
		Ad1(k)=0
	end do

	istep=1

11111	continue

C	region 1
	ido=1
	call ddaspg (m,t1,t1+delt,ido,yc1,yprc1,gcn1)
	ido=3
	call ddaspg (m,t1,t1+delt,ido,yc1,yprc1,gcn1)

	do k=1, m
		C1(k)=yc1(k)
	end do

C	region 2
	ido=1
	call ddaspg (m,t2,t2+delt,ido,yc2,yprc2,gcn2)
	ido=3
	call ddaspg (m,t2,t2+delt,ido,yc2,yprc2,gcn2)

C     region 3
	ido=1
	call ddaspg (m,t3,t3+delt,ido,yc3,yprc3,gcn3)
	ido=3
      call ddaspg (m,t3,t3+delt,ido,yc3,yprc3,gcn3)

	do k=1,m
	
C	calculate the Rp, R1, Ad, Ad1, CA, dimensionless form
	Rp(k)=ml*uf2*Pfr*yc2(k)/vset0	! from flux definition
	Rp(k)=Rp(k)*(1-(yc3(k)*a+CA(k)*vset0)*Cb
     &      *3.14*a**2/fmax)			! adjust because of limit of surface
	R1(k)=a*ur*Prf*yc3(k)/vset0
	Ad(k)=a*ur*Pra*yc3(k)/vset0
	Ad1(k)=CA(k)*Par
	TM(k)=CA(k)*Pat
	CA(k)=CA(k)+(Ad(k)-Ad1(k)-TM(k))*delt/vset0*2*b
	CTM(k)=CTM(k)+TM(k)*delt/vset0*2*b
	Zcc1(k)=(C1(k)*Cb)*(4*a*a*rw*(3.14*(a+ml)/2+4/3*a))*0.947765
     &	*(yc3(k)*a+CA(k)*vset0)*Cb
     			! 0.947765: prefactor to adjust based on goldman's eq.
	Rs1(k)=Pcc*Zcc1(k)/Cb/vset0
	Rs1(k)=Rs1(k)*(1-(yc3(k)*a+CA(k)*vset0)*Cb
     &      *3.14*a**2/fmax)			! adjust because of limit of surface
	
	Zcc2(k)=uf2*(yc2(k)*Cb)
     &       *(ml*4*a)*(yc3(k)*a+CA(k)*vset0)*Cb
	Rs2(k)=Pcc*Zcc2(k)/Cb/vset0
	Rs2(k)=Rs2(k)*(1-(yc3(k)*a+CA(k)*vset0)*Cb
     &      *3.14*a**2/fmax)			! adjust because of limit of surface

      Rstot(k)=Rs1(k)+Rs2(k)

	end do

CCCCC start to store the data

CC
CC	one variable is t
	if (optpara.gt.50) then			!t
C	only save ten time points in the whole time peroid
	outx1='time(s)'		! mark the axis for the time  

	if (mod(istep,int(nstep/10)).eq.0.and.
     &	(istep/int(nstep/10)).le.10) then

	kk=int(istep/int(nstep/10))
	x1(kk)=istep*2*b/vset0*delt
	if (optpara.gt.60) write(*,*) kk*10, '% completed......'

C	set x2 (either x or frequency)
	if(mod(optpara,10).eq.0) then	!t, x*(0.1, 1)
	do j=1,m-1
	xx=0.1*j
	call output('x*',j, xx, kk, j, yc3(j+1), CA(j+1), CTM(j+1)) 
	end do
	end if
C	any frequency as varialbe
	if (mod(optpara,10).eq.1)		!t, Pfr
     &    call output('theta_fr(/m)',i1,Pfr,kk,i1,yc3(6),CA(6),CTM(6))
	if (mod(optpara,10).eq.2)		!t, Pra
     &	call output('theta_ra(/m)',i2,Pra,kk,i2,yc3(6),CA(6),CTM(6))
	if (mod(optpara,10).eq.3)		!t, Pcc
     &	call output('theta_cc',i3,Pcc,kk,i3, yc3(6), CA(6), CTM(6))
 	if (mod(optpara,10).eq.4)		!t, Par
     &	call output('theta_ar(/s)',i4,Par,kk,i4,yc3(6),CA(6),CTM(6))
	if (mod(optpara,10).eq.5)		!t, Prf
     &	call output('theta_rf(/m)',i5,Prf,kk,i5,yc3(6),CA(6),CTM(6))

	end if
	end if

	istep=istep+1

	if (istep.le.nstep) then			    !A
	go to 11111
	else									!A

CC	x1: x;  x2: frequency
	if (optpara.gt.20.and.optpara.lt.50) then		!x*
	outx1='x*'			! mark axis as x
	do j=1,m-1
	x1(j)=0.1*j
C	any frequency as another variable
	if (mod(optpara,10).eq.1)		!x*, Pfr
     &	call output('theta_fr(/m)',i1,Pfr,j,i1,yc3(j+1),CA(j+1),CTM(j+1))
	if (mod(optpara,10).eq.2)		!x*, Pra
     &	call output('theta_ra(/m)',i2,Pra,j,i2,yc3(j+1),CA(j+1),CTM(j+1))
	if (mod(optpara,10).eq.3)		!x*, Pcc
     &	call output('theta_cc',i3,Pcc,j, i3, yc3(j+1), CA(j+1), CTM(j+1))
	if (mod(optpara,10).eq.4)		!x*, Par
     &	call output('theta_ar(/s)',i4,Par,j,i4,yc3(j+1),CA(j+1),CTM(j+1))
	if (mod(optpara,10).eq.5)		!x*, Prf
     &	call output('theta_rf(/m)',i5,Prf,j,i5,yc3(j+1),CA(j+1),CTM(j+1))
	end do
	end if


CC	only frequencies as variables
	if (optpara.lt.10) then			!B

C	theta_fr
	if (istepPfr.gt.1) then			!C
	x1(i1)=Pfr
	outx1='theta_fr(/m)'					! mark axis for x1
	if (istepPra.gt.1)				!Pfr, Pra
     &	call output('theta_ra(/m)',i2,Pra,i1,i2,yc3(6),CA(6),CTM(6))
	if (istepPcc.gt.1)				!Pfr, Pcc
     &	call output('theta_cc',i3,Pcc,i1,i3,yc3(6),CA(6),CTM(6))
	if (istepPar.gt.1)				!Pfr, Par
     &	call output('theta_ar(/s)',i4,Par,i1,i4,yc3(6),CA(6),CTM(6))
	if (istepPrf.gt.1)				!Pfr, Prf	
     &	call output('theta_rf(/m)',i5,Prf,i1,i5,yc3(6),CA(6),CTM(6))
	
	else							!C

C	pra
	if (istepPra.gt.1) then			!D
	x1(i2)=Pra
	outx1='theta_ra(/m)'					!mark axis for x1
	if (istepPcc.gt.1)				!Pra, Pcc
     &	call output('theta_cc',i3, Pcc, i2, i3, yc3(6), CA(6), CTM(6))
	if (istepPar.gt.1)				!Pra, Par
     &	call output('theta_ar(/s)',i4,Par,i2,i4,yc3(6), CA(6), CTM(6))
	if (istepPrf.gt.1)				!Pra, Prf
     &	call output('theta_rf(/m)',i5,Prf,i2,i5,yc3(6), CA(6), CTM(6))

	else							!D

C	pcc
	if (istepPcc.gt.1) then			!E	
	x1(i3)=Pcc
	outx1='theta_cc'						!mark axis as x1
	if (istepPar.gt.1)				!Pcc, Par
     &	call output('theta_ar(/s)',i4,Par,i3,i4,yc3(6),CA(6),CTM(6))
	if (istepPrf.gt.1)				!Pcc, Prf
     &	call output('theta_rf(/m)',i5,Prf,i3,i5,yc3(6),CA(6),CTM(6))

C	par and prf
	else
	x1(i4)=Par						!Par, Prf
	outx1='theta_ar(/s)'					!mark axis for x1
	call output('theta_rf(/m)',i5, Prf, i4, i5, yc3(6), CA(6), CTM(6))
	
	end if							!E
	end if							!D
	end if							!C
	end if							!B
	end if							!A
	
	i=i+1
	if (optpara.gt.10.and. optpara.lt.60) then
	write(*,*) i*10, '% completed......'
	else
	if (mod(i,10).eq.0) write(*,*) i, '% completed......'
	end if

	if (i5.lt.istepprf) go to 99999
	if (i4.lt.isteppar) go to 9999
	if (i3.lt.isteppcc) go to 999
	if (i2.lt.isteppra) go to 99
	if (i1.lt.isteppfr) go to 9

	write(*,*) 'DONE'


C	end of calculation

C	print into file

C	for the value of x2

	write(11,*)
	write(11,*) '*Simulation Results:'
	write(11,*)


	write(11,*) 'Table I: Rolling Cell Density (cells/cm^2)'
	write(11,*) '-------------------------------------------->', outx2
	write(11,887) '|', 'p1\p2', (x2(k), k=1,10)
	do j=1, 8
	write(11,888) '|', x1(j), (output1(j,k), k=1,10)
	end do
	write(11,888) 'V', x1(9), (output1(9,k), k=1,10)
	write(11,888) outx1, x1(10), (output1(10,k), k=1,10)
	write(11,*)

	write(11,*) 'Table II: Total Adherent Cell Density (cells/cm^2)'
	write(11,*) '-------------------------------------------->', outx2
	write(11,887) '|', 'p1\p2', (x2(k), k=1,10)
	do j=1, 8
	write(11,888) '|', x1(j), (output2(j,k), k=1,10)
	end do
	write(11,888) 'V', x1(9), (output2(9,k), k=1,10)
	write(11,888) outx1, x1(10), (output2(10,k), k=1,10)
	write(11,*)
887	format(1x, 2A15, 1f15.2, 9f15.3)
888	format(1x, A15, 1f15.2, 10f15.3)

C	For the transmigration
	if (option.eq.1) then	
	write(11,*) 'Table III: Total Transmigration Cell Density (cells/cm^2)'
	write(11,*) '-------------------------------------------->', outx2
	write(11,887) '|', 'p1\p2', (x2(k), k=1,10)
	do j=1, 8
	write(11,888) '|', x1(j), (output3(j,k), k=1,10)
	end do
	write(11,888) 'V', x1(9), (output3(9,k), k=1,10)
	write(11,888) outx1, x1(10), (output3(10,k), k=1,10)
	write(11,*)
	else
	write(11,*) 'No transmigration cells since it takes infinite time'
	end if

	write(11,*) 
	write(11,*) '*End'

	close(11)

	end

CC 


CCC	function to calculate the ycrit

	FUNCTION fycrit(xstart, ystart, dely, xend)
C	xstart, ystart: the starting x and y point to get the trajectory
C	dely: y step

	real*8 fycrit, xstart, ystart, dely, xend
	real*8 x1, y1, x2, y2
	
	real*8 a, ml, b, w, L
	real*8 Q, vset0
	real*8 rw
	common/var1/a,ml
	common/var2/b,w,L
	common/var3/vset0
	common/var4/rw

	Q=2*b*b*w*rw/3

	if (xstart.eq.0) then
	fycrit=0
	go to 7
	end if

	x1=xstart
	y1=ystart

77	y2=y1+dely	

C	calculate the x2, the second point
	if (y2.ge.4*a+ml) then

	if (y2. lt.2*b) then
	x2=x1-Q*(-(y2**3-y1**3)-3/2*a*(y2**2-y1**2)+3*b*(y2**2-y1**2)
     &	+(-3*a**2+6*a*b)*(y2-y1)
     &	+(-3*a**3+6*a**2*b)*LOG((y2-a)/(y1-a)))/(4*vset0*b**3*w)
	else
	fycrit=y1
	end if

	else

	x2=x1-3*Q*(a*(y2-y1)+0.5*(y2**2-y1**2)+5/16*a*a
     &	*LOG(y2/y1)+11/16*a*a*LOG((a-y2)/(a-y1)))
     &       /(2*b**2*w*vset0)
	
	end if
	
C	check if x2 will be 0, if so, y2=yend	
	if (x2.lt.xend) then
		fycrit=y2
	else
		x1=x2
		y1=y2
		go to 77
	end if

7	end

CC	end of function fycrit

CC
CCC	function of fNcoltot to calculate the total collision number

	FUNCTION fNcoltot(yend)

	real*8 yend
	real*8 fNcoltot

	
	real*8 a, ml, b, w, L
	real*8 Cb, rw

	common/var1/a, ml
	common/var2/b,w,L
	common/var4/rw
	common/var5/Cb

	if (yend.le.(a+ml)) then
	fNcoltot=0
	else

	if (yend.ge.(4.0*a+ml)) then
	fNcoltot=Cb*w*rw*((yend**2)/2-yend**3/6/b
     &     -(4.0*a+ml)**2/2+(4.0*a+ml)**3/6/b)
C	from 4cr+ml to cr+ml
	fNcoltot=fNcoltot+Cb*w*rw*((4.0*a+ml)**2/2
     &     +5/16*a**3/(4.0*a+ml)-0.5*(a+ml)**2-5/16*a**3/(a+ml))     		 
	else 
	fNcoltot=Cb*w*rw*(0.5*yend**2+5/16*a**3*(1/yend)
     &     -0.5*(a+ml)**2-5/16*a**3/(a+ml))     		 
	end if

	end if

	end

C	end of calculation of Ncoltot


CCC	subroution for region 1,2 and 3


	subroutine gcn1 (m,t1,yc1,yprc1,gval1)

	integer m, mm
	parameter(mm=11)
	real*8 t1, yc1(m), yprc1(m), gval1(m)
	real*8 rhs1(m)

	real*8 delx
	real*8 dcdx(m)
	real*8 para1, para2
	real*8 uf, vset

	integer i

	real*8 a,ml,b,w,L
	real*8 rw
	real*8 Cb, vset0
	real*8 Ncol(mm)
	real*8 Rs1(mm)

	common/var1/a, ml
	common/var2/b,w,L
	common/var3/vset0
	common/var4/rw
	common/var5/Cb
	common/var6/Ncol
	common/var8/Rs1			


C     Defining the equation in the domain
C	Defining the definition for dc/dx, three-point formulus
C	dimensionless here
	
	uf=rw*(2*a+ml)*(1-5/16*(a/(2*a+ml))**3)
	vset=1/(1+a/(a+ml-a))

	para1=-2*uf*b/(L*vset0)
	para2=-2*b/(2*a)

	delx=1.0/(m-1)

	dcdx(2)=(yc1(3)-yc1(1))/2/delx

	do i=3,m
	dcdx(i)=(yc1(i-2)-4*yc1(i-1)+3*yc1(i))/2/delx
	end do	

	do i=2,m
		rhs1(i)=para1*dcdx(i)
		rhs1(i)=rhs1(i)+para2*(yc1(i)*vset+Rs1(i))
		rhs1(i)=rhs1(i)+2*b*(m-1)*Ncol(i)/(L*(2*a)*w*vset0*Cb)
	end do
  

	do i=1, m
	  if (i.eq.1) then
	     gval1(i)=1.0-yc1(i)
	  else
	     gval1(i)=rhs1(i)-yprc1(i)
	  end if
	end do
 
	end 	


	subroutine gcn2 (m,t2,yc2,yprc2,gval2)

	integer m, mm
	parameter(mm=11)
	real*8 t2, yc2(m), yprc2(m), gval2(m)
	real*8 rhs2(m)

	real*8 delx
	real*8 dcdx(m)
	real*8 Rp(mm), R1(mm), Rs2(mm)
	real*8 C1(mm)

	integer i

	real*8 uf, vset
	real*8 a,ml,b,w,L
	real*8 rw
	real*8 Cb, vset0

	real*8 para1, para2, para3

	common/var1/a, ml
	common/var2/b,w,L
	common/var3/vset0
	common/var4/rw
	common/var5/Cb
	common/var7/Rp, R1
	common/var8a/Rs2
	common/var11/C1

C     Defining the equation in the domain
C	Defining the definition for dc/dx, three-point formulus
C	dimensionless here
	
	uf=rw*(a+ml/2)*(1-5/16*(a/(a+ml/2))**3)
	vset=1/(1+a/(a+ml-a))

	para1=-2*uf*b/(L*vset0)
	para2=2*b/ml
	para3=-2*b/ml

	delx=1.0/(m-1)

	dcdx(2)=(yc2(3)-yc2(1))/2/delx

	do i=3,m           
	dcdx(i)=(yc2(i-2)-4*yc2(i-1)+3*yc2(i))/2/delx
	end do	

	do i=2,m
		rhs2(i)=para1*dcdx(i)
		rhs2(i)=rhs2(i)+para2*vset*C1(i)
		rhs2(i)=rhs2(i)+para3*(Rp(i)+Rs2(i)-R1(i))
	end do
  

	do i=1, m
	  if (i.eq.1) then
	     gval2(i)=1.0-yc2(i)
	  else
	     gval2(i)=rhs2(i)-yprc2(i)
	  end if
	end do
 
	end 

C     for region 3

	subroutine gcn3 (m,t3,yc3,yprc3,gval3)

	integer m, mm
	parameter(mm=11)
	real*8 t3, yc3(m), yprc3(m), gval3(m)
	real*8 rhs3(m)

	integer i
	
	real*8 delx
	real*8 dcdx(m)
	real*8 Rp(mm), R1(mm)
	real*8 Ad(mm), Ad1(mm)
	real*8 Rstot(mm)
	
	real*8 a, ml, b, w, L
	real*8 ur
	real*8 vset0
	real*8 para1, para2

	common/var1/a, ml
	common/var2/b,w,L
	common/var3/vset0
	common/var7/Rp, R1
	common/var8b/Rstot
	common/var9/Ad, Ad1
	common/var10/ur

C     Defining the equations for the domain
C	Defining the definition for dc/dx, three-point formulus
	
	para1=-ur*2*b/(vset0*L)
	para2=2*b/a

	delx=1.0/(m-1)
	dcdx(2)=(yc3(3)-yc3(1))/2/delx

	do i=3,m
	dcdx(i)=(yc3(i-2)-4*yc3(i-1)+3*yc3(i))/2/delx
	end do

	do i=2,m
	    rhs3(i)=para1*dcdx(i)
     &		+para2*(Rp(i)+Rstot(i)-R1(i)-Ad(i)+Ad1(i))
      end do
	   

	do i=1, m

	   if (i.eq.1) then
	      gval3(i)=yc3(i)-0.0
	   else
	      gval3(i)=-rhs3(i)+yprc3(i)
	   end if

	end do

	end 

C	subroutine for output 
	subroutine outword (word1, val, word2)

	character *(*)word1, word2	! character for output words
	real*8 val
	character word11*40, word22*20

	word11=word1
	word22=word2
	write(11,889) word11, val, word22
	write(*,889) word11, val, word22
889	format(1x, A40, f8.3, A20)

	end

	subroutine outword1 (word, val)

	character *(*)word	! character for output words
	integer val
	character word11*40

	word11=word

	write(11,8899) word11, val
	write(*,8899) word11, val
8899	format(1x, A40, I6)

	end

	subroutine output(outx2c, i, x2i, j, k, yc3, CA, CTM)

	character *(*)outx2c, outx2cc*15
	integer i, j, k
	real*8 x2i
	real*8 yc3, CA, CTM
	
	real*8 a, ml, Cb, vset0
	character outx2*15

	real*8 x2(10)
	real*8 output1(10,10)	! matrix used to store rolling density data
	real*8 output2(10,10)	! matrix used to store total adhesion data
	real*8 output3(10,10)	! matrix used to store total transmigration date

	common/var1/a, ml
	common/var3/vset0
	common/var5/Cb
	common/var12/outx2
	common/var13/x2
	common/var14/output1, output2, output3

	outx2cc=outx2c
	outx2=outx2cc
	x2(i)=x2i
	output1(j,k)=yc3*Cb*a*1.0d-4
	output2(j,k)=CA*Cb*vset0*1.0d-4
	output3(j,k)=CTM*Cb*vset0*1.0d-4

	end