zeus_python_old/diskwbc.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
#include "zeus2d.def"
c=====================================================================
c///////////////////  SUBROUTINE DISKWBC \\\\\\\\\\\\\\\\\\\\\\\\\\\
c
      subroutine diskwbc
c
c  SETS BOUNDARY CONDITIONS AT BASE OF a DISK WIND
c
c     written by: Daniel Proga
c     date:       May 2008
c     modified1:
c
c  INPUT PARAMETERS:
c
c----------------------------------------------------------------------
      implicit NONE
#include "cons.h"
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "scratch.h"
#include "gravity.h"
#include "bndry.h"
#include "diskw.h"


      integer i,j
	  real*8 tfloor
c
      external bvalv1,bvalv2,bvalv3,bvald,bvale
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////////////
c=======================================================================
c
c  set d, e, velocity at last cell in j
c
       dfloor=1d-6
c	print *,dfloor*d0
       do  i=is,ie
        do  j=is,je
        if (d (i,j).le.1e-22) d(i,j)=1e-22
		tfloor=((3.0*ptmass*mdot*guniv)/  (8.0*pi*sb*
     &   (x1b(i)*sin(x2b(j)))**3.0) )**0.25
		 if (e(i,j) .le. (3.0*d(i,j)*boltz*tfloor)/(2.0*mp*mu)) then
			 e(i,j)=(3.0*d(i,j)*boltz*tfloor)/(2.0*mp*mu)
		 endif
        enddo
       enddo

       do  i=is,ie
        d(i,je)=djop1(i)
        v1(i,je)=0.0
        v3(i,je)=v3jop1(i)
c         if (e(i,je).le.ejop1(i)) e(i,je)=ejop1(i)
       enddo

       call bvald
       call bvale
       call bvalv1
       call bvalv2

      return
      end