Subroutine amf
Use comm

Implicit None
Real(Kind=8), Parameter :: ethresh=1d-18, ethreshl=5d-20
Real(Kind=8) :: eav,ecy,ecz,temp0,temp1,temp2,dt
Real(Kind=8), Dimension(0:ny+1,0:nz+1) :: eps
Integer :: j,k
alpha=0d0
eta=0d0
dalphay=0d0
dalphaz=0d0
detay=0d0
detaz=0d0
sig=0d0
sigy=0d0
sigz=0d0
rney=0d0
rnez=0d0
rnpy=0d0
rnpz=0d0
rnny=0d0
rnnz=0d0
ey=0d0
ez=0d0
emoby=0d0
emobz=0d0
pmoby=0d0
pmobz=0d0
nmoby=0d0
nmobz=0d0
epsy=0d0
epsz=0d0
eps=0d0

temp_sol=0d0
Do k=0,nz+1
  Do j=0,ny+1
    Select Case(ian(j,k))
    Case(0)
      sig(j,k) = elec_charge*(pmob*rnp(j,k)-(emob_air*rne(j,k)+nmob*rnn(j,k)))
      eps(j,k) = eps_air
    Case(1)
      sig(j,k) = sig_cond
      eps(j,k) = eps_cond
    Case(2)
      sig(j,k) = -elec_charge*emob_ins*rne(j,k)
      eps(j,k) = eps_ins
    Case Default
      Write(*,*) "amf: ERROR, ian can only take values 0-2.", ian(j,k)
      Write(*,*) "j,k", j,k
    End Select
  End Do
End Do

