text:
      subroutine wtransform_wt(p,q,pbDpg,ptDpg,ptDpb)
      implicit none
      include 'constants.f'
      double precision p(mxpart,4),pw(4),pt(4),lDt(5:6),lDw(5:6)
      double precision ptDpt,pwDpw,ptDpw,q(mxpart,4),root,hsin,hcos,a,b
      double precision ptDpg,pbDpg,ptDpb

      integer j,nu

      do nu=1,4
         pw(nu)=p(5,nu)+p(6,nu)
         pt(nu)=pw(nu)+p(7,nu)+p(8,nu)
         do j=1,4
            q(j,nu)=p(j,nu)
         enddo
      enddo

      pbDpg=p(7,4)*p(8,4)-p(7,1)*p(8,1)-p(7,2)*p(8,2)-p(7,3)*p(8,3)
      ptDpg=pt(4)*p(8,4)-pt(1)*p(8,1)-pt(2)*p(8,2)-pt(3)*p(8,3)
      ptDpb=pt(4)*p(7,4)-pt(1)*p(7,1)-pt(2)*p(7,2)-pt(3)*p(7,3)
      ptDpw=pt(4)*pw(4)-pt(1)*pw(1)-pt(2)*pw(2)-pt(3)*pw(3)
      ptDpt=pt(4)**2-pt(1)**2-pt(2)**2-pt(3)**2
      pwDpw=pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2

      root=dsqrt(ptDpw**2-ptDpt*pwDpw)
      hsin=0.5d0/(ptDpt*pwDpw)*(-(ptDpt-pwDpw)*ptDpw+(ptDpt+pwDpw)*root)
      hcos=0.5d0/(ptDpt*pwDpw)*(+(ptDpt+pwDpw)*ptDpw-(ptDpt-pwDpw)*root)

C---calculate coefficients of lorentz transformation
      a=hsin/root
      b=(hcos-1d0)/root**2

c---dot t and w into decay products of w
      do j=5,6
         lDt(j)=p(j,4)*pt(4)-p(j,1)*pt(1)-p(j,2)*pt(2)-p(j,3)*pt(3)
         lDw(j)=p(j,4)*pw(4)-p(j,1)*pw(1)-p(j,2)*pw(2)-p(j,3)*pw(3)
      enddo

      do nu=1,4
      do j=5,6
      q(j,nu)=p(j,nu)+a*(pt(nu)*lDw(j)-pw(nu)*lDt(j))
     . +b*(ptDpw*(pt(nu)*ldw(j)+pw(nu)*ldt(j))
     .  -pwDpw*lDt(j)*pt(nu)-ptDpt*lDw(j)*pw(nu))
      enddo
      q(7,nu)=-q(1,nu)-q(2,nu)-q(3,nu)-q(4,nu)-q(5,nu)-q(6,nu)
      q(8,nu)=0d0
      enddo

      call storeptilde(1,q)

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
	SUBROUTINE CROSSP(A,B,C)
	DIMENSION A(3),B(3),C(3)
	A(1) = B(2)*C(3) - B(3)*C(2)
	A(2) = B(3)*C(1) - B(1)*C(3)
	A(3) = B(1)*C(2) - B(2)*C(1)
	RETURN
	END

ext: for


text:
      subroutine dkqqb_hww_gs(p,msqc)
      implicit none
************************************************************************
*     Author: R.K. Ellis                                               *
*     May, 2012.                                                       *
*     calculate the subtraction term for radiation in the              *
*     W^- decay for the process                                        *
*                                                                      *
*     q(-p1)+qbar(-p2) = H-> W^+(nu(p3)+e+(p4))                        *
*                       +W^-(e-(p5)+nubar(p6)+g(p7))                   *
*                                                                      *
************************************************************************
      include 'constants.f'
      include 'ptilde.f'
      include 'qqgg.f'
      include 'incldip.f'
      double precision msqc(maxd,-nf:nf,-nf:nf),p(mxpart,4)
      double precision
     & msq57_6(-nf:nf,-nf:nf),msq67_5(-nf:nf,-nf:nf),
     & dummyv(-nf:nf,-nf:nf),sub57_6(4),sub67_5(4),dsubv
      integer j,k,nd
      external qqb_hww,donothing_gvec

      ndmax=2

      do j=-nf,nf
      do k=-nf,nf
      do nd=1,ndmax
        msqc(nd,j,k)=0d0
        incldip(nd)=.true.
      enddo
      enddo
      enddo

      call dips(1,p,5,7,6,sub57_6,dsubv,msq57_6,dummyv,
     & qqb_hww,donothing_gvec)
      call dips(2,p,6,7,5,sub67_5,dsubv,msq67_5,dummyv,
     & qqb_hww,donothing_gvec)

      msqc(1,0,0)=sub57_6(qq)*msq57_6(0,0)*2d0*cf
      msqc(2,0,0)=sub67_5(qq)*msq67_5(0,0)*2d0*cf

      return
      end


text:
!*
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!*
!     test interface sub with subroutine sub
	module interf
	  public  sub
	  interface  sub
	    module procedure  sub, subr
	  end interface

	  contains
	    subroutine  sub(II)
	      integer  II
	      II = 122
	    end subroutine  sub
	    subroutine  subr(rr)
	      real  rr
	      rr = 99
	    end subroutine  subr
	end module interf

	program  testmod
	  use  interf
	  implicit none
	  integer  II
	  integer,dimension(2) :: result,expect
	  data expect/122,99/
	  real  rr
	  ii = 0
	  rr = 0

	  call sub(II)
	  call sub(rr)
	  result(1) = ii
	  result(2) = rr
	  call check(result,expect,2)
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
! “ButterflyPACK” Copyright (c) 2018, The Regents of the University of California, through
! Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from the
! U.S. Dept. of Energy). All rights reserved.

! If you have questions about your rights to use or distribute this software, please contact
! Berkeley Lab's Intellectual Property Office at  IPO@lbl.gov.

! NOTICE.  This Software was developed under funding from the U.S. Department of Energy and the
! U.S. Government consequently retains certain rights. As such, the U.S. Government has been
! granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable
! worldwide license in the Software to reproduce, distribute copies to the public, prepare
! derivative works, and perform publicly and display publicly, and to permit other to do so.

! Developers: Yang Liu
!             (Lawrence Berkeley National Lab, Computational Research Division).

! This exmple works with double precision data
#define DAT 1

#include "ButterflyPACK_config.fi"

module APPLICATION_MODULE
use BPACK_DEFS
implicit none

	!**** define your application-related variables here
	type quant_app
		real(kind=8) sigma, lambda ! Kernel Regression: RBF parameters
		integer dimn  ! dimension of the data sets
		integer ntrain,ntest ! size of training points and test points
		character(LEN=500)::trainfile_p,trainfile_tree,trainfile_l,testfile_p,testfile_l !Kernel Regression: file pointers to train data, preordered tree, train labels, test data, and test labels
		CHARACTER (LEN=1000) DATA_DIR
		integer Nunk ! size of the matrix
		real(kind=8),allocatable:: xyz(:,:)   ! coordinates of the points

	end type quant_app

contains

	!**** cutoff distance for gaussian kernel
	real(kind=8) function arg_thresh_Zmn(quant)
		use BPACK_DEFS
		implicit none

		type(quant_app)::quant
		arg_thresh_Zmn=-log(BPACK_SafeUnderflow)*2.0*quant%sigma**2
	end function arg_thresh_Zmn

	!**** user-defined subroutine to sample Z_mn
	subroutine Zelem_RBF(m,n,value_e,quant)
		use BPACK_DEFS
		implicit none

		class(*),pointer :: quant
		integer, INTENT(IN):: m,n
		real(kind=8)::value_e

		real(kind=8) r_mn
		integer dimn

		select TYPE(quant)
		type is (quant_app)
			dimn = size(quant%xyz,1)
			value_e=0
			r_mn=sum((quant%xyz(1:dimn,m)-quant%xyz(1:dimn,n))**2)
			if(r_mn<arg_thresh_Zmn(quant))then
				value_e = exp(-r_mn/2.0/quant%sigma**2)
			endif
			if(r_mn==0)value_e = value_e + quant%lambda
		class default
			write(*,*)"unexpected type"
			stop
		end select
	end subroutine Zelem_RBF

end module APPLICATION_MODULE

