Subroutine readin
Use comm
Implicit None
Real(Kind=8), Parameter :: cexpo=5d-3, gexp=1d0+cexpo
Real(Kind=8), Parameter :: cexpo2=5d-3, gexp2=1d0+cexpo2
Real(Kind=8), Parameter :: cexpo3=7.5d-1, gexp3=1d0+cexpo3
Integer :: j,k,ny1,ny2,nyedge,nzbot,nztop,nzl,nzl2,nzbotc,nztopc,nzf
Real(Kind=8) :: r1,r2,l1,l2,yedge,zbot,ztop,zbotc,ztopc,zlf,zlc
Real(Kind=8) :: ayin, ayout, azbotc, azinc, aztopc, azf,sheathc

maxtit = 40
modrec = 10
tstep_max = 1d-9
custom_grid = .FALSE.
output_name = "output"
phi_boundary_mode = 0
sig_threshold=1d-3
initial_charge = 0d0
convection_mode = 3
surface_mode = 0
directCGS = .True.
electron_temperature=3d4
rwcount=1
t=0d0
high_charge_dens=.TRUE.
directCGS=.TRUE.
metastable_mode=0
sig_time=1d-6
ifread=.False.
Print *, "readin: Opening input file: ", Trim(AdjustL(input_file))
Open(21, File = Trim(AdjustL(input_file)), Action = "Read", Status = "Old", IOStat=error)
If (error/=0) Then
  Print *, "READIN: Error opening input file: ", Trim(AdjustL(input_file)), error
  Stop
End If
Read(21, NML=inputs)
Close(21, IOStat=error)
If (error/=0) Print *, "READIN: Error closing input file: ", Trim(AdjustL(input_file)), error

If (ifread) Then
  ny=nyo-2
  nz=nzo-2
  mmax=Max(ny,nz)
  Call zero
  Call rwdata(1)
Else
  If (custom_grid) Then
    ny=nyo-2
    nz=nzo-2
    mmax=Max(ny,nz)
    Call zero
    Call rwdata(0)
    Write(*,*) "main:j, dy,y"
    Do j=0,ny+1
      Write(*,*) j, dy(j), y(j)
    End Do
    Write(*,*) "main:k, dz,z"
    Do k=0,nz+1
      Write(*,*) k, dz(k), z(k)
    End Do
    
  Else
    mmax=Max(ny,nz)
    r1 = Min(yl1,yl2)
    r2 = Max(yl1,yl2)
    yedge = ydim-r2
      ny1=0d0
      ny2 = Int(Log(1d0+r2*cexpo2/sheath)/Log(gexp2))+1
      Call multcoef(gexp2,r2,sheath,ny2,ayin)
    nyedge = Int(Log(1d0+yedge*cexpo3/sheath)/Log(gexp3))+1
    Call multcoef(gexp3,yedge,sheath,nyedge,ayout)
    ny = ny2 + nyedge - 2! ny1 +
    Write(*,*) "readin: Debug!"
    Write(*,*) "ny,ny2,nyedge", ny, ny2, nyedge
    l1 = Min(zl1,zl2)
    l2 = Max(zl1,zl2)
    zbot = zh1-l2/2d0
    ztop = zdim-zh1-l2/2d0
    zlf=r2
    zbotc=zbot-zlf
    ztopc=ztop-zlf
    zlc=5d-1*l2-zlf
    nzf = Int(Log(1d0+zlf*cexpo/sheath)/Log(gexp))+1
    Call multcoef(gexp,zlf,sheath,nzf,azf)
    sheathc=sheath*azf**(nzf-1)
    nzbotc=Int(Log(1d0+zbotc*cexpo3/sheathc)/Log(gexp3))+1
    Call multcoef(gexp3,zbotc,sheathc,nzbotc,azbotc)
    nztopc=Int(Log(1d0+ztopc*cexpo3/sheathc)/Log(gexp3))+1
    Call multcoef(gexp3,ztopc,sheathc,nztopc,aztopc)
    nzl2=Int(Log(1d0+zlc*cexpo3/sheathc)/Log(gexp3))+1
    Call multcoef(gexp3,zlc,sheathc,nzl2,azinc)
    nzbot = nzf + nzbotc
    nztop = nzf + nztopc
    nzl = 2*(nzf + nzl2)
    nz = nzbot+nzl+nztop-2
    Write(*,*) "readin: Debug!"
    Write(*,*) "nz,nzbot,nzl,nztop", nz,nzbot,nzl,nztop
    Call zero
    
    ForAll(k=1:ny2) dy(ny2-k)= sheath*(ayin**(k-1))
    ForAll(k=1:nyedge) dy(ny2-1+k)= sheath*(ayout**(k-1))

    y(0) = 0d0
    y(ny+1) = ydim
    ForAll(j=1:ny) y(j) = Sum(dy(0:j-1)) + 5d-1*dy(j) + y(0)
    dy(0)=2d0*dy(0)
    dy(ny+1)=2d0*dy(ny+1)

    ForAll(k=1:nzf) 
      dz(nzbot-k)= sheath*(azf**(k-1))
      dz(nzbot+k-1)= sheath*(azf**(k-1))
      dz(nzbot+nzl-k)= sheath*(azf**(k-1))
      dz(nzbot+nzl+k-1)= sheath*(azf**(k-1))
    End ForAll
    ForAll(k=1:nzbotc) dz(nzbotc-k)= sheathc*(azbotc**(k-1))
    ForAll(k=1:nztopc) dz(nzbot+nzl+nzf+k-1)= sheathc*(aztopc**(k-1))
    ForAll(k=1:nzl2) 
      dz(nzbot+nzf+k-1)= sheathc*(azinc**(k-1))
      dz(nzbot+nzl-nzf-k)= sheathc*(azinc**(k-1))
    End ForAll

    z(0) = 0d0
    z(nz+1) = zdim
    ForAll(k=1:nz) z(k) = Sum(dz(0:k-1)) + 5d-1*dz(k) + z(0)
    dz(0)=2d0*dz(0)
    dz(nz+1)=2d0*dz(nz+1)
    Write(*,*) "main:j, dy,y"
    Do j=0,ny+1
      Write(*,*) j, dy(j), y(j)
    End Do
    Write(*,*) "main:k, dz,z"
    Do k=0,nz+1
      Write(*,*) k, dz(k), z(k)
    End Do
    ian=0
    Do j=0,ny2-1
      Do k=nzbot,nzbot+nzl+1
        If (Sqrt(r2**2-y(j)**2)>=z(k)-zh1-l2/2d0+r2.AND.-Sqrt(r2**2-y(j)**2)<=z(k)-zh1+l2/2d0-r2) ian(j,k)=1
      End Do
    End Do
  End If
