heatcool.src

#include "zeus2d.def"


	subroutine hc_sub(rho,xi0,e_0,e_1,i,j)

       implicit NONE

#include "cons.h"
#include "param.h"
#include "grid.h"
#include "field.h"
#include "gravity.h"

#include "root.h"
#include "diskw.h"
        REAL    xi0,temp,hc_tot,rho,rho1,xi1,hc_tot1,e_0,e_1
	REAL    t_0,t_1,hc_01
	REAL 	e_01,t_01,sqsqxi1
	INTEGER i1,j1,i,j
	double precision t_l,t_u,tol,nh1,ne1,tfloor
	REAL heatcool,zbrent,zfunc,n1
	real test
	common /hc_data/ xi1,e_01,hc_01,sqsqxi1,n1,i1,j1,nh1,ne1
	external zfunc
	
c   first calculate the current temperature of the cell	
	
	
	t_01=e_0*(2.0/3.0)/((rho/(mu*mp))*boltz)
	xi1=xi0
	i1=i
	j1=j
	e_01=e_0
	
c   compute xi to the quater power, this is used a lot and doing this now saves time	
	sqsqxi1=sqrt(sqrt(xi0))

c   n1 is the particle density (including electrons)
	n1=rho/(mu*mp)
c   nh1 is the hydrogen density (based upon solar abundances)	
	nh1=rho/(1.43*mp)
c   ne1 is the electron density 	
	ne1=nh1*1.21


c stage one, bracket to find out where the zero point is


c we use the heating rate at the   originaltemperature a lot,
c lets compute it first
	hc_01=heatcool(t_01)


c   the next few lines bracket the solution temperature
	t_l=t_01*0.9
	t_u=t_01*1.1	
c	print *,"Bracketing temp=",t_01,i,j,rho,e_0
	test=zfunc(t_l)*zfunc(t_u)
	tol=1.0
c	if (test.gt.0) then
c		print *,"PANIC"
c	endif
c	print *,"New run"
	  do while (test > 0)
		t_l=t_l*0.9
		t_u=t_u*1.1
    		test=zfunc(t_l)*zfunc(t_u)
  	end do
c now we have worked out a bracket, find the correct value of new temperature
c we search for the zero point of the equation e_old - e_old - dt(de(t(1/2)))	
c	print *,"XXXXXX calling zbrent",t_l,t_u
	     t_1=zbrent(zfunc,t_l,t_u,tol)
c we will call heatcool one last time to compute the heating and cooling rates for the new temperature	
		 
       tfloor=((3.0*ptmass*mdot*guniv)/  (8.0*pi*sb*
     &   (x1b(i)*sin(x2b(j)))**3.0) )**0.25
	  
	  if (t_1.lt.tfloor) then
c	 print *,"temperature returned from cell",i,j," too low=",t_1
	  t_1=tfloor
      endif
	
	hc_01=heatcool(t_1)




c now we have a new temperature, we will return the new internal energy. 
	e_1=t_1*1.5*n1*boltz

c	call exit()

	end



      double precision function heatcool(temp)
      implicit NONE
#include "cons.h"
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "diskw.h"
      integer i,j
      REAL    div1,div2,divv,del,pnew,pav,xi0,temp,hc_tot,rho,e_0,t_0
      REAL    gamma_comp,gamma_x,l_b,l_l,hc_0,sqsqxi,st,n,nh,ne
	common /hc_data/ xi0,e_0,hc_0,sqsqxi,n,i,j,nh,ne

	st=sqrt(temp)




	h_comp(i,j)=comp_h_pre(i,j)*8.9e-36*xi0*tx
	c_comp(i,j)=comp_c_pre(i,j)*8.9e-36*xi0*(4.0*temp)
	gamma_comp=h_comp(i,j)-c_comp(i,j)
c	if (i.eq.3 .and. j.eq.90) then
c		print *,ne,nh,xi0,temp,c_comp(i,j),comp_c_pre(i,j)
c	endif

c	gamma_x=sqrt(xi)
c	gamma_x=sqrt(gamma_x)
c	gamma_x=gamma_x*1.5e-21
c	gamma_x=gamma_x/sqrt(temp)
c	gamma_x=gamma_x*(1-temp/tx)
	h_xray(i,j)=xray_pre(i,j)*1.5e-21*(sqsqxi/st)*(1-(temp/tx))
	gamma_x=h_xray(i,j)

	l_b=brem_pre(i,j)*3.3e-27*st
	c_brem(i,j)=l_b