PROGRAM ButterflyPACK_KRR
    use BPACK_DEFS
    use APPLICATION_MODULE

	use BPACK_structure
	use BPACK_factor
	use BPACK_constr
	use omp_lib
	use MISC_Utilities

    implicit none

	! include "mkl_vml.fi"

    real(kind=8) para
    real(kind=8) tolerance
    integer Primary_block, nn, mm,kk,mn,rank,ii,jj
    integer i,j,k, threads_num,iii
	integer times(8)
	real(kind=8) t1,t2,x,y,z,r,theta,phi

	character(len=:),allocatable  :: string
	character(len=1024)  :: strings,strings1
	character(len=6)  :: info_env
	integer :: length
	integer :: ierr
	integer*8 oldmode,newmode
	type(Hoption)::option
	type(Hstat)::stats
	type(mesh)::msh
	type(kernelquant)::ker
	type(quant_app),target::quant
	type(Bmatrix)::bmat
	integer,allocatable:: groupmembers(:)
	integer nmpi
	type(proctree)::ptree
	integer MPI_thread
	integer,allocatable::Permutation(:)
	integer Nunk_loc
	integer nargs,flag
	integer v_major,v_minor,v_bugfix

	call MPI_Init(ierr)
	call MPI_Comm_size(MPI_Comm_World,nmpi,ierr)
	allocate(groupmembers(nmpi))
	do ii=1,nmpi
		groupmembers(ii)=(ii-1)
	enddo

	!**** create the process tree
	call CreatePtree(nmpi,groupmembers,MPI_Comm_World,ptree)
	deallocate(groupmembers)

	if(ptree%MyID==Main_ID)then
    write(*,*) "-------------------------------Program Start----------------------------------"
    write(*,*) "ButterflyPACK_KRR"
	call BPACK_GetVersionNumber(v_major,v_minor,v_bugfix)
	write(*,'(A23,I1,A1,I1,A1,I1,A1)') " ButterflyPACK Version:",v_major,".",v_minor,".",v_bugfix
    write(*,*) "   "
	endif

	!**** initialize stats and option
	call InitStat(stats)
	call SetDefaultOptions(option)

	!**** intialize the user-defined derived type quant
	quant%DATA_DIR='../EXAMPLE/KRR_DATA/susy_10Kn'
    quant%dimn=8
    quant%ntrain=10000
    quant%ntest=1000
	quant%sigma=0.1
	quant%lambda=10.0

	nargs = iargc()
	ii=1
	do while(ii<=nargs)
		call getarg(ii,strings)
		if(trim(strings)=='-quant')then ! user-defined quantity parameters   !**** read data file directory, data set dimension, size of training and testing sets, RBF parameters and solver parameters
			flag=1
			do while(flag==1)
				ii=ii+1
				if(ii<=nargs)then
					call getarg(ii,strings)
					if(strings(1:2)=='--')then
						ii=ii+1
						call getarg(ii,strings1)
						if(trim(strings)=='--dimn')then
							read(strings1,*)quant%dimn
						else if	(trim(strings)=='--data_dir')then
							quant%data_dir=trim(strings1)
						else if	(trim(strings)=='--ntrain')then
							read(strings1,*)quant%ntrain
						else if	(trim(strings)=='--ntest')then
							read(strings1,*)quant%ntest
						else if	(trim(strings)=='--sigma')then
							read(strings1,*)quant%sigma
						else if	(trim(strings)=='--lambda')then
							read(strings1,*)quant%lambda
						else
							if(ptree%MyID==Main_ID)write(*,*)'ignoring unknown quant: ', trim(strings)
						endif
					else
						flag=0
					endif
				else
					flag=0
				endif
			enddo
		else if(trim(strings)=='-option')then ! options of ButterflyPACK
			call ReadOption(option,ptree,ii)
		else
			if(ptree%MyID==Main_ID)write(*,*)'ignoring unknown argument: ',trim(strings)
			ii=ii+1
		endif
	enddo

	call PrintOptions(option,ptree)

	quant%trainfile_p=trim(quant%DATA_DIR)//'_train.csv'
	quant%trainfile_l=trim(quant%DATA_DIR)//'_train_label.csv'
	quant%testfile_p=trim(quant%DATA_DIR)//'_test.csv'
	quant%testfile_l=trim(quant%DATA_DIR)//'_test_label.csv'
	quant%Nunk = quant%ntrain
	if(ptree%MyID==Main_ID)write(*,*)'training set: ',trim(quant%trainfile_p)

	t1 = OMP_get_wtime()
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "geometry modeling......"
	open (90,file=quant%trainfile_p)
	allocate (quant%xyz(quant%dimn,1:quant%Nunk))
	do ii=1,quant%Nunk
		read (90,*) quant%xyz(1:quant%dimn,ii)
	enddo
    close(90)
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "modeling finished"
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "    "
	t2 = OMP_get_wtime()

	!**** register the user-defined function and type in ker
	ker%QuantApp => quant
	ker%FuncZmn => Zelem_RBF

    !**** initialization of the construction phase
    allocate(Permutation(quant%Nunk))
	call BPACK_construction_Init(quant%Nunk,Permutation,Nunk_loc,bmat,option,stats,msh,ker,ptree,Coordinates=quant%xyz)
	deallocate(Permutation) ! caller can use this permutation vector if needed

	!**** computation of the construction phase
    call BPACK_construction_Element(bmat,option,stats,msh,ker,ptree)

	!**** factorization phase
    call BPACK_Factorization(bmat,option,stats,ptree,msh)

	!**** solve phase
    call RBF_solve(bmat,option,msh,quant,ptree,stats)

	!**** deletion of quantities
	if(allocated(quant%xyz))deallocate(quant%xyz)
	call delete_proctree(ptree)
	call delete_Hstat(stats)
	call delete_mesh(msh)
	call delete_kernelquant(ker)
	call BPACK_delete(bmat)

    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "-------------------------------program end-------------------------------------"

	call blacs_exit(1)
	call MPI_Finalize(ierr)

end PROGRAM ButterflyPACK_KRR

subroutine RBF_solve(bmat,option,msh,quant,ptree,stats)

    use BPACK_DEFS
	use APPLICATION_MODULE
	use omp_lib
	use BPACK_Solve_Mul

    implicit none

    integer i, j, ii, jj, iii, jjj, ierr, ntest,Dimn,edge_m,edge_n,ncorrect
    integer level, blocks, edge, patch, node, group
    integer rank, index_near, m, n, length, flag, num_sample, n_iter_max, iter, N_unk, N_unk_loc
    real(kind=8) theta, phi, dphi, rcs_V, rcs_H,rate
    real(kind=8) T0,T1
    real(kind=8) n1,n2,rtemp
    real(kind=8) value_Z
    real(kind=8),allocatable:: Voltage_pre(:),x(:,:),b(:,:),vout(:,:),vout_tmp(:,:)
	real(kind=8):: rel_error
	type(Hoption)::option
	type(mesh)::msh
	type(quant_app)::quant
	type(proctree)::ptree
	type(Bmatrix)::bmat
	type(Hstat)::stats
	real(kind=8),allocatable:: current(:),voltage(:)
	real(kind=8), allocatable:: labels(:)
	real(kind=8),allocatable:: xyz_test(:,:)
	real(kind=8) r_mn
	integer label

	if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "Solve and Prediction......"

	N_unk=msh%Nunk
	Dimn=quant%dimn
	N_unk_loc = msh%idxe-msh%idxs+1

	if(option%ErrSol==1)then
		call BPACK_Test_Solve_error(bmat,N_unk_loc,option,ptree,stats)
	endif

	!**** read training label as local RHS and solve for the weights

	allocate(labels(N_unk))
	allocate (x(N_unk_loc,1))
	x=0
	allocate (b(N_unk_loc,1))
	open (91,file=quant%trainfile_l)
	! read(91,*)N_unk,Dimn
	do ii=1,N_unk
		read(91,*)label
		labels(ii)=label
	enddo
	do ii=1,N_unk_loc
		b(ii,1) = labels(msh%new2old(ii-1+msh%idxs))
	enddo
	deallocate(labels)
	close(91)

	n1 = OMP_get_wtime()

	call BPACK_Solution(bmat,x,b,N_unk_loc,1,option,ptree,stats)

	n2 = OMP_get_wtime()
	stats%Time_Sol = stats%Time_Sol + n2-n1
	call MPI_ALLREDUCE(stats%Time_Sol,rtemp,1,MPI_DOUBLE_PRECISION,MPI_MAX,ptree%Comm,ierr)
	if(ptree%MyID==Main_ID .and. option%verbosity>=0)write (*,*) 'Solving:',rtemp,'Seconds'
	call MPI_ALLREDUCE(stats%Flop_Sol,rtemp,1,MPI_DOUBLE_PRECISION,MPI_SUM,ptree%Comm,ierr)
	if(ptree%MyID==Main_ID .and. option%verbosity>=0)write (*,'(A13Es14.2)') 'Solve flops:',rtemp

	!**** prediction on the test sets

	ntest=quant%ntest
	T0 = OMP_get_wtime()
	open (92,file=quant%testfile_p)
	allocate (xyz_test(Dimn,ntest))
	do edge=1,ntest
		read (92,*) xyz_test(1:Dimn,edge)
	enddo
	close(92)

	allocate (vout(ntest,1))
	allocate (vout_tmp(ntest,1))
	vout_tmp = 0
	do edge=1, N_unk_loc
		do edge_m=1,ntest
			r_mn=sum((xyz_test(1:dimn,edge_m)-quant%xyz(1:dimn,msh%new2old(edge+msh%idxs-1)))**2)
			value_Z = exp(-r_mn/2.0/quant%sigma**2)
			vout_tmp(edge_m,1) = vout_tmp(edge_m,1) + value_Z*x(edge,1)
		enddo
	enddo

	call MPI_REDUCE(vout_tmp, vout, ntest,MPI_double_precision, MPI_SUM, Main_ID, ptree%Comm,ierr)
	T1 = OMP_get_wtime()
	if (ptree%MyID==Main_ID) then
		do ii=1,ntest
			if(dble(vout(ii,1))>0)then
				vout(ii,1)=1
			else
				vout(ii,1)=-1
			endif
		enddo

		open (93,file=quant%testfile_l)
		do edge=1,ntest
			read (93,*) label
			vout_tmp(edge,1)=label
		enddo
		close(93)

		ncorrect=0
		do edge=1,ntest
			if(dble(vout_tmp(edge,1))*dble(vout(edge,1))>0)then
				ncorrect = ncorrect + 1
			endif
		enddo

		rate = dble(ncorrect)/dble(ntest)

		write (*,*) ''
		write (*,*) 'Prediction time:',T1-T0,'Seconds'
		write (*,*) 'Success rate:',rate
		write (*,*) ''

	endif

	deallocate (vout)
	deallocate (vout_tmp)

	deallocate(x)
	deallocate(b)
	deallocate(xyz_test)

	call MPI_barrier(ptree%Comm,ierr)

    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "Solve and Prediction finished"
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "    "

    return