Do k=0,nz+1
  Do j=0,ny+1
    If (j<ny+1) Then
      sigy(j,k) = sig(j,k)*sig(j+1,k)
      If (sigy(j,k)>0d0) sigy(j,k) = sigy(j,k)/(g(j)*sig(j,k) + (1d0-g(j))*sig(j+1,k))
      epsy(j,k) = eps(j,k)*eps(j+1,k)/(g(j)*eps(j,k) + (1d0-g(j))*eps(j+1,k))
      
      ey(j,k) = -2d0*(phi(j+1,k)-phi(j,k))/(dy(j+1)+dy(j))

      If (convection_mode==0) Then
        neycase(j,k) = 0
        npycase(j,k) = 0
        nnycase(j,k) = 0
      Else If (convection_mode==1) Then
        neycase(j,k) = 1
        npycase(j,k) = 1
        nnycase(j,k) = 1
      Else If (convection_mode==2) Then
        neycase(j,k) = 2
        npycase(j,k) = 2
        nnycase(j,k) = 2
      Else
        neycase(j,k) = 0
        npycase(j,k) = 0
        nnycase(j,k) = 0
        If ((ey(j,k)>=0d0.AND.rne(j+1,k)>rne(j,k)).OR.(ey(j,k)<=0d0.AND.rne(j+1,k)<rne(j,k))) neycase(j,k) = 1
        If ((ey(j,k)<=0d0.AND.rnp(j+1,k)>rnp(j,k)).OR.(ey(j,k)>=0d0.AND.rnp(j+1,k)<rnp(j,k))) npycase(j,k) = 1
        If ((ey(j,k)>=0d0.AND.rnn(j+1,k)>rnn(j,k)).OR.(ey(j,k)<=0d0.AND.rnn(j+1,k)<rnn(j,k))) nnycase(j,k) = 1
        If (ian(j,k)/=ian(j+1,k)) Then
          neycase(j,k) = 0
          npycase(j,k) = 0
          nnycase(j,k) = 0
        End If
      End If

      If (high_charge_dens.AND.(ian(j,k)/=0.OR.ian(j+1,k)/=0)) Then
        If (ian(j,k)==1.AND.ian(j+1,k)==1) Then
          rney(j,k)=cp_ion
        Else If (ian(j,k)==1.AND.ian(j+1,k)==0) Then
          rney(j,k) = cp_ion*rne(j+1,k)
          If (rney(j,k)>0d0) Then
            rney(j,k) = rney(j,k)/(g(j)*cp_ion + (1d0-g(j))*rne(j+1,k))
          Else
            rney(j,k) = 0d0
          End If
        Else If (ian(j,k)==0.AND.ian(j+1,k)==1) Then
          rney(j,k) = cp_ion*rne(j,k)
          If (rney(j,k)>0d0) Then
            rney(j,k) = rney(j,k)/(g(j)*rne(j,k) + (1d0-g(j))*cp_ion)
          Else
            rney(j,k) = 0d0
          End If
        Else
          Stop
        End If
      Else
        Select Case(neycase(j,k))
        Case (0)
          rney(j,k) = rne(j,k)*rne(j+1,k)
          If (rney(j,k)>0d0) Then
            rney(j,k) = rney(j,k)/(g(j)*rne(j,k) + (1d0-g(j))*rne(j+1,k))
          Else
            rney(j,k) = 0d0
          End If
        Case (1)
          If (((rne(j,k)<0d0.AND.ian(j,k)==0).OR.(ian(j+1,k)==0.AND.rne(j+1,k)<0d0)) &
               .OR. ((rne(j,k)<0d0.OR.rne(j+1,k)<0d0).AND..NOT.high_charge_dens)) Then
            Write(*,*) j, k
            Write(*,*) rne(j,k), rne(j+1,k)
            Write(*,*) ian(j,k), ian(j+1,k)
            Stop
          End If
          rney(j,k) = (rne(j,k)**(1d0-g(j)))*(rne(j+1,k)**g(j))
        Case(2)
          rney(j,k) = (1d0-g(j))*rne(j,k)+g(j)*rne(j+1,k)
        Case(3)
          rney(j,k) = Min(rne(j,k),rne(j+1,k))
        End Select
      End If

      Select Case(npycase(j,k))
      Case (0)
        rnpy(j,k) = rnp(j,k)*rnp(j+1,k)
        If (rnpy(j,k)>0d0)  Then
          rnpy(j,k) = rnpy(j,k)/(g(j)*rnp(j,k) + (1d0-g(j))*rnp(j+1,k))
        Else
          rnpy(j,k) = 0d0
        End If
      Case (1)
        rnpy(j,k) = (rnp(j,k)**(1d0-g(j)))*(rnp(j+1,k)**g(j))
      Case(2)
        rnpy(j,k) = (1d0-g(j))*rnp(j,k)+g(j)*rnp(j+1,k)
      Case(3)
        rnpy(j,k) = Min(rnp(j,k),rnp(j+1,k))
      End Select

      Select Case(nnycase(j,k))
      Case (0)
        rnny(j,k) = rnn(j,k)*rnn(j+1,k)
        If (rnny(j,k)>0d0)  Then
          rnny(j,k) = rnny(j,k)/(g(j)*rnn(j,k) + (1d0-g(j))*rnn(j+1,k))
        Else
          rnny(j,k) = 0d0
        End If
      Case (1)
        rnny(j,k) = (rnn(j,k)**(1d0-g(j)))*(rnn(j+1,k)**g(j))
      Case(2)
        rnny(j,k) = (1d0-g(j))*rnn(j,k)+g(j)*rnn(j+1,k)
      Case(3)
        rnny(j,k) = Min(rnn(j,k),rnn(j+1,k))
      End Select
      
      If (ian(j,k)==ian(j+1,k)) Then
        Select Case(ian(j,k))
        Case(0)
          emoby(j,k)=emob_air
          pmoby(j,k)=pmob
          nmoby(j,k)=nmob
        Case(1)
          emoby(j,k)=emob_cond
          pmoby(j,k)=0d0
          nmoby(j,k)=0d0
        Case(2)
          emoby(j,k)=emob_ins
          pmoby(j,k)=0d0
          nmoby(j,k)=0d0
        End Select
      Else
        pmoby(j,k)=0d0
        nmoby(j,k)=0d0
        If (surface_mode==0) Then
          emoby(j,k)=0d0
        Else If (surface_mode==1) Then
          emoby(j,k)=emob(j,k)*emob(j+1,k)/(g(j)*emob(j,k) + (1d0-g(j))*emob(j+1,k))
        Else
          If ((ey(j,k)>=0d0.AND.ian(j+1,k)>0).OR.(ey(j,k)<=0d0.AND.ian(j,k)>0)) Then
            emoby(j,k) = 0d0
          Else
            emoby(j,k)=emob(j,k)*emob(j+1,k)/(g(j)*emob(j,k) + (1d0-g(j))*emob(j+1,k))
          End If
        End If
      End If
    End If
    If (k<nz+1) Then
      sigz(j,k) = sig(j,k)*sig(j,k+1)
      If (sigz(j,k)>0d0) sigz(j,k) = sigz(j,k)/(h(k)*sig(j,k) + (1d0-h(k))*sig(j,k+1))
      epsz(j,k) = eps(j,k)*eps(j,k+1)/(h(k)*eps(j,k) + (1d0-h(k))*eps(j,k+1))
      
      ez(j,k) = -2d0*(phi(j,k+1)-phi(j,k))/(dz(k+1)+dz(k))
      If (convection_mode==0) Then
        nezcase(j,k) = 0
        npzcase(j,k) = 0
        nnzcase(j,k) = 0
      Else If (convection_mode==1) Then
        nezcase(j,k) = 1
        npzcase(j,k) = 1
        nnzcase(j,k) = 1
      Else If (convection_mode==2) Then
        nezcase(j,k) = 2
        npzcase(j,k) = 2
        nnzcase(j,k) = 2
      Else
        nezcase(j,k) = 0
        npzcase(j,k) = 0
        nnzcase(j,k) = 0
        If ((ez(j,k)>=0d0.AND.rne(j,k+1)>rne(j,k)).OR.(ez(j,k)<=0d0.AND.rne(j,k+1)<rne(j,k))) nezcase(j,k) = 1
        If ((ez(j,k)<=0d0.AND.rnp(j,k+1)>rnp(j,k)).OR.(ez(j,k)>=0d0.AND.rnp(j,k+1)<rnp(j,k))) npzcase(j,k) = 1
        If ((ez(j,k)>=0d0.AND.rnn(j,k+1)>rnn(j,k)).OR.(ez(j,k)<=0d0.AND.rnn(j,k+1)<rnn(j,k))) nnzcase(j,k) = 1
        If (ian(j,k)/=ian(j,k+1)) Then
          nezcase(j,k) = 0
          npzcase(j,k) = 0
          nnzcase(j,k) = 0
        End If
      End If

      If (high_charge_dens.AND.(ian(j,k)/=0.OR.ian(j,k+1)/=0)) Then
        If (ian(j,k)==1.AND.ian(j,k+1)==1) Then
          rnez(j,k)=cp_ion
        Else If (ian(j,k)==1.AND.ian(j,k+1)==0) Then
          rnez(j,k) = cp_ion*rne(j,k+1)
          If (rnez(j,k)>0d0) Then
            rnez(j,k) = rnez(j,k)/(h(k)*cp_ion + (1d0-h(k))*rne(j,k+1))
          Else
            rnez(j,k) = 0d0
          End If
        Else If (ian(j,k)==0.AND.ian(j,k+1)==1) Then
          rnez(j,k) = cp_ion*rne(j,k)
          If (rnez(j,k)>0d0) Then
            rnez(j,k) = rnez(j,k)/(h(k)*rne(j,k) + (1d0-h(k))*cp_ion)
          Else
            rnez(j,k) = 0d0
          End If
        Else
          Stop
        End If
      Else
        Select Case(nezcase(j,k))
        Case (0)
          rnez(j,k) = rne(j,k)*rne(j,k+1)
          If (rnez(j,k)>0d0) Then
            rnez(j,k) = rnez(j,k)/(h(k)*rne(j,k) + (1d0-h(k))*rne(j,k+1))
          Else
            rnez(j,k) = 0d0
          End If
        Case (1)
          rnez(j,k) = (rne(j,k)**(1d0-h(k)))*(rne(j,k+1)**h(k))
        Case(2)
          rnez(j,k) = (1d0-h(k))*rne(j,k)+h(k)*rne(j,k+1)
        Case(3)
          rnez(j,k) = Min(rne(j,k),rne(j,k+1))
        End Select
      End If

      Select Case(npzcase(j,k))
      Case (0)
        rnpz(j,k) = rnp(j,k)*rnp(j,k+1)
        If (rnpz(j,k)>0d0)  Then
          rnpz(j,k) = rnpz(j,k)/(h(k)*rnp(j,k) + (1d0-h(k))*rnp(j,k+1))
        Else
          rnpz(j,k) = 0d0
        End If
      Case (1)
        rnpz(j,k) = (rnp(j,k)**(1d0-h(k)))*(rnp(j,k+1)**h(k))
      Case(2)
        rnpz(j,k) = (1d0-h(k))*rnp(j,k)+h(k)*rnp(j,k+1)
      Case(3)
        rnpz(j,k) = Min(rnp(j,k),rnp(j,k+1))
      End Select

      Select Case(nnzcase(j,k))
      Case (0)
        rnnz(j,k) = rnn(j,k)*rnn(j,k+1)
        If (rnnz(j,k)>0d0)  Then
          rnnz(j,k) = rnnz(j,k)/(h(k)*rnn(j,k) + (1d0-h(k))*rnn(j,k+1))
        Else
          rnnz(j,k) = 0d0
        End If
      Case (1)
        rnnz(j,k) = (rnn(j,k)**(1d0-h(k)))*(rnn(j,k+1)**h(k))
      Case(2)
        rnnz(j,k) = (1d0-h(k))*rnn(j,k)+h(k)*rnn(j,k+1)
      Case(3)
        rnnz(j,k) = Min(rnn(j,k),rnn(j,k+1))
      End Select

      If (ian(j,k)==ian(j,k+1)) Then
        Select Case(ian(j,k))
        Case(0)
          emobz(j,k)=emob_air
          pmobz(j,k)=pmob
          nmobz(j,k)=nmob
        Case(1)
          emobz(j,k)=emob_cond
          pmobz(j,k)=0d0
          nmobz(j,k)=0d0
        Case(2)
          emobz(j,k)=emob_ins
          pmobz(j,k)=0d0
          nmobz(j,k)=0d0
        End Select
      Else
        pmobz(j,k)=0d0
        nmobz(j,k)=0d0
        If (surface_mode==0) Then
          emobz(j,k)=0d0
        Else If (surface_mode==1) Then
          emobz(j,k)=emob(j,k)*emob(j,k+1)/(h(k)*emob(j,k) + (1d0-h(k))*emob(j,k+1))
        Else
          If ((ez(j,k)>=0d0.AND.ian(j,k+1)>0).OR.(ez(j,k)<=0d0.AND.ian(j,k)>0)) Then
            emobz(j,k) = 0d0
          Else
            emobz(j,k)=emob(j,k)*emob(j,k+1)/(h(k)*emob(j,k) + (1d0-h(k))*emob(j,k+1))
          End If
        End If
      End If
    End If
  End Do