c	l_l=exp(-1.3e5/temp)
c	l_l=l_l/xi
c	l_l=l_l/sqrt(temp)
c	l_l=l_l*1.7e-18
c	l_l=l_l+1e-24
	l_l=line_pre(i,j)*(1.7e-18*exp(-1.3e5/temp)/(xi0*st)+1e-24)
	c_line(i,j)=l_l
c	print *,temp,xi,gamma_comp,gamma_x,l_b,l_l
	hc_tot=ne*nh*(gamma_comp+gamma_x-l_b-l_l)
c	print *,ne,nh,gamma_comp,gamma_x,l_b,l_l
c	print *,ne*nh*(gamma_comp+gamma_x-l_b-l_l),hc_tot
	
	
	heatcool=hc_tot
	
c	print *,xi0/1e15,temp/1e8,heatcool,(gamma_comp+gamma_x-l_b-l_l),
c    + ne,nh
	
	end function



	double precision function zfunc(temp) 
	implicit NONE
#include "cons.h"
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "diskw.h"
	REAL rho,xi0,e_0,t_0,temp,heatcool,hc_0,logxi,sqsqxi,n,nh,ne
	INTEGER i,j
	common /hc_data/ xi0,e_0,hc_0,sqsqxi,n,i,j,nh,ne
	external heatcool
	zfunc=(1.5*temp*n*boltz)-e_0-
     &  dt*(hc_0+heatcool(temp))/2.0
c	print *,temp,zfunc
	end function


c     this file contains the quasi-vanilla and slightly modified
c     versions of FORTRAN Numerical Recipe routines

c     DOUBLE PRECISION FUNCTION ZBRENT(FUNC,X1,X2,TOL)

c     MODIFICATIONS from the vanilla Numerical Recipes: all of the
c     following routines have been trivially modified in two ways from
c     the original vanilla; (1) the declarations are now IMPLICIT NONE
c     or IMPLICIT DOUBLE PRECISION; (2) communication of errors which
c     may occur in the form of WRITE statements with hard STOPs; some of
c     the routines have been further modified (see comments in each
c     routine for details)

c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

      REAL FUNCTION ZBRENT(FUNC,X1,X2,TOL)

c     return the root of the function FUNC using Brents method where
c     the root of the function is assumed lie between X1 and X2; the
c     root is refined until its accuracy is TOL and is returned as
c     ZBRENT

c..............................................................................

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ITMAX=100,EPS=3.E-14)
      A=X1
      B=X2
      FA=FUNC(A)
      FB=FUNC(B)
c     error: Root not bracketed 
      IF(FB*FA.GT.0.) THEN
       WRITE(6,*) 'ERROR(zbrent): root of FUNC not bracketed.'
       WRITE(6,*) 'lower bracket A = ',A,' FUNC(A) = ',FUNC(A)
       WRITE(6,*) 'upper bracket B = ',B,' FUNC(B) = ',FUNC(B)
       WRITE(6,*) 'this can be remedied by changing A and/or B'
       call exit(5)
      END IF
      FC=FB
      DO 11 ITER=1,ITMAX
        IF(FB*FC.GT.0.) THEN
          C=A
          FC=FA
          D=B-A
          E=D
        ENDIF
        IF(ABS(FC).LT.ABS(FB)) THEN
          A=B
          B=C
          C=A
          FA=FB
          FB=FC
          FC=FA
        ENDIF
        TOL1=2.*EPS*ABS(B)+0.5*TOL
        XM=.5*(C-B)
        IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
          ZBRENT=B
c	print *,"zbrent took ",ITER," iterations"
          RETURN
        ENDIF
        IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
          S=FB/FA
          IF(A.EQ.C) THEN
            P=2.*XM*S
            Q=1.-S
          ELSE
            Q=FA/FC
            R=FB/FC
            P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
            Q=(Q-1.)*(R-1.)*(S-1.)
          ENDIF
          IF(P.GT.0.) Q=-Q
          P=ABS(P)
          IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
            E=D
            D=P/Q
          ELSE
            D=XM
            E=D
          ENDIF
        ELSE
          D=XM
          E=D
        ENDIF
        A=B
        FA=FB
        IF(ABS(D) .GT. TOL1) THEN
          B=B+D
        ELSE
          B=B+SIGN(TOL1,XM)
        ENDIF
        FB=FUNC(B)