End If
Write(*,*) "readin: sum ian", Sum(ian)
Call setco

cp_ion=Abs(sig_cond/(emob_cond*elec_charge))
cp_ion2=Abs(sig_air/(elec_charge*pmob))
If (high_charge_dens) Then
  Where (ian==0) 
    rne = 1d-9*cp_ion2
    rnn = 9.99999999d-1*cp_ion2
    emob = emob_air
  ElseWhere (ian==1) 
    rne = 0d0
    emob = emob_cond
  ElseWhere (ian==2) 
    rne = 0d0
    emob = emob_ins
  End Where
Else
  Where (ian==0) 
    rne = 1d-9*cp_ion2
    rnn = 9.99999999d-1*cp_ion2
    emob = emob_air
  ElseWhere (ian==1) 
    rne = cp_ion
    emob = emob_cond
  ElseWhere (ian==2) 
    rne = 0d0
    emob = emob_ins
  End Where
End If
rnp = rne+rnn

phi=0d0
pc=-1d8

ForAll (j=0:ny+1,k=1:nz+1) phi(j,k) = (pc/z(nz+1))*z(k) !Max(tstep_max/rise_time,1d0)*
Write(*,*) "readin: initial Max rne, rnp", MaxVal(rne), MaxVal(rnp)
Write(*,*) "readin: initial Min rne, rnp", MinVal(rne), MinVal(rnp)

End Subroutine readin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine multcoef(pro,z,fst,n,aa)
Implicit None
Integer, Intent(In) :: n
Real(Kind=8), Intent(In) :: z, fst, pro
Real(Kind=8) :: alp, y1, y2, x, y
Real(Kind=8), Intent(Out) :: aa
Real(Kind=8) :: fun, fun2
Integer :: i
      
alp=z/fst
i=0

If (alp>n) Then
  x = pro 
Else
  x=1d0
End If
y = fun(x,z,n,fst)
y2=y
Do While(Abs(y)>1d-12)
  i=i+1
  y1 = fun2(x,n,fst)
  x = x - y/y1
  y = fun(x,z,n,fst)
  If (y-y2==0d0.OR.i>=128) Exit
  y2=y
End Do

aa=x
End Subroutine multcoef
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Function fun(x,z,n,fst) Result(res_fun)
Implicit None
Integer, Intent(In) :: n
Real(Kind=8), Intent(In) :: x, z, fst
Real(Kind=8) :: res_fun,y
Integer :: i

If (Abs(x-1d0)>0d0) Then
  y=(1d0-x**n)/(1d0-x)
Else
  y=0d0
  If (x>=1d0) Then
    Do i=0,n-1
      y = y + x**i
    End Do
  Else
    Do i=n-1,0,-1
      y = y + x**i
    End Do
  End If
End If
res_fun = fst*y-z
End Function fun
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Function fun2(x,n,fst) Result(res_fun)
Implicit None
Integer, Intent(In) :: n
Real(Kind=8), Intent(In) :: x, fst
Real(Kind=8) :: res_fun,y
Integer :: i

If (Abs(x-1d0)>0d0) Then
  y=(x + (n*(x - 1d0) - x)*x**n)/(x*(x - 1d0)**2)
Else
  y=0d0
  If (x>=1d0) Then
    Do i=1,n-1
      y = y + i*x**(i-1)
    End Do
  Else
    Do i=n-1,1,-1
      y = y + i*x**(i-1)
    End Do
  End If
End If
res_fun = fst*y

End Function fun2