End Do

Do k=1,nz
  Do j=1,ny
    temp0 = -2d0*(dy(j)+dy(j+1))/((dy(j-1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
    temp1 = 2d0*(dy(j+1)-dy(j-1))/((dy(j-1)+dy(j))*(dy(j+1)+dy(j)))
    temp2 = 2d0*(dy(j)+dy(j-1))/((dy(j+1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
    ecy = -(temp0*phi(j-1,k)+temp1*phi(j,k)+temp2*phi(j+1,k))
    temp0 = -2d0*(dz(k)+dz(k+1))/((dz(k-1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
    temp1 = 2d0*(dz(k+1)-dz(k-1))/((dz(k-1)+dz(k))*(dz(k+1)+dz(k)))
    temp2 = 2d0*(dz(k)+dz(k-1))/((dz(k+1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
    ecz = -(temp0*phi(j,k-1)+temp1*phi(j,k)+temp2*phi(j,k+1))
      eav = Sqrt(ecy**2+ecz**2)

    If (ian(j,k)==0) Then
      If ((7.668926684612879d-19*rn)/(eav+1d-60)>6.5d2) Then
        alpha(j,k)=0d0
        temp0=0d0
      Else
        alpha(j,k) = Abs(emob_air)*eav*(2.420131208148039d-20*rn)/Exp((7.668926684612879d-19*rn)/eav)
        temp0=(emob_air*rn*(2.4201312081480392d-20*eav**2 + 1.8559808802430903d-38*eav*rn))/(Exp((7.668926684612879d-19*rn)/eav)*eav**3)
      End If
      eta(j,k) = Abs(emob_air)*eav*(3.1622776601683793d-37*rn)/(1d-80+eav/rn)**7.5d-1
      temp1=(7.905694150420948d-38*emob_air)/(1d-80+eav/rn)**1.75d0

      dalphay(j,k) = ecy*temp0
      dalphaz(j,k) = ecz*temp0
      detay(j,k) = ecy*temp1
      detaz(j,k) = ecz*temp1
    End If
      temp_sol(ijk(j,k))=phi(j,k)
      temp_sol(ijk(j,k)+nn)=Max(rne(j,k),0d0)
      temp_sol(ijk(j,k)+2*nn)=Max(rnp(j,k),0d0)
      temp_sol(ijk(j,k)+3*nn)=Max(rnn(j,k),0d0)
  End Do
End Do
If (metastable_mode/=0) Then
  If (t>=meta_time0) Then
    If (meta_time0>=meta_time1) Then
      Write(*,*) "amf: meta_time0>meta_time1",meta_time0,meta_time1
      Stop
    Else
      dt=(t-meta_time0)/(meta_time1-meta_time0)
    End If
    If(dt<1d0) Then
      Select Case(metastable_mode)
      Case(1)
        Where(ian(1:ny,1:nz)==0)
          alpha = alpha*(1d0+dt*(meta_alpha-1d0))
          dalphay = dalphay*(1d0+dt*(meta_alpha-1d0))
          dalphaz = dalphaz*(1d0+dt*(meta_alpha-1d0))
        End Where
      Case(2)
        Where(ian(1:ny,1:nz)==0)
          eta = eta*(1d0+dt*(meta_eta-1d0))
          detay = detay*(1d0+dt*(meta_eta-1d0))
          detaz = detaz*(1d0+dt*(meta_eta-1d0))
        End Where
      Case(3)
        Where(ian(1:ny,1:nz)==0)
          alpha = alpha*(1d0+dt*(meta_alpha-1d0))
          eta = eta*(1d0+dt*(meta_eta-1d0))
          dalphay = dalphay*(1d0+dt*(meta_alpha-1d0))
          dalphaz = dalphaz*(1d0+dt*(meta_alpha-1d0))
          detay = detay*(1d0+dt*(meta_eta-1d0))
          detaz = detaz*(1d0+dt*(meta_eta-1d0))
        End Where
      Case Default
        Write(*,*) "amf: invalid metastable_mode, reverting to default", metastable_mode
        metastable_mode = 0
      End Select
    Else
      Select Case(metastable_mode)
      Case(1)
        Where(ian(1:ny,1:nz)==0)
          alpha = alpha*meta_alpha
          dalphay = dalphay*meta_alpha
          dalphaz = dalphaz*meta_alpha
        End Where
      Case(2)
        Where(ian(1:ny,1:nz)==0)
          eta = eta*meta_eta
          detay = detay*meta_eta
          detaz = detaz*meta_eta
        End Where
      Case(3)
        Where(ian(1:ny,1:nz)==0)
          alpha = alpha*meta_alpha
          eta = eta*meta_eta
          dalphay = dalphay*meta_alpha
          dalphaz = dalphaz*meta_alpha
          detay = detay*meta_eta
          detaz = detaz*meta_eta
        End Where
      Case Default
        Write(*,*) "amf: invalid metastable_mode, reverting to default", metastable_mode
        metastable_mode = 0
      End Select
    End If
  End If
End If

End Subroutine amf
!************************************************************
!************************************************************

Subroutine amft

Use comm
Implicit None
Real(Kind=8) :: tfact
Integer ::  j, k, l
Logical :: check_time
Write(*,*) "amft: intializing"
Call fill
Write(*,*) "amft: begin fill"

check_time=.True.
If (t<=rise_time) Then
  tstep = Max(Min(tstep_max,5d-2*rise_time,sig_time/MaxVal(sig)),1d-6*tstep_max)
Else
  tstep = Max(Min(2d0*tstep,tstep_max,sig_time/MaxVal(sig)),1d-6*tstep_max)
End If
ttest: Do While (check_time)
check_time=.False.
Write(*,*) "amft: tstep", tstep
  Do k=1,nz
    Do j=1,ny
      tfact = 2d0*vol(j,k)/tstep
      orne(j,k,1) = -tfact
      ornp(j,k,1) = -tfact
      ornn(j,k,1) = -tfact
      orne(j,k,2) = af(ijk(j,k)+nn) - tfact*rne(j,k)
      ornp(j,k,2) = af(ijk(j,k)+2*nn) - tfact*rnp(j,k)
      ornn(j,k,2) = af(ijk(j,k)+3*nn) - tfact*rnn(j,k)
      Do l=1,20
        If (connect(ijk(j,k),l)>0) Then
          orne(j,k,2) = orne(j,k,2) - ap(ijk(j,k)+nn,l)*temp_sol(connect(ijk(j,k),l))
          ornp(j,k,2) = ornp(j,k,2) - ap(ijk(j,k)+2*nn,l)*temp_sol(connect(ijk(j,k),l))
          ornn(j,k,2) = ornn(j,k,2) - ap(ijk(j,k)+3*nn,l)*temp_sol(connect(ijk(j,k),l))
        End If
      End Do
      If ((orne(j,k,2)>0d0.AND.(ian(j,k)==0.OR..NOT.high_charge_dens)) .OR.&
          ornp(j,k,2)>0d0.OR.ornn(j,k,2)>0d0.OR.&
          4d0*Abs(ey(j,k)*emoby(j,k))*tstep>dy(j)+dy(j+1).OR.&
          4d0*Abs(ez(j,k)*emobz(j,k))*tstep>dz(k)+dz(k+1).OR.&
          4d0*Abs(ey(j-1,k)*emoby(j-1,k))*tstep>dy(j-1)+dy(j).OR.&
          4d0*Abs(ez(j,k-1)*emobz(j,k-1))*tstep>dz(k-1)+dz(k).OR.&
          ((tstep*Abs(orne(j,k,2) - orne(j,k,1)*rne(j,k))>rne(j,k)*vol(j,k).AND.(.NOT.high_charge_dens.OR.ian(j,k)==0)).OR.&
          (tstep*Abs(orne(j,k,2) - orne(j,k,1)*rne(j,k))>(rne(j,k)+cp_ion)*vol(j,k).AND.(high_charge_dens.AND.ian(j,k)==1))).OR.&
          tstep*Abs(ornp(j,k,2) - ornp(j,k,1)*rnp(j,k))>rnp(j,k)*vol(j,k).OR.&
          tstep*Abs(ornn(j,k,2) - ornn(j,k,1)*rnn(j,k))>rnn(j,k)*vol(j,k).OR.&
         (orne(j,k,1)+ap(ijk(j,k)+nn,8))>=0d0.OR.&
         (ornp(j,k,1)+ap(ijk(j,k)+2*nn,13))>=0d0.OR.&
         (ornn(j,k,1)+ap(ijk(j,k)+3*nn,18))>=0d0) Then
         
        check_time=.True.
        tstep=5d-1*tstep
        If (tstep<1d-24) Stop
        Write(*,*) "amft: reducing timestep", tstep
        Cycle ttest
      End If
    End Do
  End Do
End Do ttest
End Subroutine amft