11    CONTINUE
c     error message for too many iterations
      WRITE(6,*) 'ERROR(zbrent): number of iterations required to'
      WRITE(6,*) 'converge on the root of FUNC exceeded maximum'
      WRITE(6,*) 'allowed.  This can be remedied by increasing'
      WRITE(6,*) 'parameter ITMAX in the file or by relaxing'
      WRITE(6,*) 'the convrgence tolerance TOL'
      WRITE(6,*) '*** terminating program  ***'
      call exit(5)
      ZBRENT=B
      RETURN
      END
	  
	  
	  subroutine hc_file_input()
#include "cons.h"
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "diskw.h"
#include "control.h"
		  integer i,j
		  
		  
		  LOGICAL THERE 
   
		  real r,theta,c_heat,x_heat,b_heat,l_heat,c_cool,x_cool,py_xi
          real b_cool,l_cool,n1,py_temp,py_ne,py_vol,ne,nh
          real comp_c_pre_old,comp_h_pre_old,xray_pre_old,line_pre_old,
     +     brem_pre_old,ratio,py_rho,py_nh
	 integer pre_limit,pre_limit_tmp
	 
	 pre_limit_tmp=0
	 pre_limit=0
		  
c		  print *,"Got here OK"
         INQUIRE( FILE='py_heatcool.dat', EXIST=THERE ) 
          IF ( THERE ) then
		 	 open(66,file='py_heatcool.dat',status='old')
	  	   	read (66,*,end=20)
   		  	 do 10, n=1,10000000
 
		  	   	read (66,*,end=20) i,j,r,theta,py_vol,py_temp,py_xi,py_ne,
     &			x_heat,c_heat,l_heat,
     &   b_heat,c_cool,l_cool,b_cool,py_rho,py_nh

c				print *,i,j,r,theta,x_heat,c_heat,l_heat,b_heat,x_cool,
c     &   c_cool,l_cool,b_cool,py_xi	   
				
				n1=d(i+2,j+3)/(mu*mp)
				
				nh=d(i+2,j+3)/(1.43*mp)
				ne=nh*1.21
				
c		  print *,hc_xray(i+2,j+3),x_heat/(n1*n1*1.2),x_cool/(n1*n1*1.2),
c     &   x_heat/(n1*n1*1.2)-x_cool/(n1*n1*1.2)
c		  do 10 i=is,ie
c			  do 20 j=js,je
c		  print *,is,i+2,ie,x1b(i+2)/r,js,j+3,je,x2b(j+3)/
c     &  (theta)
				comp_c_pre_old=comp_c_pre(i+2,j+3) 
                comp_h_pre_old=comp_h_pre(i+2,j+3)
                xray_pre_old=xray_pre(i+2,j+3)
                line_pre_old=line_pre(i+2,j+3)
                brem_pre_old=brem_pre(i+2,j+3)      
                        
                        
                        
                        
                         
		 comp_h_pre(i+2,j+3)=comp_h_pre(i+2,j+3)*(c_heat)/
     &   (ne*nh*h_comp(i+2,j+3))+1e-3
		 comp_c_pre(i+2,j+3)=comp_c_pre(i+2,j+3)*(c_cool)/
     &   (ne*nh*c_comp(i+2,j+3))+1e-3
		 
c	 if (i.eq.1) then
c	  print *,"COMP test",i,j,comp_c_pre(i+2,j+3),
c     &   c_comp(i+2,j+3),ne*nh*c_comp(i+2,j+3),c_cool,ne,nh,
c     &    xi(i+2,j+3),d(i+2,j+3),
c     &   e(i+2,j+3)*(2.0/3.0)/((d(i+2,j+3)/(mu*mp))*boltz)
c        endif
		 xray_pre(i+2,j+3)=xray_pre(i+2,j+3)*(x_heat)/
     &   (ne*nh*h_xray(i+2,j+3))+1e-3
 		 line_pre(i+2,j+3)=line_pre(i+2,j+3)*l_cool/
     &   (ne*nh*(c_line(i+2,j+3)+c_xray(i+2,j+3)))+1e-3
 		 brem_pre(i+2,j+3)=brem_pre(i+2,j+3)*b_cool/
     &   (ne*nh*c_brem(i+2,j+3))+1e-3
        if (hc_lim.gt.0.0) then
        ratio=hc_lim
	  
	  
c hc_lim is negative, which means we dont want to limit hc changes	  
c	  print *,"CHECKING ",i,j
      if (comp_h_pre(i+2,j+3)/comp_h_pre_old .lt. ratio) then
c		  print *,"comp_h_pre=",comp_h_pre(i+2,j+3),' LIMITING LO'
           comp_h_pre(i+2,j+3)=comp_h_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
      endif
      if (comp_c_pre(i+2,j+3)/comp_c_pre_old .lt. ratio) then
c		  print *,"comp_c_pre=",comp_c_pre(i+2,j+3),' LIMITING LO'
		  
           comp_c_pre(i+2,j+3)=comp_c_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
		   
      endif
      if (xray_pre(i+2,j+3)/xray_pre_old .lt. ratio) then
c		  print *,"xray_pre=",xray_pre(i+2,j+3),' LIMITING LO'
		  
           xray_pre(i+2,j+3)=xray_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
		   
      endif
      if (line_pre(i+2,j+3)/line_pre_old .lt. ratio) then
c		  print *,"line_pre=",line_pre(i+2,j+3),' LIMITING LO'
		  
            line_pre(i+2,j+3)=line_pre_old*ratio
 		   pre_limit_tmp=pre_limit_tmp+1
			
      endif
      if (brem_pre(i+2,j+3)/brem_pre_old .lt. ratio) then
c		  print *,"brem_pre=",brem_pre(i+2,j+3),' LIMITING LO'
		  
            brem_pre(i+2,j+3)=brem_pre_old*ratio
 		   pre_limit_tmp=pre_limit_tmp+1
			
      endif
      
      
      ratio=1.0/hc_lim
     
      if (comp_h_pre(i+2,j+3)/comp_h_pre_old .gt. ratio) then
c		  print *,"comp_h_pre=",comp_h_pre(i+2,j+3),' LIMITING HI'
		  
           comp_h_pre(i+2,j+3)=comp_h_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
		   
      endif
      if (comp_c_pre(i+2,j+3)/comp_c_pre_old .gt. ratio) then
c		  print *,"comp_h_pre=",comp_c_pre(i+2,j+3),' LIMITING HI'
		  
           comp_c_pre(i+2,j+3)=comp_c_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
		   
      endif
      if (xray_pre(i+2,j+3)/xray_pre_old .gt. ratio) then
		   xray_pre(i+2,j+3)=xray_pre_old*ratio
		   pre_limit_tmp=pre_limit_tmp+1
		   
      endif
      if (line_pre(i+2,j+3)/line_pre_old .gt. ratio) then
c		  print *,"line_pre=",line_pre(i+2,j+3),' LIMITING HI'
		  
            line_pre(i+2,j+3)=line_pre_old*ratio
 		   pre_limit_tmp=pre_limit_tmp+1
			
      endif
      if (brem_pre(i+2,j+3)/brem_pre_old .gt. ratio) then
c		  print *,"brem_pre=",brem_pre(i+2,j+3),' LIMITING HI'
		  
            brem_pre(i+2,j+3)=brem_pre_old*ratio		
			   pre_limit_tmp=pre_limit_tmp+1

      endif
      else
		  print *,"We are not limiting changes in HC rates"
	  endif
     
c	  print *,'COMP_C_PRE    old=',comp_c_pre_old,
c     c' new=',comp_c_pre(i+2,j+3) 
c      print *,'COMP_H_PRE    old=',comp_h_pre_old,
c     c' new=',comp_h_pre(i+2,j+3)
c      print *,'XRAY_PRE    old=',xray_pre_old,
c     c' new=',xray_pre(i+2,j+3)
c      print *,'LINE_PRE    old=',line_pre_old,
c     c' new=',line_pre(i+2,j+3)
c      print *,'BREM_PRE    old=',brem_pre_old,
c     c' new=',brem_pre(i+2,j+3)  
	 
	 
     
     
     
c      print *,i,j,ne,py_ne,comp_cool_pre_old,comp_c_pre(i+2,j+3),
c     +  ne*nh*c_comp(i+2,j+3),c_cool
c      print *,x_heat,c_heat,l_heat,
c     &   b_heat,c_cool,l_cool,b_cool
c      print *,i,j,py_temp,e(i+2,j+3)/(1.5*n1*boltz),
c     + comp_cool_pre_old,comp_c_pre(i+2,j+3)

c		  print *,i+3,j+3,r,theta
c 20 		continue
		if (pre_limit_tmp .gt. 0) then
			pre_limit=pre_limit+1
			pre_limit_tmp=0
		endif
10		continue
		
20		print *,"There are ",n," records"	
		print *,pre_limit," cells had HC prefactor change limited"
		
		  close unit=66
	  else
		  print *,"There is no heatcool file - continuuing normally"
	  endif
		  
		  
		  RETURN
	  END