zeus_python_old/rhs.src at master · Higginbottom/zeus_python_old · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include "zeus2d.def"
c=======================================================================
c/////////////////////////  SUBROUTINE RHS  \\\\\\\\\\\\\\\\\\\\\\\\\\\\
c
      subroutine rhs(fn1,fn2)
#ifdef RAD
c
c  PURPOSE: Computes the finite differenced radiation and gas energy
c  equations as fn1(e,en,er,ern) and fn2(e,en,er,ern). These functions
c  are used as the RHS of the Newton-Raphson matrix eqn.  Note the -ve
c  sign appearing in the definition of the Jacobian is added when
c  the moment eqn matrix is set up (RICCG), ie Ax = -rhs.
c  The N-R scheme will continue to iterate en and ern until fn1=fn2=0.
c
c  INPUT ARGUMENTS:[none]
c
c  OUTPUT ARGUMENTS:
c   fn1 = radiation energy equation
c   fn2 = gas       energy equation
c
c  EXTERNALS:[none]
c
c  LOCALS:
c-----------------------------------------------------------------------
      implicit NONE
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "scratch.h"
#include "radiation.h"
      REAL fn1(in,jn),fn2(in,jn)
c
      integer i,j
      REAL erb(in,jn)
      REAL qa,qb,qc,qd,qe
c=======================================================================
c  INLINE FUNCTIONS
c   TH BAR: time centers a variable according to the value of radth
c=======================================================================
      REAL thbar,xn,xnp1
      th bar(xn,xnp1) = radth *xnp1 + (1.0-radth )*xn
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\////////////////////////////////
c=======================================================================
c  Compute erb = thbar(er)
c
      do 20 j=js-1,je+1
	do 10 i=iim1(j),iop1(j)
	  erb(i,j) = thbar(er(i,j),ern(i,j))
10      continue
20    continue
c  Compute the GRAD(v):P term
c
      do 40 j=js,je
	do 30 i=ii(j),io(j)
          fn1(i,j) = (dv11(i,j)*f11(i,j) + dv22(i,j)*f22(i,j)
     &    + dv33(i,j)*(1.0-f11(i,j)-f22(i,j))
     &    + (dv12(i,j) + dv21(i,j))*f12(i,j))*erb(i,j)
30      continue
40    continue
c
c  Now sum DIV(F) terms into moment eqns.
c
      do 60 j=js,je
	do 50 i=ii(j),io(j)
          qa = g2a(i+1)*g31a(i+1)*dx1bi(i+1)*dr1(i+1,j)
          qb = g2a(i  )*g31a(i  )*dx1bi(i  )*dr1(i  ,j)
          qc = g2b(i)*g2b(i)
          qd = g32a(j+1)*dx2bi(j+1)*dr2(i,j+1)
          qe = g32a(j  )*dx2bi(j  )*dr2(i,j  )
	  fn1(i,j) = fn1(i,j)
     &  -(qa*(erb(i+1,j)-erb(i,j))-qb*(erb(i,j)-erb(i-1,j)))/dvl1a(i)
     &  -(qd*(erb(i,j+1)-erb(i,j))-qe*(erb(i,j)-erb(i,j-1)))/qc/dvl2a(j)
50      continue
60    continue
c
c  Now sum remaining terms into moment eqns.
c
      do 80 j=js,je
        do 70 i=ii(j),io(j)
	  qa = 4.0*pi*thbar(kap(i,j),kapn(i,j))*thbar(bb(i,j),bbn(i,j))
     &            - c*thbar(kap(i,j),kapn(i,j))*erb(i,j)
	  fn1(i,j) = ern(i,j) - er(i,j) - dt*(qa - fn1(i,j))
	  fn2(i,j) = en(i,j) - e(i,j)
     &             + dt*(qa + thbar(p(i,j),pn(i,j))*divv(i,j))
70      continue
80    continue
#endif
      return
      end