Merge pull request #6 from starkiller-astro/scalar_vector_interfaces

Scalar/vector interfaces

source/model_input_ascii.f90

@@ -161,13 +161,17 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in )
 
     ! Initialize
     ierr = 0
-    ! Allow auxiliary nucleus
-    ! Set to 0 to force normalization of abundances
-    ! and exclude aux nuclear species
-    iaux = 1
     Do izb = zb_lo, zb_hi
       yestart(izb) = 0.0     
       ystart(:,izb) = 0.0
+
+      ! Allow auxiliary nucleus
+      ! Set to 0 to force normalization of abundances
+      ! and exclude aux nuclear species
+      iaux(izb) = 1
+      xext(izb) = 0.0
+      aext(izb) = 1.0
+      zext(izb) = 0.0
     EndDo
 
     Do izb = zb_lo, zb_hi
@@ -190,20 +194,17 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in )
 
         If ( t9start(izb) <= t9nse .or. yestart(izb) <= 0.0 .or. yestart(izb) >= 1.0 ) Then
           Call read_inab_file(inab_file(izone),abund_desc(izb),yein,yin,xext_loc,aext_loc,zext_loc,xnorm,ierr)
-          If ( iaux(izb)==0 ) then ! Turn off aux nucleus
-              yin = yin/xnorm
+          If ( iaux(izb) == 0 ) then ! Turn off aux nucleus
+              yin = yin / xnorm
               Write(lun_diag,"(a)") 'Normalizing initial abundances.'
-              xext(izb) = 0.0
-              aext(izb) = 1.0
-              zext(izb) = 0.0
           Else    
               xext(izb) = xext_loc
               aext(izb) = aext_loc
               zext(izb) = zext_loc
           Endif
           ! If Ye is not provided in the initial abundance file explicitly, calculate it here
           If ( yein <= 0.0 .or. yein >= 1.0 ) &
-          &    Call y_moment(yin,yein,ytot,abar,zbar,z2bar,zibar,izb)
+          &    Call y_moment(yin,yein,ytot,abar,zbar,z2bar,zibar,xext(izb),aext(izb),zext(izb))
           yestart(izb) = yein
 
 
@@ -223,9 +224,6 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in )
           Call nse_solve(rhostart(izb),t9start(izb),yestart(izb))
           ystart(:,izb) = ynse(:)
           iaux(izb) = 0
-          xext(izb) = 0.d0
-          aext(izb) = 1.d0
-          zext(izb) = 1.d0
         Else
           ystart(:,izb) = yin(:)
         EndIf

source/xnet_abundances.f90

@@ -20,51 +20,125 @@ Module xnet_abundances
   Real(dp), Allocatable :: aext(:)   ! Mass number of auxiliary/external species
   Real(dp), Allocatable :: zext(:)   ! Charge number of auxiliary/external species
 
+  Interface y_moment
+    Module Procedure y_moment_scalar
+    Module Procedure y_moment_vector
+  End Interface y_moment
+
+  Private :: y_moment_internal
+
 Contains
 
-  Subroutine y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb)
+  Subroutine y_moment_internal(y,ye,ytot,abar,zbar,z2bar,zibar,xext,yext,zext)
     !-----------------------------------------------------------------------------------------------
     ! This routine calculates moments of the abundance distribution for the EOS.
     !-----------------------------------------------------------------------------------------------
     Use nuclear_data, Only: ny, aa, zz, zz2, zzi
-    Use xnet_controls, Only: idiag, lun_diag
     Use xnet_constants, Only: thbim1
     Use xnet_types, Only: dp
     Implicit None
 
     ! Input variables
     Real(dp), Intent(in) :: y(ny)
-    Integer, Intent(in),Optional :: izb 
+    Real(dp), Intent(in) :: xext, yext, zext
 
     ! Output variables
     Real(dp), Intent(out) :: ye, ytot, abar, zbar, z2bar, zibar
 
     ! Local variables