end subroutine RBF_solve

max_stars_repo_name: liuyangzhuan/hod-lr-bf

max_issues_repo_name: liuyangzhuan/hod-lr-bf

max_forks_repo_name: liuyangzhuan/hod-lr-bf


text:
    subroutine write_cgns
    implicit none

    !       author: Diane Poirier (diane@icemcfd.com)
    !       last revised on March 8 2000

    !       This example test the complete SIDS for multi-block data.
    !	It creates a dummy mesh composed of 2 structured blocks in 3D.
#include "cgnstypes_f.h"
#ifdef WINNT
    include 'cgnswin_f.h'
#endif
    include 'cgnslib_f.h'

    integer Ndim
    parameter (Ndim = 3)
    cgsize_t one
    parameter (one = 1)

    cgsize_t index_dim, cell_dim, phys_dim
    cgsize_t base_no, zone_no, coord_no, sol_no, discr_no, conn_no
    cgsize_t hole_no, boco_no, field_no, dset_no
    cgsize_t num, NormalIndex(Ndim)
    cgsize_t ndims, npnts,size(Ndim*3)
    cgsize_t cg, ier, zone, coord, i, j, k, n, pos, sol, field
    cgsize_t pnts(Ndim,120), donor_pnts(Ndim,120)
    cgsize_t transform(Ndim)
    cgsize_t nptsets, nrmlistflag
    real*4 data(120), normals(360)
    double precision Dxyz(120), values(120)
    character*32 zonename, solname, fieldname
    character*32 coordname(Ndim)
    character*32 donorname

    coordname(1) = 'CoordinateX'
    coordname(2) = 'CoordinateY'
    coordname(3) = 'CoordinateZ'

    ! *** initialize
    ier = 0
    index_dim=Ndim
    cell_dim=Ndim
    phys_dim=Ndim

    ! *** open CGNS file for writing
    call cg_open_f('cgtest.cgns', CG_MODE_WRITE, cg, ier)
    if (ier .eq. ERROR) call cg_error_exit_f

    ! *** base
    call cg_base_write_f(cg, 'Basename', cell_dim, phys_dim,&
                       base_no, ier)
    if (ier .eq. ERROR) call cg_error_exit_f

    ! *** zone
    do zone=1, 2
        write(zonename,'(a5,i1)') 'zone#',zone
        num = 1
        do i=1,index_dim          		! zone#1: 3*4*5, zone#2: 4*5*6
            size(i) = i+zone+1		! nr of nodes in i,j,k
            size(i+Ndim) = size(i)-1	! nr of elements in i,j,k
            size(i+2*Ndim) = 0		! nr of bnd nodes if ordered
            num = num * size(i)		! nr of nodes
        enddo
        !234567890!234567890!234567890!234567890!234567890!234567890!23456789012
        call cg_zone_write_f(cg, base_no, zonename, size,&
                           Structured, zone_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** coordinate
        do coord=1, phys_dim
            do k=1, size(3)
                do j=1, size(2)
                    do i=1, size(1)
                        pos = i + (j-1)*size(1) + (k-1)*size(1)*size(2)
                        ! * make up some dummy coordinates just for the test:
                        if (coord.eq.1) Dxyz(pos) = i
                        if (coord.eq.2) Dxyz(pos) = j
                        if (coord.eq.3) Dxyz(pos) = k
                    enddo
                enddo
            enddo

            call cg_coord_write_f(cg, base_no, zone_no, RealDouble,&
                            coordname(coord), Dxyz, coord_no, ier)
            if (ier .eq. ERROR) call cg_error_exit_f

        enddo

        ! *** solution
        do sol=1, 2
            write(solname,'(a5,i1,a5,i1)') 'Zone#',zone,' sol#',sol
            call cg_sol_write_f(cg, base_no, zone_no, solname,&
                              Vertex, sol_no, ier)
            if (ier .eq. ERROR) call cg_error_exit_f

            ! *** solution field
            do field=1, 2
                ! make up some dummy solution values
                do i=1, num
                    values(i) = i*field*sol
                enddo
                write(fieldname,'(a6,i1)') 'Field#',field
                call cg_field_write_f(cg, base_no, zone_no, sol_no,&
                  RealDouble, fieldname, values, field_no, ier)
                if (ier .eq. ERROR) call cg_error_exit_f

            enddo				! field loop
        enddo				! solution loop

        ! *** discrete data
        call cg_discrete_write_f(cg, base_no, zone_no, 'discrete#1',&
                               discr_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** discrete data arrays, defined on vertices:
        call cg_goto_f(cg, base_no, ier, 'Zone_t', zone,&
                     'DiscreteData_t', discr_no, 'end')
        if (ier .eq. ERROR) call cg_error_exit_f

        do 123 k=1, size(3)
            do 123 j=1, size(2)
                do 123 i=1, size(1)
                    pos = i + (j-1)*size(1) + (k-1)*size(1)*size(2)
                    data(pos) = pos	! * make up some dummy data
123     continue
        call cg_array_write_f('arrayname', RealSingle, index_dim,&
                             size, data, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** discrete data arrays attribute: GOTO DataArray node
        call cg_goto_f(cg, base_no, ier, 'Zone_t', zone,&
	            'DiscreteData_t', discr_no, 'DataArray_t', 1, 'end')
        if (ier .eq. ERROR) call cg_error_exit_f

        call cg_units_write_f(Kilogram, Meter, Second, Kelvin,&
                            Radian, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** overset holes
        !  create dummy data
        do i=1,3
            ! Define 2 separate PointRange, for 2 patches in the hole
            pnts(i,1)=1
            pnts(i,2)=size(i)
            ! second PointRange of hole
            pnts(i,3)=2
            pnts(i,4)=size(i)
        enddo
        ! Hole defined with 2 point set type PointRange, so 4 points:
        nptsets = 2
        npnts = 4
        call cg_hole_write_f(cg, base_no, zone_no, 'hole#1', Vertex,&
                           PointRange, nptsets, npnts, pnts,&
                           hole_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** general connectivity
        do 100 n=1, 5
            do 100 i=1,3
                pnts(i,n)=i		! * dummy data
                donor_pnts(i,n)=i*2
100     continue
        ! create a point matching connectivity
        npnts = 5
        call cg_conn_write_f(cg, base_no, zone_no, 'Connect#1',&
          Vertex, Abutting1to1, PointList, npnts, pnts, 'zone#2',&
          Structured, PointListDonor, Integer, npnts, donor_pnts,&
          conn_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** connectivity 1to1
        !  generate data
        do i=1,3
            !**make up some dummy data:
            pnts(i,1)=1
            pnts(i,2)=size(i)
            donor_pnts(i,1)=1
            donor_pnts(i,2)=size(i)
            transform(i)=i*(-1)
        enddo
        if (zone .eq. 1) then
            donorname='zone#2'
        else if (zone .eq. 2) then
            donorname='zone#1'
        endif

        call cg_1to1_write_f(cg, base_no, zone_no, '1to1_#1',&
		donorname, pnts, donor_pnts, transform, conn_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** ZoneGridConnectivity attributes:  GOTO ZoneGridConnectivity_t node
        call cg_goto_f(cg, base_no, ier, 'Zone_t', zone,&
                     'ZoneGridConnectivity_t', one, 'end')
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** ZoneGridConnectivity attributes: Descriptor_t
        !234567890!234567890!234567890!234567890!234567890!234567890!23456789012
        call cg_descriptor_write_f('DescriptorName',&
                                 'Zone Connectivity', ier)

        ! *** bocos
        npnts = 2
        call cg_boco_write_f(cg, base_no, zone_no, 'boco#1',&
           BCInflow, PointRange, npnts, pnts, boco_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! *** boco normal
        npnts = 1
        do i=1,Ndim
            NormalIndex(i)=0
            ! compute nr of points on bc patch:
            npnts = npnts * (pnts(i,2)-pnts(i,1)+1)
        enddo
        NormalIndex(1)=1
        do i=1,phys_dim*npnts
            normals(i)=i
        enddo

        nrmlistflag = 1
        call cg_boco_normal_write_f(cg, base_no, zone_no, boco_no,&
         NormalIndex, nrmlistflag, RealSingle, normals, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition attributes: GOTO BC_t node
        call cg_goto_f(cg, base_no, ier, 'Zone_t', zone, 'ZoneBC_t',&
          one, 'BC_t', boco_no, 'end')
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition attributes:  GridLocation_t
        call cg_gridlocation_write_f(Vertex, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition dataset
        call cg_dataset_write_f(cg, base_no, zone,&
         boco_no, 'DataSetName', BCInflow, dset_no, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition data:
        call cg_bcdata_write_f(cg, base_no, zone,&
         boco_no, dset_no, Neumann, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition data arrays: GOTO BCData_t node
        call cg_goto_f(cg, base_no, ier, 'Zone_t', zone_no,&
          'ZoneBC_t', one, 'BC_t', boco_no, 'BCDataSet_t',&
          dset_no, 'BCData_t', Neumann, 'end')
        if (ier .eq. ERROR) call cg_error_exit_f

        do i=1, npnts
            data(i) = i
        enddo
        ndims = 1
        call cg_array_write_f('dataset_arrayname', RealSingle,&
           ndims, npnts, data, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

        ! ** boundary condition data attributes:
        call cg_dataclass_write_f(NormalizedByDimensional, ier)
        if (ier .eq. ERROR) call cg_error_exit_f

    enddo					! zone loop

    ! *** close CGNS file
    call cg_close_f(cg, ier)
    if (ier .eq. ERROR) call cg_error_exit_f

    end subroutine write_cgns

max_stars_repo_name: nescirem/nesDockerfile

max_issues_repo_name: nescirem/nesDockerfile

max_forks_repo_name: nescirem/nesDockerfile


text:
PROGRAM meulera

! Code converted using TO_F90 by Alan Miller
! Date: 2004-06-28  Time: 12:58:11

!       ==========================================================
!       Purpose: This program computes Euler number En using
!                subroutine EULERA
!       Example: Compute Euler number En for n = 0,2,...,10
!                Computed results:

!                   n            En
!                 --------------------------
!                   0     .100000000000D+01
!                   2    -.100000000000D+01
!                   4     .500000000000D+01
!                   6    -.610000000000D+02
!                   8     .138500000000D+04
!                  10    -.505210000000D+05
!       ==========================================================

DOUBLE PRECISION :: e
DIMENSION e(0:200)
WRITE(*,*)'  Please enter Nmax '
!        READ(*,*)N
n=10
CALL eulera(n,e)
WRITE(*,*)'   n            En'
WRITE(*,*)' --------------------------'
DO  k=0,n,2
  WRITE(*,20)k,e(k)
END DO
20      FORMAT(2X,i3,d22.12)
END PROGRAM meulera

SUBROUTINE eulera(n,en)

!       ======================================
!       Purpose: Compute Euler number En
!       Input :  n --- Serial number
!       Output:  EN(n) --- En
!       ======================================

INTEGER, INTENT(IN)                      :: n
DOUBLE PRECISION, INTENT(OUT)            :: en(0:n)
IMPLICIT DOUBLE PRECISION (a-h,o-z)

en(0)=1.0D0
DO  m=1,n/2
  s=1.0D0
  DO  k=1,m-1
    r=1.0D0
    DO  j=1,2*k
      r=r*(2.0D0*m-2.0D0*k+j)/j
    END DO
    s=s+r*en(2*k)
  END DO
  en(2*m)=-s
END DO
RETURN
END SUBROUTINE eulera

max_stars_repo_name: vasaantk/bin

max_issues_repo_name: vasaantk/bin

max_forks_repo_name: vasaantk/bin


text:
c
c  L-BFGS-B is released under the “New BSD License” (aka “Modified BSD
c  License” or “3-clause license”)
c  Please read attached file License.txt
c
      subroutine dpofa(a,lda,n,info)
      integer lda,n,info
      double precision a(lda,*)
c
c     dpofa factors a double precision symmetric positive definite
c     matrix.
c
c     dpofa is usually called by dpoco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dpoco) = (1 + 18/n)*(time for dpofa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the symmetric matrix to be factored.  only the
c                diagonal and upper triangle are used.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix  r  so that  a = trans(r)*r
c                where  trans(r)  is the transpose.
c                the strict lower triangle is unaltered.
c                if  info .ne. 0 , the factorization is not complete.
c
c        info    integer
c                = 0  for normal return.
c                = k  signals an error condition.  the leading minor
c                     of order  k  is not positive definite.
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas ddot
c     fortran sqrt
c
c     internal variables
c
      double precision ddot,t
      double precision s
      integer j,jm1,k
c     begin block with ...exits to 40
c
c
         do 30 j = 1, n
            info = j
            s = 0.0d0
            jm1 = j - 1
            if (jm1 .lt. 1) go to 20
            do 10 k = 1, jm1
               t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
               t = t/a(k,k)
               a(k,j) = t
               s = s + t*t
   10       continue
   20       continue
            s = a(j,j) - s
c     ......exit
            if (s .le. 0.0d0) go to 40
            a(j,j) = sqrt(s)
   30    continue
         info = 0
   40 continue
      return
      end

c====================== The end of dpofa ===============================

      subroutine dtrsl(t,ldt,n,b,job,info)
      integer ldt,n,job,info
      double precision t(ldt,*),b(*)
c
c
c     dtrsl solves systems of the form
c
c                   t * x = b
c     or
c                   trans(t) * x = b
c
c     where t is a triangular matrix of order n. here trans(t)
c     denotes the transpose of the matrix t.
c
c     on entry
c
c         t         double precision(ldt,n)
c                   t contains the matrix of the system. the zero
c                   elements of the matrix are not referenced, and
c                   the corresponding elements of the array can be
c                   used to store other information.
c
c         ldt       integer
c                   ldt is the leading dimension of the array t.
c
c         n         integer
c                   n is the order of the system.
c
c         b         double precision(n).
c                   b contains the right hand side of the system.
c
c         job       integer
c                   job specifies what kind of system is to be solved.
c                   if job is
c
c                        00   solve t*x=b, t lower triangular,
c                        01   solve t*x=b, t upper triangular,
c                        10   solve trans(t)*x=b, t lower triangular,
c                        11   solve trans(t)*x=b, t upper triangular.
c
c     on return
c
c         b         b contains the solution, if info .eq. 0.
c                   otherwise b is unaltered.
c
c         info      integer
c                   info contains zero if the system is nonsingular.
c                   otherwise info contains the index of
c                   the first zero diagonal element of t.
c
c     linpack. this version dated 08/14/78 .
c     g. w. stewart, university of maryland, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c     fortran mod
c
c     internal variables
c
      double precision ddot,temp
      integer case,j,jj
c
c     begin block permitting ...exits to 150
c
c        check for zero diagonal elements.
c
         do 10 info = 1, n
c     ......exit
            if (t(info,info) .eq. 0.0d0) go to 150
   10    continue
         info = 0
c
c        determine the task and go to it.
c
         case = 1
         if (mod(job,10) .ne. 0) case = 2
         if (mod(job,100)/10 .ne. 0) case = case + 2
         go to (20,50,80,110), case
c
c        solve t*x=b for t lower triangular
c
   20    continue
            b(1) = b(1)/t(1,1)
            if (n .lt. 2) go to 40
            do 30 j = 2, n
               temp = -b(j-1)
               call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1)
               b(j) = b(j)/t(j,j)
   30       continue
   40       continue
         go to 140
c
c        solve t*x=b for t upper triangular.
c
   50    continue
            b(n) = b(n)/t(n,n)
            if (n .lt. 2) go to 70
            do 60 jj = 2, n
               j = n - jj + 1
               temp = -b(j+1)
               call daxpy(j,temp,t(1,j+1),1,b(1),1)
               b(j) = b(j)/t(j,j)
   60       continue
   70       continue
         go to 140
c
c        solve trans(t)*x=b for t lower triangular.
c
   80    continue
            b(n) = b(n)/t(n,n)
            if (n .lt. 2) go to 100
            do 90 jj = 2, n
               j = n - jj + 1
               b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1)
               b(j) = b(j)/t(j,j)
   90       continue
  100       continue
         go to 140
c
c        solve trans(t)*x=b for t upper triangular.
c
  110    continue
            b(1) = b(1)/t(1,1)
            if (n .lt. 2) go to 130
            do 120 j = 2, n
               b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1)
               b(j) = b(j)/t(j,j)
  120       continue
  130       continue
  140    continue
  150 continue
      return
      end

c====================== The end of dtrsl ===============================

max_stars_repo_name: datagangchen/RLSTL


text:
      SUBROUTINE PERMIT(PERIM,IGO)
*
*
*       Check on existing multiple regularization.
*       ------------------------------------------
*
      INCLUDE 'common6.h'
      PARAMETER  (NMX=10,NMX3=3*NMX,NMXm=NMX*(NMX-1)/2)
      REAL*8  M,MASS,MC,MIJ,MKK
      COMMON/CLUMP/   BODYS(NCMAX,5),T0S(5),TS(5),STEPS(5),RMAXS(5),
     &                NAMES(NCMAX,5),ISYS(5)
      COMMON/CHAIN1/  XCH(NMX3),VCH(NMX3),M(NMX),
     &                ZZ(NMX3),WC(NMX3),MC(NMX),
     &                XI(NMX3),PI(NMX3),MASS,RINV(NMXm),RSUM,MKK(NMX),
     &                MIJ(NMX,NMX),TKK(NMX),TK1(NMX),INAME(NMX),NN
*
*
*       Search any existing subsystem.
      ISUB = 0
      ICHSUB = 1
      DO 10 L = 1,NSUB
*       Identify chain pointer for possible reduction of STEPS.
          IF (ISYS(L).EQ.3) ICHSUB = ISYS(L)
*       Distinguish between triple & quad case (denoted ISUB = 1 or 2).
          IF (JCOMP.LE.N.AND.NAMES(4,L).EQ.0) THEN
              ISUB = 1
      ELSE IF (JCOMP.GT.N.AND.NAMES(4,L).GT.0) THEN
              ISUB = 2
          END IF
   10 CONTINUE
*
*       Do not allow a second regularization of the same type.
      IF (ISUB.GT.0) THEN
*       See whether the case ISUB = 1 or 2 is used already.
          DO 20 L = 1,NSUB
              IF (ISUB.EQ.ISYS(L)) IGO = 1
   20     CONTINUE
*       Enforce chain termination at next extension if new system < RSUM/2.
          IF (PERIM.LT.0.5*RSUM.AND.IGO.GT.0) THEN
              STEPS(ICHSUB) = 0.0
          END IF
      END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
       subroutine tdscal(naxes,scale)
C      subroutine tdscal  01.09.78
C  scaling at the direction of axe (if naxe=4 then at all directions)
         common/gfmatr/t(4,4),a(4,4),nt,nt1
      nt=2
       scal1=scale
       if(naxes.ge.4) scal1=1./scale
       call hcunit(a)
       a(naxes,naxes)=scal1
       call hcmult(t,a)
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


subroutine driver
    use dynamo
    implicit none
    character( len = 256 ) :: str
    integer :: i, j

    call getarg( 1, str )
    write(*,"(a,a,a)") "[", trim( str ), "]"
    open( file = trim( str ), unit = 998, action = "read", form = "formatted" )
    read( 998, "(a)" ) str
    write(*,"(a,a,a)") "[", trim( str ), "]"
    do while( trim( str ) /= "exit" )
        if( trim( str ) == "charges" ) then
            open( file = "dynamo.chrg", unit = 999, action = "read", form = "unformatted" )
            read( 999 ) atmchg
            close( 999 )
        end if
        if( trim( str ) == "coordinates" ) then
            open( file = "dynamo.crd", unit = 999, action = "read", form = "unformatted" )
            read( 999 ) atmcrd(1:3,1:natoms)
            close( 999 )
        end if
        if( trim( str ) == "energy" ) then
            call energy
            open( file = "dynamo.dat", unit = 999, action = "write", form = "unformatted" )
            write( 999 ) etotal
            close( 999 )
        end if
        if( trim( str ) == "gradient" ) then
            if( .not. allocated( atmder ) ) allocate( atmder(1:3,1:natoms ) )
            call gradient
            open( file = "dynamo.dat", unit = 999, action = "write", form = "unformatted" )
            write( 999 ) etotal
            write( 999 ) atmder
            close( 999 )
        end if
        read( 998, "(a)" ) str
        write(*,"(a,a,a)") "[", trim( str ), "]"
    end do
    close( 998 )
end subroutine driver

program slave
    use dynamo
    implicit none

    logical, dimension(:), allocatable :: flg
    integer :: i

    call dynamo_header

    call mm_file_process( "borra", "opls" )
    call mm_system_construct( "borra", "seq" )
    call coordinates_read( "crd" )

    allocate( flg(1:natoms) )
    flg = .false.
    flg = atom_selection( subsystem = (/ "W" /), residue_number = (/ 1 /) )
    call atoms_fix( flg )
    do i = 1, natoms
        if( flg(i) ) then
            atmchg(i) = 0.0_dp
            atmchg14(i) = 0.0_dp
        end if
    end do

    call energy_initialize
    call energy_non_bonding_options( &
        list_cutoff   = 18.0_dp, &
        outer_cutoff  = 16.0_dp, &
        inner_cutoff  = 14.0_dp, &
        minimum_image = .false. )

    call driver

    deallocate( flg )
end program


text:
version https://git-lfs.github.com/spec/v1
oid sha256:4f6e1d35a72e0fea2b0fcd2db89eb74c2d834d8ef6aee36c2547f4842cafa40c
size 44819

max_stars_repo_name: initialz/tensorflow-yolo-face-ios

max_issues_repo_name: initialz/tensorflow-yolo-face-ios

max_forks_repo_name: initialz/tensorflow-yolo-face-ios


text:
C ================================================================
C  All operations with colors are binary operations.
C  The same operations go with the input parameter status.
C ================================================================

C ================================================================
C     Logical Function ifXstatus(status, iXstatus)
      Logical Function ifXnode(clr, iXnode)
C ================================================================
      Integer clr, iXnode

      ifXnode = IAND(clr, iXnode) .EQ. iXnode

      Return
      End

C ================================================================
C     Subroutine addXstatus(status, iXstatus)
      Subroutine addXnode(clr, iXnode)
C ================================================================
      Integer clr, iXnode

      clr = IOR(clr, iXnode)

      Return
      End

C ================================================================
C     Subroutine delXstatus(status, iXstatus)
      Subroutine delXnode(clr, iXnode)
C ================================================================
      Integer clr, iXnode

      clr = clr - IAND(clr, iXnode)

      Return
      End

C ================================================================
      Integer Function minClr(clr1, clr2)
C ================================================================
C  The function returns common color for both clr1 and clr2.
C ================================================================
      Integer clr1, clr2

      minClr = IAND(clr1, clr2)

      Return
      End

C ================================================================
      Integer Function maxClr(clr1, clr2)
C ================================================================
C Routine returns minimal color containing both clr1 and clr2.
C ================================================================
      Integer clr1, clr2

      maxClr = IOR(clr1, clr2)

      Return
      End

C ================================================================
      Subroutine setStatus(flagAuto, status, iPrint)
C ================================================================
      include 'status.fd'
C ================================================================
C Routine adds additional properties to the variable status
C ================================================================
      Logical  flagAuto
      Integer  status, iPrint

      Logical  ifXnode

C ================================================================
      status = max(0, status)

c ... remove obsolete and not-implemented input features
      If(status.GT.0) Then
         If(ifXnode(status, ANISmoothMesh)) Then
            Call delXnode(status, ANISmoothMesh)
            If(iPrint.GE.1) Write(*, 5001)
         End if

         If(ifXnode(status, ANIMultiConnectedGeometry)) Then
           Call delXnode(status, ANIMultiConnectedGeometry)
            If(iPrint.GE.1) Write(*, 5002)
         End if
      End if

c ... inform the user about requested features
      If(iPrint.GE.1) Then
         If(ifXnode(status, ANIForbidBoundaryElements)) Write(*, 5003)
         If(ifXnode(status, ANIFixBoundaryEdges))       Write(*, 5004)
         If(ifXnode(status, ANIFixBoundaryPoints))      Write(*, 5007)

         If(ifXnode(status, ANIUse2ArmRule))
     &      Write(*, 5005) '[user]'
         If(ifXnode(status, ANIDeleteTemporaryEdges))
     &      Write(*, 5006) '[user]'
         If(ifXnode(status, ANIUntangleMesh))
     &      Write(*, 5008) '[user]'
      End if

c ... set up default features
      If(flagAuto) Then
         If(.NOT.ifXnode(status, ANIUse2ArmRule)) Then
            If(iPrint.GE.2) Write(*, 5005) '[system]'

            Call addXnode(status, ANIUse2ArmRule)
         End if

         If(.NOT.ifXnode(status, ANIDeleteTemporaryEdges)) Then
            If(iPrint.GE.2) Write(*, 5006) '[system]'

            Call addXnode(status, ANIDeleteTemporaryEdges)
         End if
      End if

      If(iPrint.GE.1) Write(*,*)

      Return

 5001 Format('status.fd: -256  [ANISmoothMesh]             [obsolete]')
 5002 Format('status.fd: -64   [ANIMultiConnectedGeometry] [not used]')
 5003 Format('status.fd: +1    [ANIForbidBoundaryElements] [user]')
 5004 Format('status.fd: +4    [ANIFixBoundaryEdges]       [user]')
 5005 Format('status.fd: +2    [ANIUse2ArmRule]            ', A)
 5006 Format('status.fd: +8    [ANIDeleteTemporaryEdges]   ', A)
 5007 Format('status.fd: +16   [ANIFixBoundaryPoints]      [user]')
 5008 Format('status.fd: +512  [ANIUntangleMesh]           ', A)

      Return
      End


text:
! Compute the source term to the divergence constraint, S.
!
! We also construct the delta_gamma1_term here which is used to
! correct for the approximation that Gamma_1 -> \bar{Gamma_1} in the
! constraint.  First we compute the main part of the term and then
! we average it.  Finally we correct the full version of the term to
! account for the psi term that comes from the base state.

module make_S_cc_module

  use bl_types
  use multifab_module

  implicit none

  private

  public :: make_S_cc

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine make_S_cc(S_cc,delta_gamma1_term,delta_gamma1, &
                       state,u, &
                       normal, &
                       rho_omegadot,rho_Hnuc,rho_Hext,thermal, &
                       p0,gamma1bar,delta_gamma1_termbar,psi, &
                       dx,mla,the_bc_level)

    use bl_constants_module
    use bl_prof_module
    use probin_module, only: use_delta_gamma1_term
    use ml_layout_module
    use average_module
    use ml_restrict_fill_module
    use variables, only: foextrap_comp
    use geometry, only: spherical, nr_fine, dr
    use fill_3d_module, only : put_1d_array_on_cart

    type(multifab) , intent(inout) :: S_cc(:)
    type(multifab) , intent(inout) :: delta_gamma1_term(:)
    type(multifab) , intent(inout) :: delta_gamma1(:)
    type(multifab) , intent(in   ) :: state(:)
    type(multifab) , intent(in   ) :: u(:)
    type(multifab) , intent(in   ) :: normal(:)
    type(multifab) , intent(in   ) :: rho_omegadot(:)
    type(multifab) , intent(in   ) :: rho_Hnuc(:)
    type(multifab) , intent(in   ) :: rho_Hext(:)
    type(multifab) , intent(in   ) :: thermal(:)
    real(kind=dp_t), intent(in   ) :: p0(:,0:)
    real(kind=dp_t), intent(in   ) :: gamma1bar(:,0:)
    real(kind=dp_t), intent(inout) :: delta_gamma1_termbar(:,0:)
    real(kind=dp_t), intent(in   ) :: psi(:,0:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(ml_layout), intent(in   ) :: mla
    type(bc_level) , intent(in   ) :: the_bc_level(:)

    real(kind=dp_t), pointer:: srcp(:,:,:,:),dgtp(:,:,:,:),sp(:,:,:,:),up(:,:,:,:)
    real(kind=dp_t), pointer:: tp(:,:,:,:),dgp(:,:,:,:)
    real(kind=dp_t), pointer:: omegap(:,:,:,:), hep(:,:,:,:), hnp(:,:,:,:)
    real(kind=dp_t), pointer:: np(:,:,:,:)
    real(kind=dp_t), pointer:: gp0p(:,:,:,:), p0p(:,:,:,:), g1p(:,:,:,:), psp(:,:,:,:)

    integer :: lo(MAX_SPACEDIM),hi(MAX_SPACEDIM),dm,nlevs
    integer :: i,n, r
    integer :: ng_sr,ng_dt,ng_dg,ng_s,ng_u,ng_rw,ng_he,ng_hn,ng_th
    integer :: ng_gp, ng_p0, ng_g1, ng_n, ng_ps

    real(kind=dp_t) :: gradp0_rad(1,0:nr_fine-1)
    type(multifab) :: gradp0_cart(mla%nlevel)
    type(multifab) :: p0_cart(mla%nlevel)
    type(multifab) :: psi_cart(mla%nlevel)
    type(multifab) :: gamma1bar_cart(mla%nlevel)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "make_S_cc")

    dm = mla%dim
    nlevs = mla%nlevel

    if (spherical == 1) then
       do n = 1, nlevs
          call build(gradp0_cart(n), mla%la(n), 1, 0)
          call setval(gradp0_cart(n), ZERO, all=.true.)

          call build(p0_cart(n), mla%la(n), 1, 0)
          call setval(p0_cart(n), ZERO, all=.true.)

          call build(psi_cart(n), mla%la(n), 1, 0)
          call setval(psi_cart(n), ZERO, all=.true.)

          call build(gamma1bar_cart(n), mla%la(n), 1, 0)
          call setval(gamma1bar_cart(n), ZERO, all=.true.)
       enddo

       if (use_delta_gamma1_term) then
          ! compute gradp0 and put it on a cart
          do r = 0, nr_fine-1
             if (r == 0) then
                gradp0_rad(1,r) = (p0(1,r+1) - p0(1,r))/dr(1)
             else if (r == nr_fine-1) then
                gradp0_rad(1,r) = (p0(1,r) - p0(1,r-1))/dr(1)
             else
                gradp0_rad(1,r) = HALF*(p0(1,r+1) - p0(1,r-1))/dr(1)
             endif
          enddo

          call put_1d_array_on_cart(gradp0_rad,gradp0_cart,foextrap_comp, &
                                    .false.,.false.,dx,the_bc_level,mla)

          call put_1d_array_on_cart(p0,p0_cart,foextrap_comp, &
                                    .false.,.false.,dx,the_bc_level,mla)

          call put_1d_array_on_cart(psi,psi_cart,foextrap_comp, &
                                    .false.,.false.,dx,the_bc_level,mla)

          call put_1d_array_on_cart(gamma1bar,gamma1bar_cart,foextrap_comp, &
                                    .false.,.false.,dx,the_bc_level,mla)
       endif
    endif

    ng_sr = nghost(S_cc(1))
    ng_dt = nghost(delta_gamma1_term(1))
    ng_dg = nghost(delta_gamma1(1))
    ng_s  = nghost(state(1))
    ng_u  = nghost(u(1))
    ng_rw = nghost(rho_omegadot(1))
    ng_he = nghost(rho_Hext(1))
    ng_hn = nghost(rho_Hnuc(1))
    ng_th = nghost(thermal(1))

    do n = 1, nlevs

       lo(:) = 1; hi(:) = 1

       do i = 1, nfabs(state(n))
          srcp => dataptr(S_cc(n), i)
          dgtp     => dataptr(delta_gamma1_term(n), i)
          dgp    => dataptr(delta_gamma1(n), i)
          sp     => dataptr(state(n), i)
          up     => dataptr(u(n), i)
          omegap => dataptr(rho_omegadot(n), i)
          hep    => dataptr(rho_Hext(n), i)
          hnp    => dataptr(rho_Hnuc(n), i)
          tp     => dataptr(thermal(n), i)
          lo(1:dm) = lwb(get_box(state(n), i))
          hi(1:dm) = upb(get_box(state(n), i))

          if (spherical == 1) then

             gp0p => dataptr(gradp0_cart(n), i)
             ng_gp = nghost(gradp0_cart(1))

             p0p => dataptr(p0_cart(n), i)
             ng_p0 = nghost(p0_cart(1))

             g1p => dataptr(gamma1bar_cart(n), i)
             ng_g1 = nghost(gamma1bar_cart(1))

             np => dataptr(normal(n), i)
             ng_n = nghost(normal(1))

             call make_S_cc_3d_sphr(lo, hi, srcp(:,:,:,1), ng_sr, dgtp(:,:,:,1), ng_dt, &
                                    dgp(:,:,:,1), ng_dg, &
                                    sp(:,:,:,:), ng_s, up(:,:,:,:), ng_u, &
                                    omegap(:,:,:,:), ng_rw, hnp(:,:,:,1), ng_hn, &
                                    hep(:,:,:,1), ng_he, &
                                    tp(:,:,:,1), ng_th, &
                                    gp0p(:,:,:,1), ng_gp, &
                                    p0p(:,:,:,1), ng_p0, &
                                    g1p(:,:,:,1), ng_g1, &
                                    np(:,:,:,:), ng_n)
          else
             call make_S_cc_cart(dm, n, lo, hi, &
                                 srcp(:,:,:,1), lbound(srcp), &
                                 dgtp(:,:,:,1), lbound(dgtp), &
                                 dgp(:,:,:,1), lbound(dgp), &
                                 sp, lbound(sp), &
                                 up, lbound(up), &
                                 omegap, lbound(omegap), &
                                 hnp(:,:,:,1), lbound(hnp), &
                                 hep(:,:,:,1), lbound(hep), &
                                 tp(:,:,:,1), lbound(tp), &
                                 p0(n,:), gamma1bar(n,:), dx(n,:))
          end if

       end do
    enddo

    ! restrict data (has no ghost cells)
    call ml_restrict_and_fill(nlevs,S_cc,mla%mba%rr,the_bc_level, &
                              icomp=1, &
                              bcomp=foextrap_comp, &
                              nc=1, &
                              ng=S_cc(1)%ng)

    ! restrict data (has no ghost cells)
    call ml_restrict_and_fill(nlevs,delta_gamma1_term,mla%mba%rr,the_bc_level, &
                              icomp=1, &
                              bcomp=foextrap_comp, &
                              nc=1, &
                              ng=delta_gamma1_term(1)%ng)

    if (use_delta_gamma1_term) then

       call average(mla,delta_gamma1_term,delta_gamma1_termbar,dx,1)

       do n = 1, nlevs
          do i = 1, nfabs(state(n))
             dgtp   => dataptr(delta_gamma1_term(n), i)
             dgp    => dataptr(delta_gamma1(n), i)

             lo(1:dm) = lwb(get_box(state(n), i))
             hi(1:dm) = upb(get_box(state(n), i))

             select case (dm)
             case (1)
                call correct_delta_gamma1_term_1d(lo,hi,dgtp(:,1,1,1),ng_dt, &
                                                  dgp(:,1,1,1),ng_dg, &
                                                  gamma1bar(n,:),psi(n,:),p0(n,:))
             case (2)
                call correct_delta_gamma1_term_2d(lo,hi,dgtp(:,:,1,1),ng_dt, &
                                                  dgp(:,:,1,1),ng_dg, &
                                                  gamma1bar(n,:),psi(n,:),p0(n,:))
             case (3)
                if (spherical == 1) then
                   g1p => dataptr(gamma1bar_cart(n), i)
                   ng_g1 = nghost(gamma1bar_cart(1))

                   p0p => dataptr(p0_cart(n), i)
                   ng_p0 = nghost(p0_cart(1))

                   psp => dataptr(psi_cart(n), i)
                   ng_ps = nghost(psi_cart(1))

                   call correct_delta_gamma1_term_3d_sphr(lo,hi,dgtp(:,:,:,1),ng_dt, &
                                                          dgp(:,:,:,1),ng_dg, &
                                                          g1p(:,:,:,1),ng_g1, &
                                                          psp(:,:,:,1),ng_ps, &
                                                          p0p(:,:,:,1),ng_p0)
                else
                   call correct_delta_gamma1_term_3d(lo,hi,dgtp(:,:,:,1),ng_dt, &
                                                     dgp(:,:,:,1),ng_dg, &
                                                     gamma1bar(n,:),psi(n,:),p0(n,:))
                end if
             end select
          end do
       enddo

    end if

    if (spherical == 1) then
       do n = 1, nlevs
          call destroy(gradp0_cart(n))
          call destroy(p0_cart(n))
          call destroy(psi_cart(n))
          call destroy(gamma1bar_cart(n))
       enddo
    endif

    call destroy(bpt)

  end subroutine make_S_cc

  subroutine make_S_cc_cart(dm, n, lo, hi, &
                            S_cc, srlo, &
                            delta_gamma1_term, dtlo, &
                            delta_gamma1, dglo, &
                            s, slo, &
                            u, ulo, &
                            rho_omegadot, rwlo, &
                            rho_Hnuc, hnlo, &
                            rho_Hext, helo, &
                            thermal, thlo, &
                            p0,gamma1bar,dx)

    use bl_constants_module
    use eos_module, only: eos, eos_input_rt
    use eos_type_module
    use network, only: nspec
    use variables, only: rho_comp, temp_comp, spec_comp
    use probin_module, only: use_delta_gamma1_term
    use geometry, only: anelastic_cutoff_coord, nr

    integer         , intent(in   ) :: dm, n, lo(:),hi(:)
    integer         , intent(in   ) :: srlo(4), dtlo(4), dglo(4)
    integer         , intent(in   ) :: slo(4), ulo(4), rwlo(4)
    integer         , intent(in   ) :: hnlo(4), helo(4), thlo(4)

    real (kind=dp_t), intent(  out) ::              S_cc(srlo(1): ,srlo(2): ,srlo(3): )
    real (kind=dp_t), intent(  out) :: delta_gamma1_term(dtlo(1): ,dtlo(2): ,dtlo(3): )
    real (kind=dp_t), intent(  out) ::      delta_gamma1(dglo(1): ,dglo(2): ,dglo(3): )
    real (kind=dp_t), intent(in   ) ::                 s( slo(1): , slo(2): , slo(3): , slo(4): )
    real (kind=dp_t), intent(in   ) ::                 u( ulo(1): , ulo(2): , ulo(3): , ulo(4): )
    real (kind=dp_t), intent(in   ) ::      rho_omegadot(rwlo(1): ,rwlo(2): ,rwlo(3): ,rwlo(4): )
    real (kind=dp_t), intent(in   ) ::          rho_Hnuc(hnlo(1): ,hnlo(2): ,hnlo(3): )
    real (kind=dp_t), intent(in   ) ::          rho_Hext(helo(1): ,helo(2): ,helo(3): )
    real (kind=dp_t), intent(in   ) ::           thermal(thlo(1): ,thlo(2): ,thlo(3): )
    real (kind=dp_t), intent(in   ) :: p0(0:)
    real (kind=dp_t), intent(in   ) :: gamma1bar(0:)
    real (kind=dp_t), intent(in   ) :: dx(:)

    !     Local variables
    integer         :: i, j, k, r, comp
    real(kind=dp_t) :: sigma, xi_term, pres_term, gradp0

    integer :: pt_index(MAX_SPACEDIM)
    type(eos_t) :: eos_state

    S_cc = zero
    r = -1

    !$OMP PARALLEL DO PRIVATE(i,j,k,r,comp,sigma,xi_term,pres_term,gradp0,eos_state,pt_index)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             eos_state%rho   = s(i,j,k,rho_comp)
             eos_state%T     = s(i,j,k,temp_comp)
             eos_state%xn(:) = s(i,j,k,spec_comp:spec_comp+nspec-1)/eos_state%rho

             pt_index(:) = (/i, j, k/)

             ! dens, temp, and xmass are inputs
             call eos(eos_input_rt, eos_state, pt_index)

             sigma = eos_state%dpdt / &
                  (eos_state%rho * eos_state%cp * eos_state%dpdr)

             xi_term = ZERO
             pres_term = ZERO
             do comp = 1, nspec
                xi_term = xi_term - &
                     eos_state%dhdX(comp)*rho_omegadot(i,j,k,comp)/eos_state%rho

                pres_term = pres_term + &
                     eos_state%dpdX(comp)*rho_omegadot(i,j,k,comp)/eos_state%rho
             enddo

             S_cc(i,j,k) = (sigma/eos_state%rho) * &
                  ( rho_Hext(i,j,k) + rho_Hnuc(i,j,k) + thermal(i,j,k) ) &
                  + sigma*xi_term &
                  + pres_term/(eos_state%rho*eos_state%dpdr)

             select case (dm)
             case (1)
                r = i
             case (2)
                r = j
             case (3)
                r = k
             end select

             if (use_delta_gamma1_term .and. r < anelastic_cutoff_coord(n)) then
                if (r .eq. 0) then
                   gradp0 = (p0(r+1) - p0(r))/dx(dm)
                else if (r .eq. nr(n)-1) then
                   gradp0 = (p0(r) - p0(r-1))/dx(dm)
                else
                   gradp0 = HALF*(p0(r+1) - p0(r-1))/dx(dm)
                endif

                delta_gamma1(i,j,k) = eos_state%gam1 - gamma1bar(r)

                delta_gamma1_term(i,j,k) = (eos_state%gam1 - gamma1bar(r))*u(i,j,k,dm)* &
                     gradp0/(gamma1bar(r)*gamma1bar(r)*p0(r))
             else
                delta_gamma1_term(i,j,k) = ZERO
                delta_gamma1(i,j,k) = ZERO
             endif

          enddo
       enddo
    enddo
    !$OMP END PARALLEL DO

  end subroutine make_S_cc_cart

  subroutine make_S_cc_3d_sphr(lo,hi,S_cc,ng_sr,dg1_term,ng_dt,delta_gamma1, &
                               ng_dg,s,ng_s,u,ng_u,rho_omegadot,ng_rw,rho_Hnuc,ng_hn, &
                               rho_Hext,ng_he,thermal,ng_th, &
                               gradp0_cart,ng_gp,p0_cart,ng_p0, &
                               gamma1bar_cart,ng_g1,normal,ng_n)

    use bl_constants_module
    use eos_module, only: eos, eos_input_rt
    use eos_type_module
    use network, only: nspec
    use variables, only: rho_comp, temp_comp, spec_comp
    use probin_module, only: use_delta_gamma1_term

    integer         , intent(in   ) :: lo(:),hi(:)
    integer         , intent(in   ) :: ng_sr,ng_dt,ng_dg,ng_s,ng_u, &
                                       ng_rw,ng_he, &
                                       ng_hn,ng_th,ng_gp,ng_p0,ng_g1,ng_n

    real (kind=dp_t), intent(  out) ::           S_cc(lo(1)-ng_sr:,lo(2)-ng_sr:,lo(3)-ng_sr:)
    real (kind=dp_t), intent(  out) ::       dg1_term(lo(1)-ng_dt:,lo(2)-ng_dt:,lo(3)-ng_dt:)
    real (kind=dp_t), intent(  out) ::   delta_gamma1(lo(1)-ng_dg:,lo(2)-ng_dg:,lo(3)-ng_dg:)
    real (kind=dp_t), intent(in   ) ::              s(lo(1)-ng_s :,lo(2)-ng_s :,lo(3)-ng_s :,:)
    real (kind=dp_t), intent(in   ) ::              u(lo(1)-ng_s :,lo(2)-ng_s :,lo(3)-ng_s :,:)
    real (kind=dp_t), intent(in   ) ::   rho_omegadot(lo(1)-ng_rw:,lo(2)-ng_rw:,lo(3)-ng_rw:,:)
    real (kind=dp_t), intent(in   ) ::       rho_Hnuc(lo(1)-ng_hn:,lo(2)-ng_hn:,lo(3)-ng_hn:)
    real (kind=dp_t), intent(in   ) ::       rho_Hext(lo(1)-ng_he:,lo(2)-ng_he:,lo(3)-ng_he:)
    real (kind=dp_t), intent(in   ) ::        thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real (kind=dp_t), intent(in   ) ::    gradp0_cart(lo(1)-ng_gp:,lo(2)-ng_gp:,lo(3)-ng_gp:)
    real (kind=dp_t), intent(in   ) ::        p0_cart(lo(1)-ng_p0:,lo(2)-ng_p0:,lo(3)-ng_p0:)
    real (kind=dp_t), intent(in   ) :: gamma1bar_cart(lo(1)-ng_g1:,lo(2)-ng_g1:,lo(3)-ng_g1:)
    real (kind=dp_t), intent(in   ) ::         normal(lo(1)-ng_n:,lo(2)-ng_n:,lo(3)-ng_n:,:)

    !     Local variables
    integer         :: i, j, k, comp
    real(kind=dp_t) :: sigma, xi_term, pres_term

    integer :: pt_index(MAX_SPACEDIM)
    type(eos_t) :: eos_state

    real(kind=dp_t) :: Ut_dot_er

    S_cc = zero

    !$OMP PARALLEL DO PRIVATE(i,j,k,comp,sigma,xi_term,pres_term,eos_state,pt_index,Ut_dot_er)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             eos_state%rho   = s(i,j,k,rho_comp)
             eos_state%T     = s(i,j,k,temp_comp)
             eos_state%xn(:) = s(i,j,k,spec_comp:spec_comp+nspec-1)/eos_state%rho

             pt_index(:) = (/i, j, k/)

             ! dens, temp, and xmass are inputs
             call eos(eos_input_rt, eos_state, pt_index)

             sigma = eos_state%dpdt / &
                  (eos_state%rho * eos_state%cp * eos_state%dpdr)

             xi_term = ZERO
             pres_term = ZERO
             do comp = 1, nspec
                xi_term = xi_term - &
                     eos_state%dhdX(comp)*rho_omegadot(i,j,k,comp)/eos_state%rho

                pres_term = pres_term + &
                     eos_state%dpdX(comp)*rho_omegadot(i,j,k,comp)/eos_state%rho
             enddo

             S_cc(i,j,k) = (sigma/eos_state%rho) * &
                  ( rho_Hext(i,j,k) + rho_Hnuc(i,j,k) + thermal(i,j,k) ) &
                  + sigma*xi_term &
                  + pres_term/(eos_state%rho*eos_state%dpdr)

             if (use_delta_gamma1_term) then
                delta_gamma1(i,j,k) = eos_state%gam1 - gamma1bar_cart(i,j,k)

                Ut_dot_er = &
                     u(i,j,k,1)*normal(i,j,k,1) + &
                     u(i,j,k,2)*normal(i,j,k,2) + &
                     u(i,j,k,3)*normal(i,j,k,3)

                dg1_term(i,j,k) = delta_gamma1(i,j,k)*Ut_dot_er* &
                     gradp0_cart(i,j,k)/ &
                          (gamma1bar_cart(i,j,k)**2*p0_cart(i,j,k))

             else
                dg1_term(i,j,k) = ZERO
                delta_gamma1(i,j,k) = ZERO
             end if

          enddo
       enddo
    enddo
    !$OMP END PARALLEL DO

  end subroutine make_S_cc_3d_sphr

  subroutine correct_delta_gamma1_term_1d(lo,hi,delta_gamma1_term,ng_dt, &
                                          delta_gamma1,ng_dg, &
                                          gamma1bar,psi,p0)

    integer         , intent(in   ) :: lo(:), hi(:), ng_dt, ng_dg
    real (kind=dp_t), intent(inout) :: delta_gamma1_term(lo(1)-ng_dt:)
    real (kind=dp_t), intent(in   ) ::      delta_gamma1(lo(1)-ng_dg:)
    real (kind=dp_t), intent(in   ) :: gamma1bar(0:)
    real (kind=dp_t), intent(in   ) :: psi(0:)
    real (kind=dp_t), intent(in   ) :: p0(0:)

    integer :: i

    do i = lo(1), hi(1)

       delta_gamma1_term(i) = delta_gamma1_term(i) &
            + delta_gamma1(i)*psi(i)/(gamma1bar(i)**2*p0(i))

    end do

  end subroutine correct_delta_gamma1_term_1d

  subroutine correct_delta_gamma1_term_2d(lo,hi,delta_gamma1_term,ng_dt, &
                                          delta_gamma1,ng_dg, &
                                          gamma1bar,psi,p0)

    integer         , intent(in   ) :: lo(:), hi(:), ng_dt, ng_dg
    real (kind=dp_t), intent(inout) :: delta_gamma1_term(lo(1)-ng_dt:,lo(2)-ng_dt:)
    real (kind=dp_t), intent(in   ) ::      delta_gamma1(lo(1)-ng_dg:,lo(2)-ng_dg:)
    real (kind=dp_t), intent(in   ) :: gamma1bar(0:)
    real (kind=dp_t), intent(in   ) :: psi(0:)
    real (kind=dp_t), intent(in   ) :: p0(0:)

    integer :: i, j

    do j = lo(2), hi(2)
       do i = lo(1), hi(1)

          delta_gamma1_term(i,j) = delta_gamma1_term(i,j) &
               + delta_gamma1(i,j)*psi(j)/(gamma1bar(j)**2*p0(j))

       end do
    end do

  end subroutine correct_delta_gamma1_term_2d

  subroutine correct_delta_gamma1_term_3d(lo,hi,delta_gamma1_term,ng_dt, &
                                          delta_gamma1,ng_dg, &
                                          gamma1bar,psi,p0)

    integer         , intent(in   ) :: lo(:), hi(:), ng_dt, ng_dg
    real (kind=dp_t), intent(inout) :: delta_gamma1_term(lo(1)-ng_dt:,lo(2)-ng_dt:,lo(3)-ng_dt:)
    real (kind=dp_t), intent(in   ) ::    delta_gamma1(lo(1)-ng_dg:,lo(2)-ng_dg:,lo(3)-ng_dg:)
    real (kind=dp_t), intent(in   ) :: gamma1bar(0:)
    real (kind=dp_t), intent(in   ) :: psi(0:)
    real (kind=dp_t), intent(in   ) :: p0(0:)

    integer :: i, j, k

    !$OMP PARALLEL DO PRIVATE(i,j,k)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             delta_gamma1_term(i,j,k) = delta_gamma1_term(i,j,k) &
                  + delta_gamma1(i,j,k)*psi(k)/(gamma1bar(k)**2*p0(k))

          end do
       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine correct_delta_gamma1_term_3d

  subroutine correct_delta_gamma1_term_3d_sphr(lo,hi,delta_gamma1_term,ng_dt, &
                                               delta_gamma1,ng_dg, &
                                               gamma1bar_cart,ng_g1, &
                                               psi_cart,ng_ps, &
                                               p0_cart,ng_p0)

    integer         , intent(in   ) :: lo(:), hi(:), ng_dt, ng_dg, ng_g1, ng_ps, ng_p0
    real (kind=dp_t), intent(inout) :: delta_gamma1_term(lo(1)-ng_dt:,lo(2)-ng_dt:,lo(3)-ng_dt:)
    real (kind=dp_t), intent(in   ) ::    delta_gamma1(lo(1)-ng_dg:,lo(2)-ng_dg:,lo(3)-ng_dg:)
    real (kind=dp_t), intent(in   ) ::  gamma1bar_cart(lo(1)-ng_g1:,lo(2)-ng_g1:,lo(3)-ng_g1:)
    real (kind=dp_t), intent(in   ) ::         p0_cart(lo(1)-ng_p0:,lo(2)-ng_p0:,lo(3)-ng_p0:)
    real (kind=dp_t), intent(in   ) ::        psi_cart(lo(1)-ng_ps:,lo(2)-ng_ps:,lo(3)-ng_ps:)

    integer :: i, j, k

    !$OMP PARALLEL DO PRIVATE(i,j,k)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             delta_gamma1_term(i,j,k) = delta_gamma1_term(i,j,k) &
                  + delta_gamma1(i,j,k)*psi_cart(i,j,k)/(gamma1bar_cart(i,j,k)**2*p0_cart(i,j,k))

          end do
       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine correct_delta_gamma1_term_3d_sphr

end module make_S_cc_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO