text:
!----------------------------------------------------------------------
! Openpipeflow.org.  If used in your work, please cite
! Willis, A. (2017) SoftwareX 6, 124-127.
! https://doi.org/10.1016/j.softx.2017.05.003 (open access)
!                                      Thanks in advance! Ashley 2019.
!----------------------------------------------------------------------
! Find root f(x) = 0 by Newton method with gmres-hookstep:
!    invert  df(x_n)/dx . s = f(x_n)
! for s subject to constraint |s| < del  (size of 'trust region'),
!    then     x_{n+1} = x_n - s
! To define subroutines f,df below, write as
!    df . s = f
! where f is a function to be minimised and operator df provides the
! action of the Jacobian, but may also contain any constraints on the
! update s (e.g. no shifts in a homogenous direction), arranged such
! that the rhs of the constraint (in f) is 0 when converged.
!----------------------------------------------------------------------
!
! Allocate(new_x(1:n),new_fx(1:n)), then put initial guess in new_x(1:n).
! This variable is updated with the current best solution each iteration.
!
! f       evaluates y=f(x):  call f(x, y)
! df      evaluates y=Jacobian.x, see GMRESm.f90:  call df(x, y)
! sub     parameterless subroutine called at end of each iteration,
!         so that data could be saved etc.
! inprod  dot product:  d = inprod(a,b)
! m	  gmres dimension (also max num gmres its)
! n 	  dimension of x
! gtol    tolerence for gmres, typically 0.001
! tol     request |f(x)|<tol
! del     initial size of trust region for usefulness of df
!         set = 0d0 for no hookstep, <0d0 for del = |f(x0)|/10.
! mndl    min size of trust region
! mxdl    max size of trust region, if del<0 set no upper limit
! nits	  max num Newton its
! info	  on input if =1 print* details
!         on exit if =0 sucessful
!                    =1 gmres failed
!                    =2 reached max iterations, nits
!                    =3 trust region got too small
!                    =4 unable to allocate memory
!							A.P.Willis 2008
!----------------------------------------------------------------------

 subroutine newtonhook(f, df, sub, inprod, m, n, gtol, &
                       tol, del, mndl, mxdl, nits, info,out)

   external                        :: f, df, sub
   real(dp), external      :: inprod
   integer(i4),          intent(in)    :: m, n
   real(dp), intent(in)    :: gtol
   real(dp), intent(inout) :: tol, del, mndl, mxdl
   integer(i4),          intent(in)    :: nits
   integer(i4),          intent(inout) :: info
   integer(i4), intent(in) :: out
   real(dp), allocatable :: v(:,:)
   real(dp) :: tol_,x_(n),fx_(n), tol__,x__(n),fx__(n),del__
   real(dp) :: mxdl_, ared, pred, snrm, s(n)
   real(dp) :: gres, gdel, h((m+1)*m)
   integer(i4) :: ginfo

   real(dp) :: normOld, normNew

   write(out, *) 'beginning newton'
   flush(out)
   allocate(v(n,m+1), stat=ginfo)
   if(ginfo/=0) then
      if(info==1) write(out, *)  'newton: unable to allocate memory'
      info = 4
      return
   end if

   new_nits = 0
   new_gits = 0
   new_del  = del
   mxdl_    = mxdl
   ginfo    = info
   call f(new_x, new_fx)
   new_tol  = sqrt(inprod(new_fx,new_fx))
   if(del<0.0_dp)  new_del = new_tol / 10.0_dp
   if(del<0.0_dp)  mxdl_   = 1.0e99_dp
    if(info==1)  then
        write(out, *) 'newton: nits=',new_nits,' res=', new_tol
    end if
   call sub()
   x_   = new_x
   fx_  = new_fx
   tol_ = new_tol
   tol__ = 1.0e99_dp

   if(new_tol<tol) then
    if(info==1)  then
        write(out, *) 'newton: input already converged'
        flush(out)
    end if
      info = 0
      deallocate(v)
      return
   end if

   normOld = sqrt(inprod(new_x,new_x))
   del__ = new_del
    ! - - - - Start main loop - - - - -  -
    do while(.true.)

        if(new_del<mndl) then
            if(info==1)  then
                write(out, *)  'newton: trust region too small'
                flush(out)
            end if
            info = 3
            deallocate(v)
            return
        end if
                     ! find hookstep s and update x
      s        = 0
      gres     = gtol * new_tol
      gdel     = new_del
      if(ginfo/=2) new_gits = m
      if(abs(del)<small) new_gits = 9999
      call gmresm(m,n,s,fx_,df,inprod,h,v,gres,gdel,new_gits,ginfo,out)
      ginfo = info
      new_x = x_ - s
                     ! calc new norm, compare with prediction
      call f(new_x, new_fx)
      new_tol = sqrt(inprod(new_fx,new_fx))
      snrm = sqrt(inprod(s,s))
      ared = tol_ - new_tol
      pred = tol_ - gdel

      if(info==1) then
         write(out, *) 'newton: nits=',new_nits,' res=', new_tol
         write(out, *) 'newton: gits=',new_gits,' del=', new_del
         write(out, *) 'newton: |s|=', snrm,' pred=', pred
         write(out, *) 'newton: ared/pred=', ared/pred
      end if

      if(abs(del)<small) then
        if(info==1) write(out, *) 'newton: took full newton step'
      else if(new_tol>tol__) then
         if(info==1) write(out, *)  'newton: accepting previous step'
         new_x   = x__
         new_fx  = fx__
         new_tol = tol__
         new_del = del__
      else if(ared<0.0_dp) then
         if(info==1) write(out, *)  'newton: norm increased, try smaller step'
         new_del = snrm * 0.5_dp
         ginfo   = 2
      else if(ared/pred<0.75_dp) then
         if(info==1)  write(out, *)  'newton: step ok, trying smaller step'
         x__     = new_x
         fx__    = new_fx
         tol__   = new_tol
         if(ared/pred> 0.1_dp)  del__ = snrm
         if(ared/pred<=0.1_dp)  del__ = snrm*0.5_dp
         new_del = snrm * 0.7_dp
         ginfo   = 2
      else if(snrm<new_del*0.9_dp) then
         if(info==1) write(out, *)  'newton: step good, took full newton step'
         new_del = min(mxdl_,snrm*2.0_dp)
      else if(new_del<mxdl_*0.9_dp) then
         if(info==1) write(out, *)  'newton: step good, trying larger step'
         x__     = new_x
         fx__    = new_fx
         tol__   = new_tol
         del__   = new_del
         new_del = min(mxdl_,snrm*2.0_dp)
         ginfo   = 2
      end if
                         ! check if need to try another s
      if(ginfo==2) cycle
                         ! end of iteration
      new_nits = new_nits + 1
      call sub()
      x_   = new_x
      fx_  = new_fx
      tol_ = new_tol
      tol__ = 1.0e99_dp

      ! Re-compute parameters after each call of sub()
      ! Otherwise when Newton stops is not consistent with the relative error
      ! requirement in general
      normNew = sqrt(inprod(new_x,new_x))
      if(abs(normNew)<small) normNew=1.0_dp
      tol = (tol / normOld) * normNew
      del = (del / normOld) * normNew
      mndl = (mndl / normOld) * normNew
      mxdl = (mxdl / normOld) * normNew
      normOld = normNew

      new_del = del
      mxdl_ = mxdl
      if(del<0.0_dp)  new_del = new_tol / 10.0_dp
      if(del<0.0_dp)  mxdl_   = 1.0e99_dp

      if(new_tol<tol) then
         if(info==1) write(out, *)  'newton: converged'
         info = 0
         deallocate(v)
         return
      else if(new_nits==nits) then
         if(info==1) write(out, *)  'newton: reached max its'
         info = 2
         deallocate(v)
         return
      end if

   end do

 end subroutine newtonhook

max_stars_repo_name: gokhanyalniz/dnsbox

max_issues_repo_name: gokhanyalniz/dnsbox

max_forks_repo_name: gokhanyalniz/dnsbox


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!-----------------------------------------------------------------------
!
!     This module defines data structures for ten Tusscher-Panfilov
!     epicardial cellular activation model for cardiac electrophysiology
!
!-----------------------------------------------------------------------

      MODULE TTPMOD
      USE TYPEMOD
      IMPLICIT NONE

      CONTAINS
!-----------------------------------------------------------------------
      SUBROUTINE TTP_INIT(imyo, nX, nG, X, Xg, X0, Xg0)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX, nG
      REAL(KIND=RKIND), INTENT(OUT) :: X(nX), Xg(nG)
      REAL(KIND=RKIND), INTENT(IN), OPTIONAL :: X0(nX), Xg0(nG)

      SELECT CASE (imyo)
      CASE (1) ! EPI !
!        Initialize state variables
         X(1)   = -85.23_RKIND      ! V      (units: mV)
         X(2)   =  136.89_RKIND     ! K_i    (units: mM)
         X(3)   =  8.6040_RKIND     ! Na_i   (units: mM)
         X(4)   =  1.26E-4_RKIND    ! Ca_i   (units: mM)
         X(5)   =  3.6E-4_RKIND     ! Ca_ss  (units: mM)
         X(6)   =  3.64_RKIND       ! Ca_sr  (units: mM)
         X(7)   =  0.9073_RKIND     ! R'     (dimensionless)

!        Initialize gating variables
         Xg(1)  =  6.21E-3_RKIND    ! x_r1   (dimensionless)
         Xg(2)  =  0.4712_RKIND     ! x_r2   (dimensionless)
         Xg(3)  =  9.5E-3_RKIND     ! x_s    (dimensionless)
         Xg(4)  =  1.72E-3_RKIND    ! m      (dimensionless)
         Xg(5)  =  0.7444_RKIND     ! h      (dimensionless)
         Xg(6)  =  0.7045_RKIND     ! j      (dimensionless)
         Xg(7)  =  3.373E-5_RKIND   ! d      (dimensionless)
         Xg(8)  =  0.7888_RKIND     ! f      (dimensionless)
         Xg(9)  =  0.9755_RKIND     ! f_2    (dimensionless)
         Xg(10) =  0.9953_RKIND     ! f_cass (dimensionless)
         Xg(11) =  0.999998_RKIND   ! s      (dimensionless)
         Xg(12) =  2.42E-8_RKIND    ! r      (dimensionless)

      CASE (2) ! ENDO !
!        Initialize state variables
         X(1)   = -86.709_RKIND     ! V      (units: mV)
         X(2)   =  138.4_RKIND      ! K_i    (units: mM)
         X(3)   =  10.355_RKIND     ! Na_i   (units: mM)
         X(4)   =  1.3E-4_RKIND    ! Ca_i   (units: mM)
         X(5)   =  3.6E-4_RKIND     ! Ca_ss  (units: mM)
         X(6)   =  3.715_RKIND      ! Ca_sr  (units: mM)
         X(7)   =  0.9068_RKIND     ! R'     (dimensionless)

!        Initialize gating variables
         Xg(1)  =  4.48E-3_RKIND    ! x_r1   (dimensionless)
         Xg(2)  =  0.476_RKIND      ! x_r2   (dimensionless)
         Xg(3)  =  8.7E-3_RKIND     ! x_s    (dimensionless)
         Xg(4)  =  1.55E-3_RKIND    ! m      (dimensionless)
         Xg(5)  =  0.7573_RKIND     ! h      (dimensionless)
         Xg(6)  =  0.7225_RKIND     ! j      (dimensionless)
         Xg(7)  =  3.164E-5_RKIND   ! d      (dimensionless)
         Xg(8)  =  0.8009_RKIND     ! f      (dimensionless)
         Xg(9)  =  0.9778_RKIND     ! f_2    (dimensionless)
         Xg(10) =  0.9953_RKIND     ! f_cass (dimensionless)
         Xg(11) =  0.3212_RKIND     ! s      (dimensionless)
         Xg(12) =  2.235E-8_RKIND   ! r      (dimensionless)

      CASE (3) ! MID-MYO !
!        Initialize state variables
         X(1)   = -85.423_RKIND     ! V      (units: mV)
         X(2)   =  138.52_RKIND     ! K_i    (units: mM)
         X(3)   =  10.132_RKIND     ! Na_i   (units: mM)
         X(4)   =  1.53E-4_RKIND    ! Ca_i   (units: mM)
         X(5)   =  4.2E-4_RKIND     ! Ca_ss  (units: mM)
         X(6)   =  4.272_RKIND      ! Ca_sr  (units: mM)
         X(7)   =  0.8978_RKIND     ! R'     (dimensionless)

!        Initialize gating variables
         Xg(1)  =  1.65E-2_RKIND    ! x_r1   (dimensionless)
         Xg(2)  =  0.4730_RKIND     ! x_r2   (dimensionless)
         Xg(3)  =  1.74E-2_RKIND    ! x_s    (dimensionless)
         Xg(4)  =  1.65E-3_RKIND    ! m      (dimensionless)
         Xg(5)  =  0.7490_RKIND     ! h      (dimensionless)
         Xg(6)  =  0.6788_RKIND     ! j      (dimensionless)
         Xg(7)  =  3.288E-5_RKIND   ! d      (dimensionless)
         Xg(8)  =  0.7026_RKIND     ! f      (dimensionless)
         Xg(9)  =  0.9526_RKIND     ! f_2    (dimensionless)
         Xg(10) =  0.9942_RKIND     ! f_cass (dimensionless)
         Xg(11) =  0.999998_RKIND   ! s      (dimensionless)
         Xg(12) =  2.347E-8_RKIND   ! r      (dimensionless)

      END SELECT

      IF (PRESENT(X0)) THEN
         IF (SIZE(X0,1) .NE. nX) THEN
            STOP "ERROR: inconsistent array size (TTP initialization)"
         END IF
         X(:)  = X0(:)
      END IF

      IF (PRESENT(Xg0)) THEN
         IF (SIZE(Xg0,1) .NE. nG) THEN
            STOP "ERROR: inconsistent array size (TTP initialization)"
         END IF
         Xg(:) = Xg0(:)
      END IF

      RETURN
      END SUBROUTINE TTP_INIT