-    Real(dp) :: atot, ztot, yext, xext_loc,zext_loc, aext_loc
-
-    If (present(izb)) Then
-        yext = xext(izb)/aext(izb)
-        xext_loc=xext(izb)
-        aext_loc=aext(izb)
-        zext_loc=zext(izb)
-    Else
-        yext=0.0
-        xext_loc=0.0
-        aext_loc=1.0
-        zext_loc=0.0
-    Endif
+    Real(dp) :: atot, ztot
+
     ! Calculate abundance moments
     ytot  = sum(y) + yext
-    atot  = sum(y(:) * aa(:)) + xext_loc
-    ztot  = sum(y * zz) + yext*zext_loc
+    atot  = sum(y * aa) + xext
+    ztot  = sum(y * zz) + yext * zext
     abar  = atot / ytot
     zbar  = ztot / ytot
-    z2bar = ( sum(y * zz2) + yext * zext_loc * zext_loc ) / ytot 
-    zibar = ( sum(y * zzi) + yext * zext_loc**thbim1 ) / ytot
-    ye = ztot / atot
-    If ( idiag >= 3 ) Write(lun_diag,"(a4,7es23.15)") 'YMom',ytot,abar,zbar,z2bar,zibar,ye, atot
+    z2bar = ( sum(y * zz2) + yext * zext * zext ) / ytot 
+    zibar = ( sum(y * zzi) + yext * zext**thbim1 ) / ytot
+    ye    = ztot / atot
+
+    Return
+  End Subroutine y_moment_internal
+
+  Subroutine y_moment_scalar(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
+    !-----------------------------------------------------------------------------------------------
+    ! This is the interface for the scalar version.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_controls, Only: idiag, lun_diag
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: y(ny)
+    Real(dp), Intent(in) :: xext, aext, zext
+
+    ! Output variables
+    Real(dp), Intent(out) :: ye, ytot, abar, zbar, z2bar, zibar
+
+    ! Local variables
+    Real(dp) :: yext
+
+    ! Calculate abundance moments
+    yext = xext / aext
+    Call y_moment_internal(y,ye,ytot,abar,zbar,z2bar,zibar,xext,yext,zext)
+    If ( idiag >= 3 ) Write(lun_diag,"(a4,6es23.15)") 'YMom',ytot,abar,zbar,z2bar,zibar,ye
+
+    Return
+  End Subroutine y_moment_scalar
+
+  Subroutine y_moment_vector(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This is the interface for the vector version.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: y(ny,zb_lo:zb_hi)
+    Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi)
+
+    ! Output variables
+    Real(dp), Intent(out) :: ye(zb_lo:zb_hi), ytot(zb_lo:zb_hi), abar(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: zbar(zb_lo:zb_hi), z2bar(zb_lo:zb_hi), zibar(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Real(dp) :: yext
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    ! Calculate abundance moments
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        yext = xext(izb) / aext(izb)
+        Call y_moment_internal(y(:,izb),ye(izb),ytot(izb), &
+          & abar(izb),zbar(izb),z2bar(izb),zibar(izb),xext(izb),yext,zext(izb))
+      EndIf
+    EndDo
+
+    If ( idiag >= 3 ) Then
+      Do izb = zb_lo, zb_hi
+        If ( mask(izb) ) Then
+          Write(lun_diag,"(a4,6es23.15)") 'YMom', &
+            ytot(izb),abar(izb),zbar(izb),z2bar(izb),zibar(izb),ye(izb)
+        EndIf
+      EndDo
+    EndIf
 
     Return
-  End Subroutine y_moment
+  End Subroutine y_moment_vector
 
 End Module xnet_abundances

source/xnet_conditions.f90

@@ -38,24 +38,82 @@ Module xnet_conditions
   Real(dp), Allocatable :: tstart(:),tstop(:),tdelstart(:),t9start(:),rhostart(:),yestart(:)
   Real(dp), Allocatable :: th(:,:),t9h(:,:),rhoh(:,:),yeh(:,:)
 
+  Interface t9rhofind
+    Module Procedure t9rhofind_scalar
+    Module Procedure t9rhofind_vector
+  End Interface t9rhofind
+
 Contains
 
-  Subroutine t9rhofind(kstep,tf,nf,t9f,rhof,mask_in)
+  Subroutine t9rhofind_scalar(kstep,tf,nf,t9f,rhof,ns,ts,t9s,rhos)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine calculates t9 and rho as a function of time, either via interpolation or from an
+    ! analytic expression. (scalar interface)
+    !-----------------------------------------------------------------------------------------------
+    Use, Intrinsic :: iso_fortran_env, Only: lun_stdout=>output_unit
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Integer, Intent(in) :: kstep, ns
+    Real(dp), Intent(in) :: tf, ts(:), t9s(:), rhos(:)
+
+    ! Output variables
+    Integer, Intent(out) :: nf
+    Real(dp), Intent(out) :: t9f, rhof
+
+    ! Local variables
+    Real(dp) :: dt, rdt, dt9, drho
+    Integer :: n
+
+    ! For constant conditions (ns = 1), set temperature and density
+    If ( ns == 1 ) Then
+      t9f = t9s(1)
+      rhof = rhos(1)
+      nf = 1
+
+    ! Otherwise, calculate T9 and rho by interpolation
+    Else
+      Do n = 1, ns
+        If ( tf <= ts(n) ) Exit
+      EndDo
+      nf = n
+      If ( n > 1 .and. n <= ns ) Then
+        rdt = 1.0 / (ts(n)-ts(n-1))
+        dt = tf - ts(n-1)
+        dt9 = t9s(n) - t9s(n-1)
+        drho = rhos(n) - rhos(n-1)
+        t9f = dt*rdt*dt9 + t9s(n-1)
+        rhof = dt*rdt*drho + rhos(n-1)
+      ElseIf ( n == 1 ) Then
+        t9f = t9s(1)
+        rhof = rhos(1)
+      Else
+        t9f = t9s(ns)
+        rhof = rhos(ns)
+        Write(lun_stdout,*) 'Time beyond thermodynamic range',tf,' >',ts(ns)
+      EndIf
+    EndIf
+
+    Return
+  End Subroutine t9rhofind_scalar
+
+  Subroutine t9rhofind_vector(kstep,tf,nf,t9f,rhof,mask_in)
     !-----------------------------------------------------------------------------------------------
     ! This routine calculates t9 and rho as a function of time, either via interpolation or from an
-    ! analytic expression.
+    ! analytic expression. (vector interface)
     !-----------------------------------------------------------------------------------------------
-    Use xnet_controls, Only: lun_diag, lun_stdout, nzevolve, zb_lo, zb_hi, lzactive
+    Use xnet_controls, Only: lun_diag, lun_stdout, zb_lo, zb_hi, lzactive
     Use xnet_types, Only: dp
     Implicit None
 
     ! Input variables
     Integer, Intent(in) :: kstep
-    Real(dp), Intent(in) :: tf(nzevolve)
+    Real(dp), Intent(in) :: tf(zb_lo:zb_hi)
 
     ! Input/Output variables
-    Integer, Intent(inout) :: nf(nzevolve)
-    Real(dp), Intent(inout) :: t9f(nzevolve), rhof(nzevolve)
+    Integer, Intent(inout) :: nf(zb_lo:zb_hi)
+    Real(dp), Intent(inout) :: t9f(zb_lo:zb_hi), rhof(zb_lo:zb_hi)
 
     ! Optional variables
     Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
@@ -107,6 +165,6 @@ Subroutine t9rhofind(kstep,tf,nf,t9f,rhof,mask_in)
     EndDo
 
     Return
-  End Subroutine t9rhofind
+  End Subroutine t9rhofind_vector
 
 End Module xnet_conditions