!-----------------------------------------------------------------------
!     Time integration performed using Forward Euler method
      SUBROUTINE TTP_INTEGFE(imyo, nX, nG, X, Xg, Ts, dt, Istim, Ksac,
     2   RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX, nG
      REAL(KIND=RKIND), INTENT(INOUT) :: X(nX), Xg(nG), RPAR(18)
      REAL(KIND=RKIND), INTENT(IN) :: Ts, dt, Istim, Ksac

      REAL(KIND=RKIND) :: f(nX)

!     Get time derivatives (RHS)
      CALL TTP_GETF(imyo, nX, nG, X, Xg, f, Istim, Ksac, RPAR)

!     Update gating variables
      CALL TTP_UPDATEG(imyo, dt, nX, nG, X, Xg)

!     Update state variables
      X = X + dt*f(:)

      RETURN
      END SUBROUTINE TTP_INTEGFE
!-----------------------------------------------------------------------
!     Time integration performed using 4th order Runge-Kutta method
      SUBROUTINE TTP_INTEGRK(imyo, nX, nG, X, Xg, Ts, dt, Istim, Ksac,
     2   RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX, nG
      REAL(KIND=RKIND), INTENT(INOUT) :: X(nX), Xg(nG), RPAR(18)
      REAL(KIND=RKIND), INTENT(IN) :: Ts, dt, Istim, Ksac

      REAL(KIND=RKIND) :: dt6, Xrk(nX), Xgr(nG), frk(nX,4)

      dt6 = dt/6._RKIND

!     RK4: 1st pass
      Xrk = X
      CALL TTP_GETF(imyo, nX, nG, Xrk, Xg, frk(:,1), Istim, Ksac,
     2   RPAR)

!     Update gating variables by half-dt
      Xgr = Xg
      CALL TTP_UPDATEG(imyo, 0.5_RKIND*dt, nX, nG, X, Xgr)

!     RK4: 2nd pass
      Xrk = X + 0.5_RKIND*dt*frk(:,1)
      CALL TTP_GETF(imyo, nX, nG, Xrk, Xgr, frk(:,2), Istim, Ksac,
     2   RPAR)

!     RK4: 3rd pass
      Xrk = X + 0.5_RKIND*dt*frk(:,2)
      CALL TTP_GETF(imyo, nX, nG, Xrk, Xgr, frk(:,3), Istim, Ksac,
     2   RPAR)

!     Update gating variables by full-dt
      Xgr = Xg
      CALL TTP_UPDATEG(imyo, dt, nX, nG, X, Xgr)

!     RK4: 4th pass
      Xrk = X + dt*frk(:,3)
      CALL TTP_GETF(imyo, nX, nG, Xrk, Xgr, frk(:,4), Istim, Ksac,
     2   RPAR)

      X = X + dt6*(frk(:,1) + 2._RKIND*(frk(:,2) + frk(:,3)) + frk(:,4))
      Xg = Xgr

      RETURN
      END SUBROUTINE TTP_INTEGRK
!-----------------------------------------------------------------------
!     Compute currents and time derivatives of state variables
      SUBROUTINE TTP_GETF(i, nX, nG, X, Xg, dX, I_stim, K_sac, RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: i, nX, nG
      REAL(KIND=RKIND), INTENT(IN) :: X(nX), Xg(nG), I_stim, K_sac
      REAL(KIND=RKIND), INTENT(OUT) :: dX(nX)
      REAL(KIND=RKIND), INTENT(INOUT) :: RPAR(18)

      INCLUDE "PARAMS_TTP.f"

      REAL(KIND=RKIND) :: RT, a, b, tau, sq5, e1, e2, e3, e4, n1, n2,
     2   d1, d2, d3, I_sac

!     Local copies of state variables
      V     = X(1)
      K_i   = X(2)
      Na_i  = X(3)
      Ca_i  = X(4)
      Ca_ss = X(5)
      Ca_sr = X(6)
      R_bar = X(7)

!     Local copies of gating variables
      xr1   = Xg(1)
      xr2   = Xg(2)
      xs    = Xg(3)
      m     = Xg(4)
      h     = Xg(5)
      j     = Xg(6)
      d     = Xg(7)
      f     = Xg(8)
      f2    = Xg(9)
      fcass = Xg(10)
      s     = Xg(11)
      r     = Xg(12)

!     Stretch-activated currents
      I_sac = K_sac * (Vrest - V)

!      Diff = 1._RKIND / (1.0D1 * rho * Cm * sV)
      RT   = Rc * Tc / Fc
      E_K  = RT * LOG(K_o/K_i)
      E_Na = RT * LOG(Na_o/Na_i)
      E_Ca = 0.5_RKIND * RT * LOG(Ca_o/Ca_i)
      E_Ks = RT * LOG( (K_o + p_KNa*Na_o)/(K_i + p_KNa*Na_i) )

!     I_Na: Fast sodium current
      I_Na = G_Na * (m**3._RKIND) * h * j * (V - E_Na)

!     I_to: transient outward current
      I_to = G_to(i) * r * s * (V - E_K)

!     I_K1: inward rectifier outward current
      e1   = EXP(0.06_RKIND*(V - E_K - 200._RKIND))
      e2   = EXP(2.E-4_RKIND*(V - E_K + 100._RKIND))
      e3   = EXP(0.1_RKIND*(V - E_K - 10._RKIND))
      e4   = EXP(-0.5_RKIND*(V - E_K))
      a    = 0.1_RKIND/(1._RKIND + e1)
      b    = (3._RKIND*e2 + e3) / (1._RKIND + e4)
      tau  = a / (a + b)
      sq5  = SQRT(K_o/5.4_RKIND)
      I_K1 = G_K1 * sq5 * tau * (V - E_K)

!     I_Kr: rapid delayed rectifier current
      I_Kr = G_Kr * sq5 * xr1 * xr2 * (V - E_K)

!     I_Ks: slow delayed rectifier current
      I_Ks = G_Ks(i) * (xs**2._RKIND) * (V - E_Ks)

!     I_CaL: L-type Ca current
      a     = 2._RKIND*(V-15._RKIND)/RT
      b     = 2._RKIND*a*Fc * (0.25_RKIND*Ca_ss*EXP(a) - Ca_o) /
     2   (EXP(a)-1._RKIND)
      I_CaL = G_CaL * d * f * f2 * fcass * b

!     I_NaCa: Na-Ca exchanger current
      e1     = EXP(gamma*V/RT)
      e2     = EXP((gamma-1._RKIND)*V/RT)
      n1     = e1*(Na_i**3._RKIND)*Ca_o - e2*(Na_o**3._RKIND)*Ca_i*alpha
      d1     = K_mNai**3._RKIND + Na_o**3._RKIND
      d2     = K_mCa + Ca_o
      d3     = 1._RKIND + K_sat*e2
      I_NaCa = K_NaCa * n1 / (d1*d2*d3)

!     I_NaK: Na-K pump current
      e1    = EXP(-0.1_RKIND*V/RT)
      e2    = EXP(-V/RT)
      n1    = P_NaK * K_o * Na_i
      d1    = K_o + K_mK
      d2    = Na_i + K_mNa
      d3    = 1._RKIND + 0.1245_RKIND*e1 + 0.0353_RKIND*e2
      I_NaK = n1 / (d1*d2*d3)

!     I_pCa: plateau Ca current
      I_pCa = G_pCa * Ca_i / (K_pCa + Ca_i)

!     I_pK: plateau K current
      I_pK  = G_pK * (V-E_K) /
     2   (1._RKIND + EXP((25._RKIND-V)/5.98_RKIND))

!     I_bCa: background Ca current
      I_bCa = G_bCa * (V - E_Ca)

!     I_bNa: background Na current
      I_bNa = G_bNa * (V - E_Na)

!     I_leak: Sacroplasmic Reticulum Ca leak current
      I_leak = V_leak * (Ca_sr - Ca_i)

!     I_up: Sacroplasmic Reticulum Ca pump current
      I_up  = Vmax_up / (1._RKIND + (K_up/Ca_i)**2._RKIND)

!     I_rel: Ca induced Ca current (CICR)
      k_casr = max_sr - ((max_sr-min_sr)/
     2   (1._RKIND + (EC/Ca_sr)**2._RKIND) )
      k1     = k1p / k_casr
      O      = k1 * R_bar * (Ca_ss**2._RKIND) /
     2   (k3 + k1*(Ca_ss**2._RKIND))
      I_rel  = V_rel * O * (Ca_sr - Ca_ss)

!     I_xfer: diffusive Ca current between Ca subspae and cytoplasm
      I_xfer = V_xfer * (Ca_ss - Ca_i)

!-----------------------------------------------------------------------
!     Now compute time derivatives
!     dV/dt: rate of change of transmembrane voltage
      dX(1)  = -(I_Na + I_to + I_K1 + I_Kr + I_Ks + I_CaL + I_NaCa +
     2   I_NaK + I_pCa + I_pK + I_bCa + I_bNa  + I_stim) + I_sac

!     dK_i/dt
      dX(2)  = -(Cm/(V_c*Fc)) * (I_K1 + I_to + I_Kr + I_Ks + I_pK -
     2   2._RKIND*I_NaK + I_stim)

!     dNa_i/dt
      dX(3)  = -(Cm/(V_c*Fc)) * (I_Na + I_bNa +
     2   3._RKIND*(I_NaK + I_NaCa))

!     dCa_i/dt
      n1     = (I_leak - I_up)*V_sr/V_c + I_xfer
      n2     = -(Cm/(V_c*Fc)) * (I_bCa + I_pCa - 2._RKIND*I_Naca)
     2   / 2._RKIND
      d1     = 1._RKIND + K_bufc*Buf_c/(Ca_i + K_bufc)**2._RKIND
      dX(4)  = (n1 + n2)/d1

!     dCa_ss: rate of change of Ca_ss
      n1     = (-I_CaL*Cm/(2._RKIND*Fc) + I_rel*V_sr - V_c*I_xfer)/V_ss
      d1     = 1._RKIND + K_bufss*Buf_ss/(Ca_ss + K_bufss)**2._RKIND
      dX(5)  = n1 / d1

!     dCa_sr: rate of change of Ca_sr
      n1     = I_up - I_leak - I_rel
      d1     = 1._RKIND + K_bufsr*Buf_sr/(Ca_sr + K_bufsr)**2._RKIND
      dX(6)  = n1 / d1

!     Rbar: ryanodine receptor
      k2     = k2p * k_casr
      dX(7)  = -k2*Ca_ss*R_bar + k4*(1._RKIND - R_bar)

!     Quantities to be written to file
      RPAR(3)  = I_Na
      RPAR(4)  = I_K1
      RPAR(5)  = I_to
      RPAR(6)  = I_Kr
      RPAR(7)  = I_Ks
      RPAR(8)  = I_CaL
      RPAR(9)  = I_NaCa
      RPAR(10) = I_NaK
      RPAR(11) = I_pCa
      RPAR(12) = I_pK
      RPAR(13) = I_bCa
      RPAR(14) = I_bNa
      RPAR(15) = I_leak
      RPAR(16) = I_up
      RPAR(17) = I_rel
      RPAR(18) = I_xfer

      RETURN
      END SUBROUTINE TTP_GETF
!-----------------------------------------------------------------------
!     Update all the gating variables
      SUBROUTINE TTP_UPDATEG(i, dt, n, nG, X, Xg)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: i, n, nG
      REAL(KIND=RKIND), INTENT(IN) :: dt, X(n)
      REAL(KIND=RKIND), INTENT(INOUT) :: Xg(nG)

      INCLUDE "PARAMS_TTP.f"

      REAL(KIND=RKIND) :: a, b, c, tau

      V     = X(1)
      Ca_ss = X(5)

      xr1   = Xg(1)
      xr2   = Xg(2)
      xs    = Xg(3)
      m     = Xg(4)
      h     = Xg(5)
      j     = Xg(6)
      d     = Xg(7)
      f     = Xg(8)
      f2    = Xg(9)
      fcass = Xg(10)
      s     = Xg(11)
      r     = Xg(12)

!     xr1: activation gate for I_Kr
      xr1i   = 1._RKIND/(1._RKIND + EXP(-(26._RKIND+V)/7._RKIND))
      a      = 450._RKIND/(1._RKIND + EXP(-(45._RKIND+V)/10._RKIND))
      b      = 6._RKIND/(1._RKIND + EXP((30._RKIND+V)/11.5_RKIND))
      tau    = a*b
      Xg(1)  = xr1i - (xr1i - xr1)*EXP(-dt/tau)

!     xr2: inactivation gate for I_Kr
      xr2i   = 1._RKIND /(1._RKIND + EXP((88._RKIND+V)/24._RKIND))
      a      = 3._RKIND /(1._RKIND + EXP(-(60._RKIND+V)/20._RKIND))
      b      = 1.12_RKIND/(1._RKIND + EXP(-(60._RKIND-V)/20._RKIND))
      tau    = a*b
      Xg(2)  = xr2i - (xr2i - xr2)*EXP(-dt/tau)

!     xs: activation gate for I_Ks
      xsi    = 1._RKIND/(1._RKIND + EXP(-(5._RKIND+V)/14._RKIND))
      a      = 1400._RKIND/SQRT(1._RKIND + EXP((5._RKIND-V)/6._RKIND))
      b      = 1._RKIND/(1._RKIND + EXP((V-35._RKIND)/15._RKIND))
      tau    = a*b + 80._RKIND
      Xg(3)  = xsi - (xsi - xs)*EXP(-dt/tau)

!     m: activation gate for I_Na
      mi     = 1._RKIND/( (1._RKIND +
     2   EXP(-(56.86_RKIND+V)/9.03_RKIND))**2._RKIND )
      a      = 1._RKIND/(1._RKIND + EXP(-(60._RKIND+V)/5._RKIND))
      b      = 0.1_RKIND/(1._RKIND + EXP((35._RKIND+V)/5._RKIND))
     2       + 0.1_RKIND/(1._RKIND + EXP((V-50._RKIND)/200._RKIND))
      tau    = a*b
      Xg(4)  = mi - (mi - m)*EXP(-dt/tau)

!     h: fast inactivation gate for I_Na
      hi     = 1._RKIND/( (1._RKIND
     2       + EXP((71.55_RKIND+V)/7.43_RKIND))**2._RKIND )
      IF (V .GE. -40._RKIND) THEN
         a   = 0._RKIND
         b   = 0.77_RKIND/(0.13_RKIND*(1._RKIND
     2       + EXP(-(10.66_RKIND+V)/11.1_RKIND)))
      ELSE
         a   = 5.7E-2_RKIND*EXP(-(80._RKIND+V)/6.8_RKIND)
         b   = 2.7_RKIND*EXP(0.079_RKIND*V)
     2       + 310000._RKIND*EXP(0.3485_RKIND*V)
      END IF
      tau    = 1._RKIND / (a + b)
      Xg(5)  = hi - (hi - h)*EXP(-dt/tau)

!     j: slow inactivation gate for I_Na
      ji     = 1._RKIND/( (1._RKIND
     2       + EXP((71.55_RKIND+V)/7.43_RKIND))**2._RKIND )
      IF (V .GE. -40._RKIND) THEN
         a   = 0._RKIND
         b   = 0.6_RKIND*EXP(5.7E-2_RKIND*V)
     2       / (1._RKIND + EXP(-0.1_RKIND*(V+32._RKIND)))
      ELSE
         a   = -(25428._RKIND*EXP(0.2444_RKIND*V)
     2       + 6.948E-6_RKIND*EXP(-0.04391_RKIND*V)) * (V+37.78_RKIND)
     3       / (1._RKIND + EXP(0.311_RKIND*(79.23_RKIND+V)))
         b   = 0.02424_RKIND*EXP(-0.01052_RKIND*V) /
     2         (1._RKIND + EXP(-0.1378_RKIND*(40.14_RKIND+V)))
      END IF
      tau    = 1._RKIND / (a + b)
      Xg(6)  = ji - (ji - j)*EXP(-dt/tau)

!     d: activation gate for I_CaL
      di     = 1._RKIND/(1._RKIND + EXP(-(8._RKIND+V)/7.5_RKIND))
      a      = 1.4_RKIND/(1._RKIND + EXP(-(35._RKIND+V)/13._RKIND))
     2       + 0.25_RKIND
      b      = 1.4_RKIND/(1._RKIND + EXP((5._RKIND+V)/5._RKIND))
      c      = 1._RKIND/(1._RKIND + EXP((50._RKIND-V)/20._RKIND))
      tau    = a*b + c
      Xg(7)  = di - (di - d)*EXP(-dt/tau)

!     f: slow inactivation gate for I_CaL
      fi     = 1._RKIND/(1._RKIND + EXP((20._RKIND+V)/7._RKIND))
      a      = 1102.5_RKIND*EXP(-((V+27._RKIND)**2._RKIND)/225._RKIND)
      b      = 200._RKIND/(1._RKIND + EXP((13._RKIND-V)/10._RKIND))
      c      = 180._RKIND/(1._RKIND + EXP((30._RKIND+V)/10._RKIND))
     2       + 20._RKIND
      tau    = a + b + c
c!     for spiral wave breakup
c      IF (V .GT. 0._RKIND) tau = tau*2._RKIND
      Xg(8)  = fi - (fi - f)*EXP(-dt/tau)

!     f2: fast inactivation gate for I_CaL
      f2i    = 0.67_RKIND/(1._RKIND + EXP((35._RKIND+V)/7._RKIND))
     2       + 0.33_RKIND
      a      = 562._RKIND*EXP(-((27._RKIND+V)**2._RKIND) /240._RKIND)
      b      = 31._RKIND/(1._RKIND + EXP((25._RKIND-V)/10._RKIND))
      c      = 80._RKIND/(1._RKIND + EXP((30._RKIND+V)/10._RKIND))
      tau    = a + b + c
      Xg(9)  = f2i - (f2i - f2)*EXP(-dt/tau)

!     fCass: inactivation gate for I_CaL into subspace
      c      = 1._RKIND/(1._RKIND + (Ca_ss/0.05_RKIND)**2._RKIND)
      fcassi = 0.6_RKIND*c  + 0.4_RKIND
      tau    = 80._RKIND*c + 2._RKIND
      Xg(10) = fcassi - (fcassi - fcass)*EXP(-dt/tau)

!     s: inactivation gate for I_to
      IF (i.EQ.1 .OR. i.EQ.3) THEN
         si  = 1._RKIND/(1._RKIND + EXP((20._RKIND+V)/5._RKIND))
         tau = 85._RKIND*EXP(-((V+45._RKIND)**2._RKIND) /320._RKIND)
     2       + 5._RKIND/(1._RKIND+EXP((V-20._RKIND)/5._RKIND))
     3       + 3._RKIND
      ELSE IF (i .EQ. 2) THEN
         si  = 1._RKIND/(1._RKIND + EXP((28._RKIND+V)/5._RKIND))
         tau = 1000._RKIND*EXP(-((V+67._RKIND)**2._RKIND) /1000._RKIND)
     2       + 8._RKIND
      END IF
      Xg(11) = si - (si - s)*EXP(-dt/tau)

!     r: activation gate for I_to
      ri     = 1._RKIND/(1._RKIND + EXP((20._RKIND-V)/6._RKIND))
      tau    = 9.5_RKIND*EXP(-((V+40._RKIND)**2._RKIND) /1800._RKIND)
     2       + 0.8_RKIND
      Xg(12) = ri - (ri - r)*EXP(-dt/tau)

      RETURN
      END SUBROUTINE TTP_UPDATEG
!-----------------------------------------------------------------------
!     Time integration performed using Crank-Nicholson method
      SUBROUTINE TTP_INTEGCN2(imyo, nX, nG, Xn, Xg, Ts, dt, Istim,
     2   Ksac, IPAR, RPAR)
      USE MATFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX, nG
      INTEGER(KIND=IKIND), INTENT(INOUT) :: IPAR(2)
      REAL(KIND=RKIND), INTENT(INOUT) :: Xn(nX), Xg(nG), RPAR(18)
      REAL(KIND=RKIND), INTENT(IN) :: Ts, dt, Istim, Ksac

      REAL(KIND=RKIND), PARAMETER :: eps = EPSILON(eps)

      INTEGER(KIND=IKIND) :: i, k, itMax
      LOGICAL :: l1, l2, l3
      REAL(KIND=RKIND) :: t, atol, rtol, Xk(nX), fn(nX), fk(nX), rK(nX),
     2   Im(nX,nX), JAC(nX,nX), rmsA, rmsR

      itMax = IPAR(1)
      atol  = RPAR(1)
      rtol  = RPAR(2)

      Im = MAT_ID(nX)
      CALL TTP_GETF(imyo, nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR)

      k  = 0
      Xk = Xn
      l1 = .FALSE.
      l2 = .FALSE.
      l3 = .FALSE.
      t  = Ts + dt
      DO
         k = k + 1

         CALL TTP_GETF(imyo, nX, nG, Xk, Xg, fk, Istim, Ksac, RPAR)
         rK(:) = Xk(:) - Xn(:) - 0.5_RKIND*dt*(fk(:) + fn(:))

         rmsA = 0._RKIND
         rmsR = 0._RKIND
         DO i=1, nX
            rmsA = rmsA + rK(i)**2._RKIND
            rmsR = rmsR + ( rK(i) / (Xk(i)+eps) )**2._RKIND
         END DO
         rmsA = SQRT(rmsA/REAL(nX, KIND=RKIND))
         rmsR = SQRT(rmsR/REAL(nX, KIND=RKIND))

         l1   = k .GT. itMax
         l2   = rmsA .LE. atol
         l3   = rmsR .LE. rtol
         IF (l1 .OR. l2 .OR. l3) EXIT

         CALL TTP_GETJ(imyo, nX, nG, Xk, Xg, JAC, Ksac)
         JAC = Im - 0.5_RKIND*dt*JAC
         JAC = MAT_INV(JAC, nX)
         rK  = MATMUL(JAC, rK)
         Xk  = Xk - rK
      END DO
      Xn(:) = Xk(:)
      CALL TTP_UPDATEG(imyo, dt, nX, nG, Xn, Xg)
      CALL TTP_GETF(imyo, nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR)

      IF (.NOT.l2 .AND. .NOT.l3) IPAR(2) = IPAR(2) + 1

      RETURN
      END SUBROUTINE TTP_INTEGCN2
!-----------------------------------------------------------------------
      SUBROUTINE TTP_GETJ(i, nX, nG, X, Xg, JAC, Ksac)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: i, nX, nG
      REAL(KIND=RKIND), INTENT(IN) :: X(nX), Xg(nG), Ksac
      REAL(KIND=RKIND), INTENT(OUT) :: JAC(nX,nX)

      INCLUDE "PARAMS_TTP.f"

      REAL(KIND=RKIND) :: RT, a, b, c, tau, sq5, e1, e2, e3, e4, n1, n2,
     2   d1, d2, d3

!     Local copies of state variables
      V     = X(1)
      K_i   = X(2)
      Na_i  = X(3)
      Ca_i  = X(4)
      Ca_ss = X(5)
      Ca_sr = X(6)
      R_bar = X(7)

!     Local copies of gating variables
      xr1   = Xg(1)
      xr2   = Xg(2)
      xs    = Xg(3)
      m     = Xg(4)
      h     = Xg(5)
      j     = Xg(6)
      d     = Xg(7)
      f     = Xg(8)
      f2    = Xg(9)
      fcass = Xg(10)
      s     = Xg(11)
      r     = Xg(12)

      RT    = Rc * Tc / Fc
      E_K   = RT * LOG(K_o/K_i)
      E_Na  = RT * LOG(Na_o/Na_i)
      E_Ca  = 0.5_RKIND * RT * LOG(Ca_o/Ca_i)
      E_Ks  = RT * LOG( (K_o + p_KNa*Na_o)/(K_i + p_KNa*Na_i) )

      E_K_Ki   = -RT / K_i
      E_Na_Nai = -RT / Na_i
      E_Ca_Cai = -RT / Ca_i / 2._RKIND
      E_Ks_Ki  = -RT / (K_i + p_KNa*Na_i)
      E_Ks_Nai = p_KNa * E_Ks_Ki

!     I_Na: Fast sodium current
      I_Na = G_Na * (m**3._RKIND) * h * j * (V - E_Na)
      I_Na_V   = G_Na * (m**3._RKIND) * h * j
      I_Na_Nai = I_Na_V * (-E_Na_Nai)

!     I_to: transient outward current
      I_to = G_to(i) * r * s * (V - E_K)
      I_to_V  = G_to(i) * r * s
      I_to_Ki = I_to_V * (-E_K_Ki)

!     I_K1: inward rectifier outward current
      e1   = EXP(0.06_RKIND*(V - E_K - 200._RKIND))
      e2   = EXP(2.E-4_RKIND*(V - E_K + 100._RKIND))
      e3   = EXP(0.1_RKIND*(V - E_K - 10._RKIND))
      e4   = EXP(-0.5_RKIND*(V - E_K))
      a    = 0.1_RKIND/(1._RKIND + e1)
      b    = (3._RKIND*e2 + e3) / (1._RKIND + e4)
      tau  = a / (a + b)
      sq5  = SQRT(K_o/5.4_RKIND)
      n1   = -6.E-3_RKIND*e1/(1._RKIND + e1)**2._RKIND
      n2   = (6.E-4_RKIND*e2 + 0.1_RKIND*e3 + 0.5_RKIND*b*e4)
     2     / (1._RKIND + e4)
      n1   = (a + b)*n1 - a*(n1 + n2)
      d1   = (a + b)**2._RKIND
      I_K1 = G_K1 * sq5 * tau * (V - E_K)
      I_K1_V  = G_K1 * sq5 * (tau + (V - E_K)*n1/d1)
      I_K1_Ki = I_K1_V * (-E_K_Ki)

!     I_Kr: rapid delayed rectifier current
      I_Kr = G_Kr * sq5 * xr1 * xr2 * (V - E_K)
      I_Kr_V   = G_Kr * sq5 * xr1 * xr2
      I_Kr_Ki  = I_Kr_V * (-E_K_Ki)

!     I_Ks: slow delayed rectifier current
      I_Ks = G_Ks(i) * (xs**2._RKIND) * (V - E_Ks)
      I_Ks_V   = G_Ks(i) * (xs**2._RKIND)
      I_Ks_Ki  = I_Ks_V * (-E_Ks_Ki)
      I_Ks_Nai = I_Ks_V * (-E_Ks_Nai)

!     I_CaL: L-type Ca current
      a     = 2._RKIND*(V-15._RKIND)/RT
      b     = (0.25_RKIND*Ca_ss*EXP(a) - Ca_o) / (EXP(a)-1._RKIND)
      c     = G_CaL * d * f * f2 * fcass * (2._RKIND*a*Fc)
      n1    = (EXP(a)/RT) / (EXP(a)-1._RKIND)
      I_CaL = c * b
      I_CaL_V = I_CaL/(V-15._RKIND) + n1*(c*0.5_RKIND*Ca_ss
     2        - 2._RKIND*I_CaL)
      I_CaL_Cass  = c * 0.25_RKIND * n1 * RT

!     I_NaCa: Na-Ca exchanger current
      e1     = EXP(gamma*V/RT)
      e2     = EXP((gamma-1._RKIND)*V/RT)
      n1     = e1*(Na_i**3._RKIND)*Ca_o - e2*(Na_o**3._RKIND)*Ca_i*alpha
      d1     = K_mNai**3._RKIND + Na_o**3._RKIND
      d2     = K_mCa + Ca_o
      d3     = 1._RKIND + K_sat*e2
      c      = 1._RKIND/(d1*d2*d3)
      I_NaCa = K_NaCa * n1 * c

      n1     = K_NaCa * c * ( e1*(Na_i**3._RKIND)*Ca_o*(gamma/RT) -
     2    e2*(Na_o**3._RKIND)*Ca_i*alpha*((gamma-1._RKIND)/RT) )
      n2     = I_NaCa*K_sat*((gamma-1._RKIND)/RT)*e2/d3
      I_NaCa_V   = n1 - n2
      I_NaCa_Nai =  K_NaCa * e1 * (3._RKIND*Na_i**2._RKIND) * Ca_o * c
      I_NaCa_Cai = -K_NaCa * e2 * (Na_o**3._RKIND) * alpha * c

!     I_NaK: Na-K pump current
      e1    = EXP(-0.1_RKIND*V/RT)
      e2    = EXP(-V/RT)
      n1    = P_NaK * K_o * Na_i
      d1    = K_o + K_mK
      d2    = Na_i + K_mNa
      d3    = 1._RKIND + 0.1245_RKIND*e1 + 0.0353_RKIND*e2
      I_NaK = n1 / (d1*d2*d3)
      n1    = (0.01245_RKIND*e1 + 0.0353_RKIND*e2)/RT
      I_NaK_V   = I_NaK * n1 / d3
      I_NaK_Nai = I_NaK * K_mNa/(Na_i*d2)

!     I_pCa: plateau Ca current
      d1 = (K_pCa + Ca_i)
      I_pCa = G_pCa * Ca_i / d1
      I_pCa_Cai = G_pCa * K_pCa/(d1*d1)

!     I_pK: plateau K current
      e1 = EXP((25._RKIND-V)/5.98_RKIND)
      I_pK  = G_pK * (V-E_K) / (1._RKIND + e1)
      I_pK_V  = (G_pK + I_pK*e1/5.98_RKIND) / (1._RKIND+e1)
      I_pK_Ki = G_pK * (-E_K_Ki) / (1._RKIND + e1)

!     I_bCa: background Ca current
      I_bCa = G_bCa * (V - E_Ca)
      I_bCa_V = G_bCa
      I_bCa_Cai = G_bca * (-E_Ca_Cai)

!     I_bNa: background Na current
      I_bNa = G_bNa * (X(1) - E_Na)
      I_bNa_V = G_bNa
      I_bNa_Nai = G_bNa * (-E_Na_Nai)

!     I_leak: Sacroplasmic Reticulum Ca leak current
      I_leak = V_leak * (Ca_sr - Ca_i)
      I_leak_Cai  = -V_leak
      I_leak_Casr =  V_leak

!     I_up: Sacroplasmic Reticulum Ca pump current
      d1 = 1._RKIND + (K_up/Ca_i)**2._RKIND
      I_up  = Vmax_up / d1
      I_up_Cai = (I_up / d1) * (2._RKIND*K_up**2._RKIND /Ca_i**3._RKIND)

!     I_rel: Ca induced Ca current (CICR)
      n1 = max_sr - min_sr
      d1 = 1._RKIND + (EC/Ca_sr)**2._RKIND
      k_casr = max_sr - (n1/d1)
      k1     = k1p / k_casr
      n2     = Ca_ss*2._RKIND
      d2     = k3 + k1*n2
      O      = k1 * R_bar * n2 / d2
      I_rel  = V_rel * O * (Ca_sr - Ca_ss)

      k_casr_sr = (n1 / (d1**2._RKIND)) * (2._RKIND*EC**2._RKIND
     2          / Ca_sr**3._RKIND)
      k1_casr   = -k1p * k_casr_sr / (k_casr**2._RKIND)
      O_Cass = 2._RKIND * k3 * O / (Ca_ss * d2)
      O_Casr = k1_casr * n2 * (R_bar - O) / d2
      O_Rbar = k1 * n2 / d2

      I_rel_Cass = V_rel * (O_Cass*(Ca_sr - Ca_ss) - O)
      I_rel_Casr = V_rel * (O_Casr*(Ca_sr - Ca_ss) + O)
      I_rel_Rbar = V_rel * O_Rbar *(Ca_sr - Ca_ss)

!     I_xfer: diffusive Ca current between Ca subspae and cytoplasm
      I_xfer = V_xfer * (Ca_ss - Ca_i)
      I_xfer_Cai  = -V_xfer
      I_xfer_Cass =  V_xfer

!-----------------------------------------------------------------------
!     Compute Jacobian matrix
      JAC = 0._RKIND
      c   = Cm/(V_c*Fc)

!     V
      JAC(1,1)  = -(I_Na_V + I_to_V + I_K1_V + I_Kr_V + I_Ks_V + I_CaL_V
     2   + I_NaCa_V + I_NaK_V + I_pK_V + I_bCa_V + I_bNa_V + Ksac)
      JAC(1,2)  = -(I_to_Ki + I_K1_Ki + I_Kr_Ki + I_Ks_Ki + I_pK_Ki)
      JAC(1,3)  = -(I_Na_Nai + I_Ks_Nai + I_NaCa_Nai + I_NaK_Nai
     2   + I_bNa_Nai)
      JAC(1,4)  = -(I_NaCa_Cai + I_pCa_Cai + I_bCa_Cai)
      JAC(1,5) = -I_CaL_Cass

!     K_i
      JAC(2,1)  = -c * (I_K1_V + I_to_V + I_Kr_V + I_Ks_V + I_pK_V
     2   - 2._RKIND*I_NaK_V )
      JAC(2,2)  = -c * (I_K1_Ki + I_to_Ki + I_Kr_Ki + I_Ks_Ki + I_pK_Ki)
      JAC(2,3)  = -c * (I_Ks_Nai - 2._RKIND*I_NaK_Nai)

!     Na_i
      JAC(3,1)  = -c * (I_Na_V + I_bNa_V +
     2   3._RKIND*(I_NaK_V + I_NaCa_V))
      JAC(3,3)  = -c * (I_Na_Nai + I_bNa_Nai + 3._RKIND*(I_NaK_Nai +
     2   I_NaCa_Nai))
      JAC(3,4)  = -c * (3._RKIND*I_NaCa_Cai)

!     Ca_i
      n1 = (I_leak - I_up)*V_sr/V_c + I_xfer
     2   - 0.5_RKIND*c*(I_bCa + I_pCa - 2._RKIND*I_NaCa)
      n2 = (I_leak_Cai - I_up_Cai)*V_sr/V_c + I_xfer_Cai
     2   - 0.5_RKIND*c*(I_bCa_Cai + I_pCa_Cai - 2._RKIND*I_NaCa_Cai)
      d1 = 1._RKIND + K_bufc*Buf_c/(Ca_i + K_bufc)**2._RKIND
      d2 = 2._RKIND*K_bufc*Buf_c/(Ca_i + K_bufc)**3._RKIND
      JAC(4,1)  = -c * (I_bCa_V - 2._RKIND*I_NaCa_V) / 2._RKIND / d1
      JAC(4,3)  = c * I_NaCa_Nai / d1
      JAC(4,4)  = (n2 + n1*d2/d1) / d1
      JAC(4,5) = I_xfer_Cass / d1
      JAC(4,6) = (I_leak_Casr*V_sr/V_c) / d1

!     Ca_ss
      a  = Cm/(2._RKIND*Fc*V_ss)
      b  = V_sr / V_ss
      c  = V_c / V_ss
      n1 = -a*I_CaL + b*I_rel - c*I_xfer
      n2 = -a*I_CaL_Cass + b*I_rel_Cass - c*I_xfer_Cass
      d1 = 1._RKIND + K_bufss*Buf_ss/(Ca_ss + K_bufss)**2._RKIND
      d2 = 2._RKIND*K_bufss*Buf_ss/(Ca_ss + K_bufss)**3._RKIND
      JAC(5,1)  = -a * I_CaL_V / d1
      JAC(5,4)  = -c * I_xfer_Cai / d1
      JAC(5,5) = (n2 + n1*d2/d1) / d1
      JAC(5,6) = b * I_rel_Casr / d1
      JAC(5,7) = b * I_rel_Rbar / d1

!     Ca_sr
      n1 = I_up - I_leak - I_rel
      n2 = -(I_leak_Casr + I_rel_Casr)
      d1 = 1._RKIND + K_bufsr*Buf_sr/(Ca_sr + K_bufsr)**2._RKIND
      d2 = 2._RKIND*K_bufsr*Buf_sr/(Ca_sr + K_bufsr)**3._RKIND
      JAC(6,4) = (I_up_Cai - I_leak_Cai) / d1
      JAC(6,5) = -I_rel_Cass / d1
      JAC(6,6) = (n2 + n1*d2/d1) / d1
      JAC(6,7) = -I_rel_Rbar / d1

!     Rbar: ryanodine receptor
      k2 = k2p * k_casr
      JAC(7,5) = -k2 * R_bar
      JAC(7,6) = -(k2p * k_casr_sr) * Ca_ss * R_bar
      JAC(7,7) = -(k2*Ca_ss + k4)

      RETURN
      END SUBROUTINE TTP_GETJ
!-----------------------------------------------------------------------
!     Compute activation force for electromechanics based on active
!     stress model
      SUBROUTINE TTP_ACTVSTRS(c_Ca, dt, Tact, epsX)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: c_Ca, dt
      REAL(KIND=RKIND), INTENT(OUT) :: epsX
      REAL(KIND=RKIND), INTENT(INOUT) :: Tact

      INCLUDE "PARAMS_TTP.f"

      REAL(KIND=RKIND) :: nr

      epsX = EXP(-EXP(-xi_T*(c_Ca - Ca_crit)))
      epsX = eps_0 + (eps_i - eps_0)*epsX
      nr   = Tact + epsX*dt*eta_T*(c_Ca - Ca_rest)
      Tact = nr / (1._RKIND + epsX*dt)

      RETURN
      END SUBROUTINE TTP_ACTVSTRS
!-----------------------------------------------------------------------
!     Compute macroscopic fiber strain based on sacromere force-length
!     relationship and calcium concentration
      SUBROUTINE TTP_ACTVSTRN(c_Ca, I4f, dt, gf)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: c_Ca, I4f, dt
      REAL(KIND=RKIND), INTENT(INOUT) :: gf

      INCLUDE "PARAMS_TTP.f"

      REAL(KIND=RKIND) :: SL, Fa, rtmp

!     fiber length
      SL = I4f * SL0

!     Sacromere force-length relationship
      IF (SL.GE.SLmin .AND. SL.LE.SLmax) THEN
         SL = 0.5_RKIND*f0 + fc1*COS(SL) + fs1*SIN(SL) +
     2      fc2*COS(2._RKIND*SL) + fs2*SIN(2._RKIND*SL)  +
     3      fc3*COS(3._RKIND*SL) + fs3*SIN(3._RKIND*SL)
      ELSE
         SL = 0._RKIND
      END IF

!     Active force
      Fa   = alFa * (c_Ca-c_Ca0)*(c_Ca-c_Ca0) * SL

      rtmp = 2._RKIND*I4f*(1._RKIND/(1._RKIND+gf)**3._RKIND - 1._RKIND)
      gf = gf + dt*(Fa + rtmp)/(mu_Ca * c_Ca * c_Ca)

      RETURN
      END SUBROUTINE TTP_ACTVSTRN
!-----------------------------------------------------------------------
      END MODULE TTPMOD
!#######################################################################

ext: f

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


module chimera_parser_module

  ! we read in the number of variables and their order and use this to
  ! map them into the model_state array.  We ignore anything other than
  ! density, temperature, pressure and composition.
  !
  ! composition is assumed to be in terms of mass fractions

  use amrex_paralleldescriptor_module, only: parallel_IOProcessor=>amrex_pd_ioprocessor
  use network
  use amrex_fort_module, only: rt => amrex_real
  use hdf5_read_write
  use hdf5

  implicit none

  ! array dimensions
  integer, save :: nx_chim, ny_chim, nz_chim, nnc_chim, npart_chim, npart_shell_chim

  ! radial index bounds
  integer, save :: imin_chim, imax_chim

  ! theta index bounds
  integer, save :: jmin_chim, jmax_chim

  ! phi index bounds
  integer, save :: kmin_chim, kmax_chim

  ! grid coordinates
  real (rt), allocatable, save :: x_e_chim(:), dx_e_chim(:)
  real (rt), allocatable, save :: y_e_chim(:), dy_e_chim(:)
  real (rt), allocatable, save :: z_e_chim(:), dz_e_chim(:)
  real (rt), allocatable, save :: x_c_chim(:), dx_c_chim(:)
  real (rt), allocatable, save :: y_c_chim(:), dy_c_chim(:)
  real (rt), allocatable, save :: z_c_chim(:), dz_c_chim(:)

  real (rt), allocatable, save :: volx_e_chim(:), dvolx_e_chim(:)
  real (rt), allocatable, save :: voly_e_chim(:), dvoly_e_chim(:)
  real (rt), allocatable, save :: volz_e_chim(:), dvolz_e_chim(:)
  real (rt), allocatable, save :: volx_c_chim(:), dvolx_c_chim(:)
  real (rt), allocatable, save :: voly_c_chim(:), dvoly_c_chim(:)
  real (rt), allocatable, save :: volz_c_chim(:), dvolz_c_chim(:)

  real (rt), allocatable, save :: domega_chim(:,:)
  real (rt), allocatable, save :: dvol_e_chim(:,:,:)
  real (rt), allocatable, save :: dmass_e_chim(:,:,:)
  real (rt), save :: omega_chim

  ! time data
  real (rt), save :: time_chim, t_bounce

  ! fluid velocities
  real (rt), allocatable, save :: u_c_chim(:,:,:)
  real (rt), allocatable, save :: v_c_chim(:,:,:)
  real (rt), allocatable, save :: w_c_chim(:,:,:)

  ! thermo data
  real (rt), allocatable, save :: rho_c_chim(:,:,:)
  real (rt), allocatable, save :: t_c_chim(:,:,:)
  real (rt), allocatable, save :: ye_c_chim(:,:,:)
  real (rt), allocatable, save :: p_c_chim(:,:,:)
  real (rt), allocatable, save :: ei_c_chim(:,:,:)
  real (rt), allocatable, save :: et_c_chim(:,:,:)
  real (rt), allocatable, save :: eb_c_chim(:,:,:)
  real (rt), allocatable, save :: em_c_chim(:,:,:)
  real (rt), allocatable, save :: s_c_chim(:,:,:)

  real (rt), allocatable, save :: rhobar_c_chim(:)
  real (rt), allocatable, save :: tbar_c_chim(:)

  ! gravity data
  real (rt), allocatable, save :: gravx_c_chim(:,:,:)
  real (rt), allocatable, save :: gravy_c_chim(:,:,:)
  real (rt), allocatable, save :: gravz_c_chim(:,:,:)
  real (rt), allocatable, save :: gravx_c_avg_chim(:)

  ! nuclei data
  real (rt), allocatable, save :: a_nuc_chim(:)
  real (rt), allocatable, save :: z_nuc_chim(:)
  real (rt), allocatable, save :: mex_nuc_chim(:)
  real (rt), allocatable, save :: be_nuc_chim(:)
  character (5), allocatable, save :: name_nuc_chim(:)

  ! composition data
  real (rt), allocatable, save :: xn_read(:,:,:,:)
  real (rt), allocatable, save :: xn_c_chim(:,:,:,:)
  integer,     allocatable, save :: nse_c_chim(:,:,:)

  ! representative heavy/auxiliary data
  real (rt), allocatable, save :: a_aux_c_chim(:,:,:)
  real (rt), allocatable, save :: z_aux_c_chim(:,:,:)
  real (rt), allocatable, save :: be_aux_c_chim(:,:,:)

  ! particle data
  real (rt), allocatable, save :: px_chim(:)
  real (rt), allocatable, save :: py_chim(:)
  real (rt), allocatable, save :: pz_chim(:)

  ! chimera_initialized will be .true. once the model is read in and the
  ! model data arrays are initialized and filled
  logical, save :: chimera_initialized = .false.

  ! HDF5 variables
  integer(hid_t) :: file_id
  integer(hid_t) :: group_id
  integer(hid_t) :: dataset_id
  integer(hid_t) :: dataspace_id
  integer(hid_t) :: plist_id

  integer(hsize_t) :: datasize1d(1)
  integer(hsize_t) :: datasize2d(2)
  integer(hsize_t) :: datasize3d(3)
  integer(hsize_t) :: datasize4d(4)
  integer(hsize_t) :: datasize5d(5)

  integer(hsize_t) :: slab_offset2d(2)
  integer(hsize_t) :: slab_offset3d(3)
  integer(hsize_t) :: slab_offset4d(4)
  integer(hsize_t) :: slab_offset5d(5)

  public :: open_chimera_file, read_chimera_file, close_chimera_file

contains

  subroutine open_chimera_file(fname)

    use amrex_constants_module
    use amrex_error_module

    character (*), intent(in) :: fname

    ! local variables
    integer :: i_read(1)

    integer :: ierr

    call h5open_f(ierr)
    call h5fopen_f(fname, H5F_ACC_RDONLY_F, file_id, ierr)

    if (ierr /= 0) then
      if (parallel_IOProcessor()) print *,'Could not open file: ',fname
      call bl_error('Aborting now -- please supply filename')
    end if

    call h5gopen_f(file_id, '/mesh', group_id, ierr)
    if (ierr /= 0) call bl_error('Could not open /mesh group')

    datasize1d(1) = 1
    call read_1d_slab('nz_hyperslabs', i_read, group_id, datasize1d)

    call h5gclose_f(group_id, ierr)

    return

  end subroutine open_chimera_file

  subroutine read_chimera_file

    use actual_network, only: nspec, aion, zion
    use amrex_constants_module
    use amrex_error_module
    use fundamental_constants_module
    use probdata_module

    ! local variables
    real (rt) :: d_read(1)

    integer :: array_dimensions(3)
    integer :: index_bounds(2)
    integer :: eos_dimensions(2)
    integer :: part_dimensions(2)

    integer :: i, j, k, ierr
    integer :: itmp(nspec)
    integer :: net_to_castro(nspec)
    integer, allocatable :: net_in_castro(:)

    real (rt) :: dx, dvol

    ! open mesh group
    call h5gopen_f(file_id, '/mesh', group_id, ierr)
    if (ierr /= 0) call bl_error('Could not open /mesh group')

    ! read array dimensions
    datasize1d(1) = 3
    call read_1d_slab('array_dimensions', array_dimensions, group_id, datasize1d)
    nx_chim = array_dimensions(1)
    ny_chim = array_dimensions(2)
    nz_chim = array_dimensions(3)

    ! read radial index bounds
    datasize1d(1) = 2
    call read_1d_slab('radial_index_bound', index_bounds, group_id, datasize1d)
    imin_chim = index_bounds(1)
    imax_chim = index_bounds(2)

    ! read theta index bounds
    call read_1d_slab('theta_index_bound', index_bounds, group_id, datasize1d)
    jmin_chim = index_bounds(1)
    jmax_chim = index_bounds(2)
    if ( jmin_chim /= 1 .or. jmax_chim /= ny_chim ) then
       if (parallel_IOProcessor()) print *,'Error: inner theta bounds /= (1,ny_chim), (',jmin_chim,',',jmax_chim,')'
       call bl_error('Aborting now -- check theta bounds in model file')
     end if

    ! read phi index bounds
    call read_1d_slab('phi_index_bound', index_bounds, group_id, datasize1d)
    kmin_chim = index_bounds(1)
    kmax_chim = index_bounds(2)
    if ( kmin_chim /= 1 .or. kmax_chim /= nz_chim ) then
       if (parallel_IOProcessor()) print *,'Error: inner phi bounds /= (1,nz_chim), (',kmin_chim,',',kmax_chim,')'
       call bl_error('Aborting now -- check phi bounds in model file')
     end if

    ! allocate grid variables
    allocate (x_e_chim(nx_chim+1))
    allocate (y_e_chim(-5:ny_chim+7))
    allocate (z_e_chim(-5:nz_chim+7))
    allocate (x_c_chim(nx_chim))
    allocate (y_c_chim(-5:ny_chim+6))
    allocate (z_c_chim(-5:nz_chim+6))
    allocate (dx_e_chim(nx_chim))
    allocate (dy_e_chim(-5:ny_chim+6))
    allocate (dz_e_chim(-5:nz_chim+6))
    allocate (dx_c_chim(nx_chim))
    allocate (dy_c_chim(-5:ny_chim+6))
    allocate (dz_c_chim(-5:nz_chim+6))

    ! read zone edge coordinaets
    datasize1d(1) = nx_chim+1
    call read_1d_slab('x_ef', x_e_chim, group_id, datasize1d)
    datasize1d(1) = ny_chim+1
    call read_1d_slab('y_ef', y_e_chim(1:ny_chim+1), group_id, datasize1d)
    datasize1d(1) = nz_chim+1
    call read_1d_slab('z_ef', z_e_chim(1:nz_chim+1), group_id, datasize1d)

    ! read zone center coordinates
    datasize1d(1) = nx_chim
    call read_1d_slab('x_cf', x_c_chim, group_id, datasize1d)
    datasize1d(1) = ny_chim
    call read_1d_slab('y_cf', y_c_chim(1:ny_chim), group_id, datasize1d)
    datasize1d(1) = nz_chim
    call read_1d_slab('z_cf', z_c_chim(1:nz_chim), group_id, datasize1d)

    ! read zone widths
    datasize1d(1) = nx_chim
    call read_1d_slab('dx_cf', dx_c_chim, group_id, datasize1d)
    datasize1d(1) = ny_chim
    call read_1d_slab('dy_cf', dy_c_chim(1:ny_chim), group_id, datasize1d)
    datasize1d(1) = nz_chim
    call read_1d_slab('dz_cf', dz_c_chim(1:nz_chim), group_id, datasize1d)

    ! fill grid boundaries
    call grid_bc(y_e_chim, y_c_chim, dy_c_chim, 1, ny_chim, 0, 0)
    call grid_bc(z_e_chim, z_c_chim, dz_c_chim, 1, nz_chim, 3, 3)
    dx_e_chim(1:nx_chim) = x_e_chim(2:nx_chim+1) - x_e_chim(1:nx_chim)
    dy_e_chim(-5:ny_chim+6) = y_e_chim(-4:ny_chim+7) - y_e_chim(-5:ny_chim+6)
    dz_e_chim(-5:nz_chim+6) = z_e_chim(-4:nz_chim+7) - z_e_chim(-5:nz_chim+6)

    ! read times
    datasize1d(1) = 1
    call read_1d_slab('time', d_read, group_id, datasize1d)
    time_chim = d_read(1)
    call read_1d_slab('t_bounce', d_read, group_id, datasize1d)
    t_bounce = d_read(1)

    ! solid angles and volumes
    allocate (volx_e_chim(nx_chim+1))
    allocate (voly_e_chim(ny_chim+1))
    allocate (volz_e_chim(nz_chim+1))
    allocate (volx_c_chim(nx_chim))
    allocate (voly_c_chim(ny_chim))
    allocate (volz_c_chim(nz_chim))
    allocate (dvolx_e_chim(nx_chim))
    allocate (dvoly_e_chim(ny_chim))
    allocate (dvolz_e_chim(nz_chim))
    allocate (dvolx_c_chim(nx_chim))
    allocate (dvoly_c_chim(ny_chim))
    allocate (dvolz_c_chim(nz_chim))
    allocate (domega_chim(ny_chim,nz_chim))
    allocate (dvol_e_chim(nx_chim,ny_chim,nz_chim))
    allocate (dmass_e_chim(nx_chim,ny_chim,nz_chim))

    datasize2d = (/ ny_chim, nz_chim /)
    call read_2d_slab('d_omega', domega_chim, group_id, datasize2d)

    volx_e_chim(1) = third * x_e_chim(1)**3
    do i = 1, nx_chim
      dx               = dx_e_chim(i)
      dvolx_e_chim(i)  = dx * ( x_e_chim(i) * x_e_chim(i+1) + dx * dx * third )
      volx_e_chim(i+1) = volx_e_chim(i) + dvolx_e_chim(i)
    end do
    volx_c_chim(1:nx_chim) = volx_e_chim(1:nx_chim) + half*dvolx_e_chim(1:nx_chim)
    dvolx_c_chim(1:nx_chim-1) = volx_c_chim(2:nx_chim) - volx_c_chim(1:nx_chim-1)
    dvolx_c_chim(nx_chim) = dvolx_c_chim(nx_chim-1)

    voly_e_chim(1) = one - cos( y_e_chim(1) )
    do j = 1, ny_chim
      dvoly_e_chim(j)  = cos( y_e_chim(j) ) - cos( y_e_chim(j+1) )
      voly_e_chim(j+1) = voly_e_chim(j) + dvoly_e_chim(j)
    end do
    voly_c_chim(1:ny_chim) = voly_e_chim(1:ny_chim) + half*dvoly_e_chim(1:ny_chim)
    dvoly_c_chim(1:ny_chim-1) = voly_c_chim(2:ny_chim) - voly_c_chim(1:ny_chim-1)
    dvoly_c_chim(ny_chim) = dvoly_c_chim(ny_chim-1)

    volz_e_chim(1) = z_e_chim(1)
    do k = 1, nz_chim
      dvolz_e_chim(k)  = dz_e_chim(k)
      volz_e_chim(k+1) = volz_e_chim(k) + dvolz_e_chim(k)
    end do
    volz_c_chim(1:nz_chim) = volz_e_chim(1:nz_chim) + half*dvolz_e_chim(1:nz_chim)
    dvolz_c_chim(1:nz_chim-1) = volz_c_chim(2:nz_chim) - volz_c_chim(1:nz_chim-1)
    dvolz_c_chim(nz_chim) = dvolz_c_chim(nz_chim-1)

    omega_chim = sum( domega_chim(:,:) )
    do i = 1, nx_chim
      dvol_e_chim(i,:,:) = dvolx_e_chim(i) * domega_chim(:,:)
    end do

    ! close mesh group
    call h5gclose_f(group_id, ierr)

    ! open fluid group
    call h5gopen_f(file_id, '/fluid', group_id, ierr)
    if (ierr /= 0) call bl_error('Could not open /fluid group')

    ! allocate fluid variables
    allocate (u_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (v_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (w_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (rho_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (t_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (ye_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (p_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (ei_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (et_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (eb_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (em_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (s_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (gravx_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (gravy_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (gravz_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (gravx_c_avg_chim(nx_chim))
    allocate (rhobar_c_chim(nx_chim))
    allocate (tbar_c_chim(nx_chim))

    ! read velocity variables
    datasize3d = (/ nx_chim, ny_chim, nz_chim /)
    slab_offset3d = (/ 0, 0, 0 /)
    call read_ray_hyperslab('u_c',     u_c_chim, group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('v_c',     v_c_chim, group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('w_c',     w_c_chim, group_id, datasize3d, slab_offset3d)

    ! read thermo variables
    call read_ray_hyperslab('rho_c',   rho_c_chim, group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('t_c',     t_c_chim,   group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('ye_c',    ye_c_chim,  group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('press',   p_c_chim,   group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('e_int',   ei_c_chim,  group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('entropy', s_c_chim,   group_id, datasize3d, slab_offset3d)

    ! convert units to erg/K
    s_c_chim(:,:,:) = s_c_chim(:,:,:) * k_B * n_A

    dmass_e_chim(:,:,:) = rho_c_chim(:,:,:) * dvol_e_chim(:,:,:)

    ! read gravity variables
    call read_ray_hyperslab('grav_x_c', gravx_c_chim, group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('grav_y_c', gravy_c_chim, group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('grav_z_c', gravz_c_chim, group_id, datasize3d, slab_offset3d)
    do i = imax_chim+1, nx_chim
      gravx_c_chim(i,:,:) = gravx_c_chim(imax_chim,:,:) * (volx_c_chim(imax_chim)/volx_c_chim(i))**two3rd
    end do

    do i = 1, nx_chim
      rhobar_c_chim(i) = sum( rho_c_chim(i,:,:) * domega_chim(:,:) ) / omega_chim
      tbar_c_chim(i) = sum( t_c_chim(i,:,:) * domega_chim(:,:) ) / omega_chim
      gravx_c_avg_chim(i) = sum( gravx_c_chim(i,:,:) * domega_chim(:,:) ) / omega_chim
    end do

    ! read number of nuclei
    datasize1d(1) = 2
    call read_1d_slab('dimeos', eos_dimensions, group_id, datasize1d)
    nnc_chim = eos_dimensions(1)

    ! close fluid group
    call h5gclose_f(group_id, ierr)

    ! open abundance group
    call h5gopen_f(file_id, '/abundance', group_id, ierr)
    if (ierr /= 0) call bl_error('Could not open /abundance group')

    ! allocate composition variables
    allocate (a_nuc_chim(nnc_chim-1))
    allocate (z_nuc_chim(nnc_chim-1))
    allocate (mex_nuc_chim(nnc_chim-1))
    allocate (be_nuc_chim(nnc_chim-1))
    allocate (name_nuc_chim(nnc_chim))
    allocate (xn_read(nnc_chim,nx_chim,ny_chim,nz_chim))
    allocate (xn_c_chim(nspec,nx_chim,ny_chim,nz_chim))
    allocate (nse_c_chim(nx_chim+1,ny_chim,nz_chim))
    allocate (a_aux_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (z_aux_c_chim(nx_chim,ny_chim,nz_chim))
    allocate (be_aux_c_chim(nx_chim,ny_chim,nz_chim))

    ! read nuclei values
    datasize1d(1) = nnc_chim-1
    call read_1d_slab('a_nuc',    a_nuc_chim,   group_id, datasize1d)
    call read_1d_slab('z_nuc',    z_nuc_chim,   group_id, datasize1d)
    call read_1d_slab('m_ex_nuc', mex_nuc_chim, group_id, datasize1d)

    ! binding energy
    be_nuc_chim(1:nnc_chim-1) = z_nuc_chim(1:nnc_chim-1)*7.2889705_rt + &
    &                           (a_nuc_chim(1:nnc_chim-1)-z_nuc_chim(1:nnc_chim-1))*8.0713171_rt - &
    &                           mex_nuc_chim(1:nnc_chim-1)

    ! read nuclei names
    datasize1d(1) = nnc_chim
    call read_1d_slab('a_name', name_nuc_chim, 5, group_id, datasize1d)

    ! read mass fractions
    datasize4d = (/ nnc_chim, nx_chim, ny_chim, nz_chim /)
    slab_offset4d = (/ 0, 0, 0, 0 /)
    call read_ray_hyperslab('xn_c', xn_read, group_id, datasize4d, slab_offset4d)

    ! create lookup tables for isotopes in castro net
    k = 0
    net_to_castro(:) = nnc_chim
    do i = 1, nspec
      do j = 1, nnc_chim-1
        if ( nint(aion(i)) == nint(a_nuc_chim(j)) .and. nint(zion(i)) == nint(z_nuc_chim(j)) ) then
          k = k + 1
          net_to_castro(i) = j
          itmp(k) = i
          exit
        end if
      end do
      if ( j > nnc_chim-1 ) then
        if (parallel_IOProcessor()) then
          write(*,'(2(a,i3),a)') ' could not find isotope (',nint(zion(i)),',',nint(aion(i)),') in CHIMERA net'
        end if
      end if
    end do

    if ( k > 0 ) then
      allocate (net_in_castro(k))
      net_in_castro(:) = itmp(1:k)
    else
      call bl_error("no species in chimera net in castro net")
    end if

    xn_c_chim(net_in_castro,:,:,:) = xn_read(net_to_castro(net_in_castro),:,:,:)
    if ( nspec > nspec_evolve ) then
      xn_c_chim(nspec,:,:,:) = xn_read(nnc_chim,:,:,:)
    end if

    ! read nse state
    datasize3d = (/ nx_chim+1, ny_chim, nz_chim /)
    slab_offset3d = (/ 0, 0, 0 /)
    call read_ray_hyperslab('nse_c', nse_c_chim, group_id, datasize3d, slab_offset3d)

    ! read representative heavy nucleus/auxiliary nucleus properties
    datasize3d = (/ nx_chim, ny_chim, nz_chim /)
    slab_offset3d = (/ 0, 0, 0 /)
    call read_ray_hyperslab('a_nuc_rep_c',  a_aux_c_chim,  group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('z_nuc_rep_c',  z_aux_c_chim,  group_id, datasize3d, slab_offset3d)
    call read_ray_hyperslab('be_nuc_rep_c', be_aux_c_chim, group_id, datasize3d, slab_offset3d)

    ! binding energy
    eb_c_chim(:,:,:) = zero
    do i = 1, nnc_chim-1
      eb_c_chim(:,:,:) = eb_c_chim(:,:,:) - xn_read(i,:,:,:)*be_nuc_chim(i)/a_nuc_chim(i)
    end do
    eb_c_chim(:,:,:) = eb_c_chim(:,:,:) - xn_read(nnc_chim,:,:,:)*be_aux_c_chim(:,:,:)/a_aux_c_chim(:,:,:)

    ! convert from MeV/c^2 to erg/g
    eb_c_chim(:,:,:) = eb_c_chim(:,:,:) * MeV2eV * ev2erg * n_A

    ! rest mass energy (including electrons) relative to free neutrons
    em_c_chim(:,:,:) = - ( m_n - m_p - m_e ) * ye_c_chim(:,:,:) * c_light * c_light * n_A + eb_c_chim(:,:,:)

    ! remove rest mass energy and constant offset from internal energy
    et_c_chim(:,:,:) = ei_c_chim(:,:,:) - em_c_chim(:,:,:) - 8.9_rt * MeV2eV * ev2erg * n_A

    ! close abundance group
    call h5gclose_f(group_id, ierr)

    if ( do_particles ) then
      ! open particle group
      call h5gopen_f(file_id, '/particle', group_id, ierr)
      if (ierr /= 0) then
        call bl_error('Could not open /particle group')
      else
        ! read particle dimensions
        datasize1d(1) = 2
        call read_1d_slab('partdim', part_dimensions, group_id, datasize1d)
        npart_chim = part_dimensions(1)
        npart_shell_chim = part_dimensions(2)

        ! allocate particle variables
        allocate (px_chim(npart_chim))
        allocate (py_chim(npart_chim))
        allocate (pz_chim(npart_chim))

        ! read particle variables
        datasize1d(1) = npart_chim
        call read_1d_slab('px', px_chim, group_id, datasize1d)
        call read_1d_slab('py', py_chim, group_id, datasize1d)
        call read_1d_slab('pz', pz_chim, group_id, datasize1d)

        ! close particle group
        call h5gclose_f(group_id, ierr)
      end if
    end if

    ! close hdf5 file
    call h5fclose_f(file_id, ierr)
    call h5close_f(ierr)

    chimera_initialized = .true.

    return

  end subroutine read_chimera_file

  subroutine close_chimera_file

    use probdata_module

    if (chimera_initialized) then
       deallocate(x_e_chim)
       deallocate(y_e_chim)
       deallocate(z_e_chim)
       deallocate(x_c_chim)
       deallocate(y_c_chim)
       deallocate(z_c_chim)
       deallocate(dx_e_chim)
       deallocate(dy_e_chim)
       deallocate(dz_e_chim)
       deallocate(dx_c_chim)
       deallocate(dy_c_chim)
       deallocate(dz_c_chim)
       deallocate(volx_e_chim)
       deallocate(voly_e_chim)
       deallocate(volz_e_chim)
       deallocate(volx_c_chim)
       deallocate(voly_c_chim)
       deallocate(volz_c_chim)
       deallocate(dvolx_e_chim)
       deallocate(dvoly_e_chim)
       deallocate(dvolz_e_chim)
       deallocate(dvolx_c_chim)
       deallocate(dvoly_c_chim)
       deallocate(dvolz_c_chim)
       deallocate(domega_chim)
       deallocate(dvol_e_chim)
       deallocate(dmass_e_chim)
       deallocate(u_c_chim)
       deallocate(v_c_chim)
       deallocate(w_c_chim)
       deallocate(rho_c_chim)
       deallocate(t_c_chim)
       deallocate(ye_c_chim)
       deallocate(p_c_chim)
       deallocate(ei_c_chim)
       deallocate(et_c_chim)
       deallocate(eb_c_chim)
       deallocate(em_c_chim)
       deallocate(s_c_chim)
       deallocate(gravx_c_chim)
       deallocate(gravy_c_chim)
       deallocate(gravz_c_chim)
       deallocate(gravx_c_avg_chim)
       deallocate(rhobar_c_chim)
       deallocate(tbar_c_chim)
       deallocate(a_nuc_chim)
       deallocate(z_nuc_chim)
       deallocate(mex_nuc_chim)
       deallocate(be_nuc_chim)
       deallocate(name_nuc_chim)
       deallocate(xn_read)
       deallocate(xn_c_chim)
       deallocate(nse_c_chim)
       deallocate(a_aux_c_chim)
       deallocate(z_aux_c_chim)
       deallocate(be_aux_c_chim)
       if ( do_particles ) then
         deallocate(px_chim)
         deallocate(py_chim)
         deallocate(pz_chim)
       end if
       chimera_initialized = .false.
    endif

    return

  end subroutine close_chimera_file

  subroutine grid_bc( xag, xcg, dxg, nmin, nmax, nleft, nright )

    use amrex_constants_module
    use amrex_error_module

    !  boundary condition flags : nleft, nright
    !    = 0 : reflecting
    !    = 3 : periodic (eg, u(nmin-1) = u(nmax))

    ! input variables
    integer, intent(in) :: nmin         ! minimum coordinate index
    integer, intent(in) :: nmax         ! maximum coordinate index
    integer, intent(in) :: nleft        ! left boundary type
    integer, intent(in) :: nright       ! right boundary type

    ! output variables
    real (rt), intent(inout) :: xag(nmin-6:nmax+7)    ! grid edges
    real (rt), intent(inout) :: xcg(nmin-6:nmax+6)    ! grid centers
    real (rt), intent(inout) :: dxg(nmin-6:nmax+6)    ! grid widths

    ! local variables
    integer :: n        ! zone index

    ! left (inner) boundary
    select case (nleft)

    case (0,162)        ! grid symmetric with left (inner) edge

      do n = 1, 6
        dxg(nmin - n)  = dxg(nmin + n - 1)
        xag(nmin - n)  = xag(nmin - n + 1) - dxg(nmin-n)
        xcg(nmin - n)  = xag(nmin - n) + half * dxg(nmin - n)
      end do

    case (3)            ! periodic

      do n = 1, 6
        dxg(nmin - n)  = dxg(nmax + 1 - n)
        xag(nmin - n)  = xag(nmin - n + 1) - dxg(nmin - n)
        xcg(nmin - n)  = xag(nmin - n) + half * dxg(nmin - n)
      end do

    case default

      call bl_error("invalid value for boundary condition, nleft")

    end select

    ! right (outer) boundary
    select case (nright)

    case (0,162)        ! grid symmetric with right (outer) edge

      do n = 1, 6
        dxg(nmax + n)     = dxg(nmax + 1 - n)
        xag(nmax + n + 1) = xag(nmax + n) + dxg(nmax + n)
        xcg(nmax + n)     = xag(nmax + n) + half * dxg(nmax + n)
      end do ! n = 1, 6

    case (3)            ! grid symetric with left (inner) edge

      do n = 1, 6
        dxg(nmax + n)     = dxg(nmin + n - 1)
        xag(nmax + n + 1) = xag(nmax + n) + dxg(nmax + n)
        xcg(nmax + n)     = xag(nmax + n) + half * dxg(nmax + n)
      end do ! n = 1, 6

    case default

      call bl_error("invalid value for boundary condition, nright")

    end select

    return
  end subroutine grid_bc

  subroutine interp1dvol_chimera( x_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp1d_linear, interp1d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:)
    real (rt), intent(in) :: state_chim(:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out))

    ! local variables
    integer :: ix(size(x_out))
    integer :: ix_max

    real (rt) :: dvol1, dvol2
    real (rt) :: volx_out(size(x_out))

    integer :: i, n

    volx_out = third * x_out**3
    ix_max = imax_chim

    if ( interp_method == 1 ) then
      ix(:) = 0
      do i = 1, size(x_out)
        if ( volx_out(i) <= volx_c_chim(1) ) then
          ix(i) = 0
        else if ( volx_out(i) >= volx_c_chim(ix_max) ) then
          ix(i) = ix_max
        else
          ix(i) = locate( volx_out(i), ix_max, volx_c_chim ) - 1
        end if
      end do
      call interp1d_linear( ix, ix_max, volx_c_chim, state_chim, volx_out, state_out )
    else if ( interp_method == 2 ) then
      call interp1d_spline( volx_c_chim, state_chim, volx_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp1dvol_chimera

  subroutine interp1drad_chimera( x_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp1d_linear, interp1d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:)
    real (rt), intent(in) :: state_chim(:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out))

    ! local variables
    integer :: ix(size(x_out))
    integer :: ix_max

    integer :: i, n

    ix_max = imax_chim

    if ( interp_method == 1 ) then
      ix(:) = 0
      do i = 1, size(x_out)
        if ( x_out(i) <= x_c_chim(1) ) then
          ix(i) = 0
        else if ( x_out(i) >= x_c_chim(ix_max) ) then
          ix(i) = ix_max
        else
          ix(i) = locate( x_out(i), ix_max, x_c_chim ) - 1
        end if
      end do
      call interp1d_linear( ix, ix_max, x_c_chim, state_chim, x_out, state_out )
    else if ( interp_method == 2 ) then
      call interp1d_spline( x_c_chim, state_chim, x_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp1drad_chimera

  subroutine interp2dvol_chimera( x_out, y_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp2d_linear, interp2d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:,:)
    real (rt), intent(in) :: y_out(:,:)
    real (rt), intent(in) :: state_chim(:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,1),size(x_out,2))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2))
    integer :: iy_max

    real (rt) :: volx_out(size(x_out,1),size(x_out,2))
    real (rt) :: voly_out(size(x_out,1),size(x_out,2))

    real (rt) :: dvol1, dvol2

    integer :: i, j, n

    volx_out   = third * x_out**3
    voly_out = one - cos( y_out )
    ix_max = imax_chim
    iy_max = jmax_chim

    if ( interp_method == 1 ) then

      ix(:,:) = 0
      iy(:,:) = 0
      do j = 1, size(x_out,2)
        do i = 1, size(x_out,1)
          if ( volx_out(i,j) <= volx_c_chim(1) ) then
            ix(i,j) = 0
          else if ( volx_out(i,j) >= volx_c_chim(ix_max) ) then
            ix(i,j) = ix_max
          else
            ix(i,j) = locate( volx_out(i,j), ix_max, volx_c_chim ) - 1
          end if
          if ( voly_out(i,j) <= voly_c_chim(1) ) then
            iy(i,j) = 0
          else if ( voly_out(i,j) >= voly_c_chim(iy_max) ) then
            iy(i,j) = iy_max
          else
            iy(i,j) = locate( voly_out(i,j), iy_max, voly_c_chim ) - 1
          end if
        end do
      end do

      call interp2d_linear( ix, ix_max, iy, iy_max, volx_c_chim, voly_c_chim, state_chim, &
      &                     volx_out, voly_out, state_out )

    else if ( interp_method == 2 ) then
      call interp2d_spline( volx_c_chim(1:ix_max), voly_c_chim(1:iy_max), state_chim(1:ix_max,1:iy_max), &
      &                     volx_out, voly_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp2dvol_chimera

  subroutine interp2drad_chimera( x_out, y_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp2d_linear, interp2d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:,:)
    real (rt), intent(in) :: y_out(:,:)
    real (rt), intent(in) :: state_chim(:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,1),size(x_out,2))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2))
    integer :: iy_max

    integer :: i, j, n

    ix_max = imax_chim
    iy_max = jmax_chim

    if ( interp_method == 1 ) then

      ix(:,:) = 0
      iy(:,:) = 0
      do j = 1, size(x_out,2)
        do i = 1, size(x_out,1)
          if ( x_out(i,j) <= x_c_chim(1) ) then
            ix(i,j) = 0
          else if ( x_out(i,j) >= x_c_chim(ix_max) ) then
            ix(i,j) = ix_max
          else
            ix(i,j) = locate( x_out(i,j), ix_max, x_c_chim ) - 1
          end if
          if ( y_out(i,j) <= y_c_chim(1) ) then
            iy(i,j) = 0
          else if ( y_out(i,j) >= y_c_chim(iy_max) ) then
            iy(i,j) = iy_max
          else
            iy(i,j) = locate( y_out(i,j), iy_max, y_c_chim ) - 1
          end if
        end do
      end do

      call interp2d_linear( ix, ix_max, iy, iy_max, x_c_chim, y_c_chim, state_chim, &
      &                     x_out, y_out, state_out )

    else if ( interp_method == 2 ) then
      call interp2d_spline( x_c_chim(1:ix_max), y_c_chim(1:iy_max), state_chim(1:ix_max,1:iy_max), &
      &                     x_out, y_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp2drad_chimera

  subroutine interp2dvol_quad_chimera( x_out, y_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use model_interp_module, only: interp2d_linear, interp2d_spline
    use interpolate_module, only: locate
    use probdata_module, only: interp_method, nquad
    use quadrature_module, only: wquad, quad_avg

    ! input variables
    real (rt), intent(in) :: x_out(:,:,:,:)
    real (rt), intent(in) :: y_out(:,:,:,:)
    real (rt), intent(in) :: state_chim(:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,3),size(x_out,4))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2))
    integer :: iy_max

    real (rt) :: volx_quad(size(x_out,1),size(x_out,2))
    real (rt) :: voly_quad(size(x_out,1),size(x_out,2))

    real (rt) :: dvol1, dvol2
    real (rt) :: state_quad(size(x_out,1),size(x_out,1))

    integer :: i, ii, j, jj, n

    ix_max = imax_chim
    iy_max = jmax_chim

    if ( interp_method == 1 ) then

      do j = 1, size(x_out,4)
        do i = 1, size(x_out,3)
          volx_quad   = third * x_out(:,:,i,j)**3
          voly_quad = one - cos( y_out(:,:,i,j) )

          ix(:,:) = 0
          iy(:,:) = 0
          do jj = 1, nquad
            do ii = 1, nquad
              if ( volx_quad(ii,jj) <= volx_c_chim(1) ) then
                ix(ii,jj) = 0
              else if ( volx_quad(ii,jj) >= volx_c_chim(ix_max) ) then
                ix(ii,jj) = ix_max
              else
                ix(ii,jj) = locate( volx_quad(ii,jj), ix_max, volx_c_chim ) - 1
              end if
              if ( voly_quad(ii,jj) <= voly_c_chim(1) ) then
                iy(ii,jj) = 0
              else if ( voly_quad(ii,jj) >= voly_c_chim(iy_max) ) then
                iy(ii,jj) = iy_max
              else
                iy(ii,jj) = locate( voly_quad(ii,jj), iy_max, voly_c_chim ) - 1
              end if
            end do
          end do

          call interp2d_linear( ix, ix_max, iy, iy_max, volx_c_chim, voly_c_chim, state_chim, &
          &                     volx_quad, voly_quad, state_quad )

          state_out(i,j) = quad_avg( wquad, state_quad )

        end do
      end do

    else if ( interp_method == 2 ) then
      do j = 1, size(x_out,4)
        do i = 1, size(x_out,3)
          volx_quad   = third * x_out(:,:,i,j)**3
          voly_quad = one - cos( y_out(:,:,i,j) )

          call interp2d_spline( volx_c_chim(1:ix_max), voly_c_chim(1:iy_max), state_chim(1:ix_max,1:iy_max), &
          &                     volx_quad, voly_quad, state_quad )

          state_out(i,j) = quad_avg( wquad, state_quad )

        end do
      end do
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp2dvol_quad_chimera

  subroutine interp2drad_quad_chimera( x_out, y_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use model_interp_module, only: interp2d_linear, interp2d_spline
    use interpolate_module, only: locate
    use probdata_module, only: interp_method, nquad
    use quadrature_module, only: wquad, quad_avg

    ! input variables
    real (rt), intent(in) :: x_out(:,:,:,:)
    real (rt), intent(in) :: y_out(:,:,:,:)
    real (rt), intent(in) :: state_chim(:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,3),size(x_out,4))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2))
    integer :: iy_max

    real (rt) :: x_quad(size(x_out,1),size(x_out,2))
    real (rt) :: y_quad(size(x_out,1),size(x_out,2))

    real (rt) :: state_quad(size(x_out,1),size(x_out,1))

    integer :: i, ii, j, jj, n

    ix_max = imax_chim
    iy_max = jmax_chim

    if ( interp_method == 1 ) then

      do j = 1, size(x_out,4)
        do i = 1, size(x_out,3)
          x_quad   = x_out(:,:,i,j)
          y_quad =   y_out(:,:,i,j)

          ix(:,:) = 0
          iy(:,:) = 0
          do jj = 1, nquad
            do ii = 1, nquad
              if ( x_quad(ii,jj) <= x_c_chim(1) ) then
                ix(ii,jj) = 0
              else if ( x_quad(ii,jj) >= x_c_chim(ix_max) ) then
                ix(ii,jj) = ix_max
              else
                ix(ii,jj) = locate( x_quad(ii,jj), ix_max, x_c_chim ) - 1
              end if
              if ( y_quad(ii,jj) <= y_c_chim(1) ) then
                iy(ii,jj) = 0
              else if ( y_quad(ii,jj) >= y_c_chim(iy_max) ) then
                iy(ii,jj) = iy_max
              else
                iy(ii,jj) = locate( y_quad(ii,jj), iy_max, y_c_chim ) - 1
              end if
            end do
          end do

          call interp2d_linear( ix, ix_max, iy, iy_max, x_c_chim, y_c_chim, state_chim, &
          &                     x_quad, y_quad, state_quad )

          state_out(i,j) = quad_avg( wquad, state_quad )

        end do
      end do

    else if ( interp_method == 2 ) then
      do j = 1, size(x_out,4)
        do i = 1, size(x_out,3)
          x_quad = x_out(:,:,i,j)
          y_quad = y_out(:,:,i,j)

          call interp2d_spline( x_c_chim(1:ix_max), y_c_chim(1:iy_max), state_chim(1:ix_max,1:iy_max), &
          &                     x_quad, y_quad, state_quad )

          state_out(i,j) = quad_avg( wquad, state_quad )

        end do
      end do
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp2drad_quad_chimera

  subroutine interp3dvol_chimera( x_out, y_out, z_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp3d_linear, interp3d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:,:,:)
    real (rt), intent(in) :: y_out(:,:,:)
    real (rt), intent(in) :: z_out(:,:,:)
    real (rt), intent(in) :: state_chim(:,:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,1),size(x_out,2),size(x_out,3))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: iy_max

    integer :: iz(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: iz_max

    real (rt) :: volx_out(size(x_out,1),size(x_out,2),size(x_out,3))
    real (rt) :: voly_out(size(x_out,1),size(x_out,2),size(x_out,3))
    real (rt) :: volz_out(size(x_out,1),size(x_out,2),size(x_out,3))

    real (rt) :: dvol1, dvol2

    integer :: i, j, k, n

    volx_out   = third * x_out**3
    voly_out = one - cos( y_out )
    volz_out = z_out
    ix_max = imax_chim
    iy_max = jmax_chim
    iz_max = kmax_chim

    if ( interp_method == 1 ) then

      ix = 0
      iy = 0
      iz = 0
      do k = 1, size(x_out,3)
        do j = 1, size(x_out,2)
          do i = 1, size(x_out,1)
            if ( volx_out(i,j,k) <= volx_c_chim(1) ) then
              ix(i,j,k) = 0
            else if ( volx_out(i,j,k) >= volx_c_chim(ix_max) ) then
              ix(i,j,k) = ix_max
            else
              ix(i,j,k) = locate( volx_out(i,j,k), ix_max, volx_c_chim ) - 1
            end if

            if ( voly_out(i,j,k) <= voly_c_chim(1) ) then
              iy(i,j,k) = 0
            else if ( voly_out(i,j,k) >= voly_c_chim(iy_max) ) then
              iy(i,j,k) = iy_max
            else
              iy(i,j,k) = locate( voly_out(i,j,k), iy_max, voly_c_chim ) - 1
            end if

            if ( volz_out(i,j,k) <= volz_c_chim(1) ) then
              iz(i,j,k) = 0
            else if ( volz_out(i,j,k) >= volz_c_chim(iz_max) ) then
              iz(i,j,k) = iz_max
            else
              iz(i,j,k) = locate( volz_out(i,j,k), iz_max, volz_c_chim ) - 1
            end if
          end do
        end do
      end do

      call interp3d_linear( ix, ix_max, iy, iy_max, iz, iz_max, &
      &                     volx_c_chim, voly_c_chim, volz_c_chim, state_chim, &
      &                     volx_out, voly_out, volz_out, state_out )

!   else if ( interp_method == 2 ) then
!     call interp3d_spline( volx_c_chim(1:ix_max), voly_c_chim(1:iy_max), volz_c_chim(1:iz_max), &
!     &                     state_chim(1:ix_max,1:iy_max,1:iz_max), &
!     &                     volx_out, voly_out, volz_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp3dvol_chimera

  subroutine interp3drad_chimera( x_out, y_out, z_out, state_chim, state_out )

    use amrex_constants_module
    use amrex_error_module
    use interpolate_module, only: locate
    use model_interp_module, only: interp3d_linear, interp3d_spline
    use probdata_module, only: interp_method

    ! input variables
    real (rt), intent(in) :: x_out(:,:,:)
    real (rt), intent(in) :: y_out(:,:,:)
    real (rt), intent(in) :: z_out(:,:,:)
    real (rt), intent(in) :: state_chim(:,:,:)

    ! output variables
    real (rt), intent(out) :: state_out(size(x_out,1),size(x_out,2),size(x_out,3))

    ! local variables
    integer :: ix(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: ix_max

    integer :: iy(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: iy_max

    integer :: iz(size(x_out,1),size(x_out,2),size(x_out,3))
    integer :: iz_max

    integer :: i, j, k, n

    ix_max = imax_chim
    iy_max = jmax_chim
    iz_max = kmax_chim

    if ( interp_method == 1 ) then

      ix = 0
      iy = 0
      iz = 0
      do k = 1, size(x_out,3)
        do j = 1, size(x_out,2)
          do i = 1, size(x_out,1)
            if ( x_out(i,j,k) <= x_c_chim(1) ) then
              ix(i,j,k) = 0
            else if ( x_out(i,j,k) >= x_c_chim(ix_max) ) then
              ix(i,j,k) = ix_max
            else
              ix(i,j,k) = locate( x_out(i,j,k), ix_max, x_c_chim ) - 1
            end if

            if ( y_out(i,j,k) <= y_c_chim(1) ) then
              iy(i,j,k) = 0
            else if ( y_out(i,j,k) >= y_c_chim(iy_max) ) then
              iy(i,j,k) = iy_max
            else
              iy(i,j,k) = locate( y_out(i,j,k), iy_max, y_c_chim ) - 1
            end if

            if ( z_out(i,j,k) <= z_c_chim(1) ) then
              iz(i,j,k) = 0
            else if ( z_out(i,j,k) >= z_c_chim(iz_max) ) then
              iz(i,j,k) = iz_max
            else
              iz(i,j,k) = locate( z_out(i,j,k), iz_max, z_c_chim ) - 1
            end if
          end do
        end do
      end do

      call interp3d_linear( ix, ix_max, iy, iy_max, iz, iz_max, &
      &                     x_c_chim, y_c_chim, z_c_chim, state_chim, &
      &                     x_out, y_out, z_out, state_out )

!   else if ( interp_method == 2 ) then
!     call interp3d_spline( volx_c_chim(1:ix_max), voly_c_chim(1:iy_max), volz_c_chim(1:iz_max), &
!     &                     state_chim(1:ix_max,1:iy_max,1:iz_max), &
!     &                     volx_out, voly_out, volz_out, state_out )
    else
      call bl_error("invalid value for interp_method")
    end if

    return
  end subroutine interp3drad_chimera

end module chimera_parser_module


text:
program ch1009
    ! Timing of Reading Unformatted Files
    implicit none
    integer, parameter :: n = 1000000
    integer, dimension(1:n) :: x
    real, dimension(1:n) :: y
    integer :: i
    real :: t, t1, t2, t3
    character(15) :: comment

    call cpu_time(t)
    t1 = t
    comment = ' Program starts'
    print 100, comment, t1

    open(10,file='ch0914.dat',form='unformatted',&
        status='old')
    read(10) x
    call cpu_time(t)
    t2 = t-t1
    comment = ' Integer read '
    print 100, comment, t2

    read(10) y
    call cpu_time(t)
    t3 = t-t1-t2
    comment = ' Real read '
    print 100, comment, t3
    close(10)

    do i = 1, 10
        print 110, x(i), y(i)
    end do
100 format(1x,a,2x,f7.3)
110 format(1x,i4,2x,f10.7)
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
great_ne(jj1,kk1).
great_ne(y1,b1).
great_ne(ee1,y1).
great_ne(hh1,v1).
great_ne(bb1,v1).
great_ne(m1,h1).
great_ne(t1,ii1).
great_ne(kk1,dd1).
great_ne(z1,cc1).
great_ne(dd1,x1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
SUBROUTINE inner(g)

!-------------------------------------------------------------
!
!  Directs the inner iterations
!   Calls for the mesh sweep in 'sweep'
!   Evaulates convergence
!
!-------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: i, j, k, l, it
REAL*8 :: df, dfmx

! Initialize the previous flux iterate
e = 0.0
! Start the iterations
DO it = 1, itmx
   ! Call for the mesh sweep
   CALL sweep(g)

   ! Compare new and old flux iterates for user chosen range of moments, iall
   dfmx = -1.0
   DO l = 0, iall
      DO k = 0, iall
         DO j = 1, ny
            DO i = 1, nx
               ! Compute the difference depending on 'e' value
               IF (e(i,j,k,l) >= tolr) THEN
                  df = ABS((f(i,j,k,l,g) - e(i,j,k,l))/e(i,j,k,l))
               ELSE
                  df = ABS(f(i,j,k,l,g) - e(i,j,k,l))
               END IF
               ! Find the largest value
               IF (df > dfmx) THEN
                  dfmx = df
               END IF
            END DO
         END DO
      END DO
   END DO

   ! Print whether or not convergence was reached
   IF (dfmx > err .AND. it < itmx) THEN
      ! Set previous iterate of flux equal to current iterate
      e = f(:,:,:,:,g)
   ELSE IF (dfmx < err) THEN
      cnvf(g) = 1
      EXIT
   ELSE IF (it == itmx) THEN
      cnvf(g) = 0
      EXIT
   END IF

! End the iterations
END DO

RETURN
END SUBROUTINE inner

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text: IM ANDREA BUT I GO BY ANDI! I WANT TO PARTY IN DAVIS BUT MY FRIENDS DONT WANT TO!! IM MOVING TO DAVIS IN JANUARY 2012 AND I NEED DAVIS FRIENDS BUT DONT WANT TO SEEM DESPERATE HAHA HELP ME OUT, LETS GO PARTY AND ILL DEFINITELY TEACH YOU HOW TO DOUGIE

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-options " -O -findirect-inlining" }
function more_OK (fcn)
  character(*) more_OK
  character (*), external :: fcn
  more_OK = fcn ()
end function more_OK
  character(4) :: answer
  character(4), external :: is_OK, more_OK
  answer = more_OK (is_OK)
contains
END


text:
!-----------------------------------------------------------------------
! Interface file for vpush2.f
      module vpush2_h
      implicit none
!
      interface
         subroutine DISTR2T(part,vtx,vty,vdx,vdy,npx,npy,idimp,npe,nx,ny&
     &,ipbc)
         implicit none
         integer :: npx, npy, idimp, npe, nx, ny, ipbc
         real :: vtx, vty, vdx, vdy
         real, dimension(npe,idimp) :: part
         end subroutine
      end interface
!
      interface
         subroutine GPUSH2LT(part,fxy,qbm,dt,ek,idimp,nop,npe,nx,ny,nxv,&
     &nyv,ipbc)
         implicit none
         integer :: idimp, nop, npe, nx, ny, nxv, nyv, ipbc
         real :: qbm, dt, ek
         real, dimension(npe,idimp) :: part
         real, dimension(2,nxv,nyv) :: fxy
         end subroutine
      end interface
!
      interface
         subroutine VGPUSH2LT(part,fxy,qbm,dt,ek,idimp,nop,npe,nx,ny,nxv&
     &,nyv,ipbc)
         implicit none
         integer :: idimp, nop, npe, nx, ny, nxv, nyv, ipbc
         real :: qbm, dt, ek
         real, dimension(npe,idimp) :: part
         real, dimension(2,nxv,nyv) :: fxy
         end subroutine
      end interface
!
      interface
         subroutine GPOST2LT(part,q,qm,nop,npe,idimp,nxv,nyv)
         implicit none
         integer :: nop, npe, idimp, nxv, nyv
         real :: qm
         real, dimension(npe,idimp) :: part
         real, dimension(nxv,nyv) :: q
         end subroutine
      end interface
!
      interface
         subroutine VGPOST2LT(part,q,qm,nop,npe,idimp,nxv,nyv)
         implicit none
         integer :: nop, npe, idimp, nxv, nyv
         real :: qm
         real, dimension(npe,idimp) :: part
         real, dimension(nxv,nyv) :: q
         end subroutine
      end interface
!
      interface
         subroutine DSORTP2YLT(parta,partb,npic,idimp,nop,npe,ny1)
         implicit none
         integer :: idimp, nop, npe, ny1
         real, dimension(npe,idimp) :: parta, partb
         integer, dimension(ny1) :: npic
         end subroutine
      end interface
!
      interface
         subroutine CGUARD2L(fxy,nx,ny,nxe,nye)
         implicit none
         integer :: nx, ny, nxe, nye
         real, dimension(2,nxe,nye) :: fxy
         end subroutine
      end interface
!
      interface
         subroutine AGUARD2L(q,nx,ny,nxe,nye)
         implicit none
         integer :: nx, ny, nxe, nye
         real, dimension(nxe,nye) :: q
         end subroutine
      end interface
!
      interface
         subroutine VPOIS22(q,fxy,isign,ffc,ax,ay,affp,we,nx,ny,nxvh,nyv&
     &,nxhd,nyhd)
         implicit none
         integer :: isign, nx, ny, nxvh, nyv, nxhd, nyhd
         real :: ax, ay, affp ,we
         real, dimension(2*nxvh,nyv) :: q
         real, dimension(2,2*nxvh,nyv) :: fxy
         complex, dimension(nxhd,nyhd) :: ffc
         end subroutine
      end interface
!
      interface
         subroutine WFFT2RINIT(mixup,sct,indx,indy,nxhyd,nxyhd)
         implicit none
         integer :: indx, indy, nxhyd, nxyhd
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine WFFT2RVX(f,isign,mixup,sct,indx,indy,nxhd,nyd,nxhyd,&
     &nxyhd)
         implicit none
         integer :: isign, indx, indy, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine WFFT2RV2(f,isign,mixup,sct,indx,indy,nxhd,nyd,nxhyd,&
     &nxyhd)
         implicit none
         integer :: isign, indx, indy, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2,2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine FFT2RVXX(f,isign,mixup,sct,indx,indy,nyi,nyp,nxhd,  &
     &nyd,nxhyd,nxyhd)
         implicit none
         integer :: isign, indx, indy, nyi, nyp, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine FFT2RXY(f,isign,mixup,sct,indx,indy,nxi,nxp,nxhd,nyd&
     &,nxhyd,nxyhd)
         implicit none
         integer :: isign, indx, indy, nxi, nxp, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine FFT2RV2X(f,isign,mixup,sct,indx,indy,nyi,nyp,nxhd,  &
     &nyd,nxhyd,nxyhd)
         implicit none
         integer :: isign, indx, indy, nyi, nyp, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2,2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         subroutine FFT2R2Y(f,isign,mixup,sct,indx,indy,nxi,nxp,nxhd,nyd&
     &,nxhyd,nxyhd)
         implicit none
         integer :: isign, indx, indy, nxi, nxp, nxhd, nyd, nxhyd, nxyhd
         real, dimension(2,2*nxhd,nyd) :: f
         integer, dimension(nxhyd) :: mixup
         complex, dimension(nxyhd) :: sct
         end subroutine
      end interface
!
      interface
         function ranorm()
         implicit none
         double precision :: ranorm
         end function
      end interface
!
      interface
         function randum()
         implicit none
         double precision :: randum
         end function
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20100912 16:31:08 nbsp Welcome to the Wiki Howdy Susan, and a belated welcome to the wiki!   My names Evan, pleased to meet you!  I read your really helpful comment on Yolo County SPCA and thought it was so insightful and helpful that it really should be at the top of the entry, so I moved it up there.  Keep in mind that everything on the wiki was written by people like you and me... dont feel that you are only allowed to comment.  You can edit any entry and add your voice right into the top part; it was all written by normal folk, and everybody has an equal opportunity to write and edit anything on the wiki.  Thanks again for the clarification on what is probably a misunderstood (and apparently often seen as hostile) process.  Getting the internal perspective helps people understand what is going on.  Once again, welcome to the wiki! Users/JabberWokky Evan JabberWokky Edwards


text:
	SUBROUTINE DCSELS_OUT ( iflno, dattim, rdata, iret)
C************************************************************************
C* DCSELS								*
C*									*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'BRIDGE.PRM'
C
	INTEGER		iflno,iret
	CHARACTER	dattim*(*)
	REAL		rdata(*)
C
	REAL		rdout(4)
	CHARACTER	stype,errstr*132
	INTEGER		ihour,idno
	CHARACTER	stid*8,ostid*8
	REAL		slat,slon,selv
	CHARACTER	state*2,country*2

	DATA	state,country /'--','--'/
C************************************************************************
        iret = 0
C
	IF (rdata(1).eq.1) THEN
           stype = 'T'
        ELSE IF (rdata(1).eq.2) THEN
	   stype = 'W'
        ELSE IF (rdata(1).eq.3) THEN
	   stype = 'A'
	ELSE IF (rdata(1).eq.4) THEN
	   stype = 'G'
	ELSE
	   stype = 'U'
	END IF
C
        idno = 0
        ihour = int(rdata(4))
C
C*      See if time is already in file
C
	CALL SF_STIM (iflno, dattim, ier)
C
        IF (ier.ne.0) THEN
           write(stid,1000) stype,ihour,idno
        ELSE
           DO WHILE ((ier.eq.0).and.(idno.lt.1000))
              write(stid,1000) stype,ihour,idno
1000          FORMAT(1A,I4.4,I3.3)
              CALL SF_SSTN (iflno, stid, ostid, isnum, slat, slon,
     +                      selv, ispri, ier)
              IF (ier.eq.0) THEN
                 IF ((abs(slat-rdata(5)).gt.0.012).or.
     +              (abs(slon-rdata(6)).gt.0.012)) THEN
                    idno = idno + 1
                 ELSE
                    idno = 1000
                 END IF
              END IF
           END DO
        END IF

        IF (idno.eq.1000) THEN
           write(errstr,1001) stid, rdata(5), rdata(6)
1001	   FORMAT('Duplicate report ',A8,' Lat = ',F6.2,' Lon = ',F7.2)
           CALL DC_WLOG(2,'DCSELSOUT',0,errstr,iout)
        ELSE
           write(errstr,1002) stid, rdata(5), rdata(6)
1002	   FORMAT('Storm: ',A8,' Lat = ',F6.2,' Lon = ',F7.2)
           CALL DC_WLOG(1,'DCSELSOUT',0,errstr,iout)
           slat = rdata(5)
           slon = rdata(6)
           selv = 0
           isnum = rdata(3)*10000 + rdata(4)
	   rdout(1) = rdata(1)
           rdout(2) = rdata(2)
           rdout(3) = rdata(3)
           rdout(4) = rdata(4)
           CALL SF_WSDD (iflno, dattim, stid, isnum, slat, slon,
     +                   selv,state,country,ihour,rdout,iret)
        END IF
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


program problem_4

use aoc_utilities

implicit none

type :: board
  integer,dimension(5,5) :: nums = 0
  logical,dimension(5,5) :: mask = .false.
  logical :: won = .false.
end type board

type(board),dimension(:),allocatable :: boards
integer :: iunit, n_lines, i, n_boards, n_nonblank_lines, iboard, j, k
character(len=:),allocatable :: line
logical :: status_ok
integer,dimension(:),allocatable :: draws
type(string),dimension(:),allocatable :: vals
integer :: winning_draw, winning_board, winning_draw_num, score
integer :: num_not_won

open(newunit=iunit,file='inputs/day4.txt')

! count boards
n_lines = number_of_lines_in_file(iunit)
n_nonblank_lines = 0
do i = 1, n_lines
    call read_line_from_file(iunit,line,status_ok)
    if (line/='') n_nonblank_lines = n_nonblank_lines + 1
end do
n_boards = (n_nonblank_lines - 1) / 5
! write(*,*) 'n_lines           = ', n_lines
! write(*,*) 'n_nonblank_lines  = ', n_nonblank_lines
! write(*,*) 'n_boards          = ', n_boards

allocate(boards(n_boards))
iboard = 0
j = 0
rewind(iunit)
do i = 1, n_lines
    call read_line_from_file(iunit,line,status_ok)
    if (.not. allocated(line)) cycle
    if (i==1) then
        ! numbers to be drawn:
        call split(line,',',vals)
        allocate(draws(size(vals)))
        do j = 1, size(vals)
            read(vals(j)%str, *) draws(j)
        end do
    else if (line == '') then
        iboard = iboard + 1
        j = 0
        cycle
    else
        ! read the jth line of the iboard
        j = j + 1
        read(line, '(I2,1X,I2,1X,I2,1X,I2,1X,I2))') boards(iboard)%nums(j,:)
    end if
end do

! write(*,*) boards(1)%nums

! now, play the game:
main: do i = 1, size(draws)

    ! each board:
    do j = 1, size(boards)
        where (boards(j)%nums == draws(i))
            boards(j)%mask = .true.
        end where
    end do

    ! check for any wins:
    do j = 1, size(boards)
        do k = 1, 5
            if (all(boards(j)%mask(:,k))) then ! cols
                winning_draw_num = i
                winning_draw = draws(i)
                winning_board = j
                exit main
            elseif (all(boards(j)%mask(k,:))) then ! rows
                winning_draw_num = i
                winning_draw = draws(i)
                winning_board = j
                exit main
            end if
        end do
    end do

end do main

write(*,*) 'board', winning_board, 'wins on draw', winning_draw_num, '(', winning_draw,')'

! write(*,*) boards(winning_board)%mask
!The score of the winning board can now be calculated.
!Start by finding the sum of all unmarked numbers on that board;
!in this case, the sum is 188. Then, multiply that sum by the
!number that was just called when the board won, 24, to get
!the final score, 188 * 24 = 4512.

score = 0
do i = 1, 5
    score = score + sum( boards(winning_board)%nums(:,i), mask = .not. boards(winning_board)%mask(:,i) )
end do
score = score * winning_draw

write(*,*) '4A: ', score

!**************** find the last to win:

do i = 1, size(boards)
    boards(i)%mask = .false. ! reset boards
end do

! now, play the game:
main2: do i = 1, size(draws)

    ! each board:
    do j = 1, size(boards)
        where (boards(j)%nums == draws(i))
            boards(j)%mask = .true.
        end where
    end do

    ! check for any wins:
    do j = 1, size(boards)
        do k = 1, 5
            if (all(boards(j)%mask(:,k))) then ! cols
                boards(j)%won = .true.
            elseif (all(boards(j)%mask(k,:))) then ! rows
                boards(j)%won = .true.
            end if
        end do
    end do

    ! check if only one left that hasn't won:
    if ( count(boards(:)%won) == n_boards-1 ) then ! only one left:
        do j = 1, size(boards)
            if (.not. boards(j)%won) then
                winning_board = j
            end if
        end do
    elseif (count(boards(:)%won) == n_boards) then ! last one done
        winning_draw_num = i
        winning_draw = draws(i)
        exit main2
    end if

end do main2

write(*,*) 'board', winning_board, 'last win on draw', winning_draw_num, '(', winning_draw,')'

score = 0
do i = 1, 5
    score = score + sum( boards(winning_board)%nums(:,i), mask = .not. boards(winning_board)%mask(:,i) )
end do
score = score * winning_draw

write(*,*) '4B: ', score

end program problem_4