zeus_python_old/newvg.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
91
92
93
94
95
96
97
98
99
100
101
102
#include "zeus2d.def"
c=======================================================================
c///////////////////////////  SUBROUTINE NEWVG  \\\\\\\\\\\\\\\\\\\\\\\\
c
      subroutine newvg
c
c  PURPOSE: Computes grid velocities at current timestep to be used in
c  updating the grid variables.  Velocoties are claculated over the
c  range i=is-2,ie+2 and j=js-2,je+2.  The method used to compute the
c  grid velocities depends on the flag igcon and the sign of x1fac,x2fac
c        igcon = 1 gives "lagrangian" tracking in x1 [x2] lines
c        igcon = 2 for input grid boundary speeds
c                vg1(io) = x1fac * central soundspeed
c                vg2(jo) = x2fac * central soundspeed
c        igcon = 3 for uniform translation
c
c  LOCALS:
c-----------------------------------------------------------------------
      implicit NONE
#include "param.h"
#include "grid.h"
#include "field.h"
#include "bndry.h"
#include "root.h"
c
      integer i,j
      REAL    qa,qb
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\//////////////////////////////////////
c=======================================================================
c
c  igcon = 0: "Lagrangian" tracking of grid lines
c
      if (igcon .eq. 1) then
      if (nx1z .gt. 1 .and. x1fac .lt. 0.0) then
        do 100 i=is+1,ie
          vg1(i) = v1(i,js)
100     continue
        vg1(is  ) = 0.0
        vg1(is-1) = -vg1(is+1)
        vg1(is-2) = -vg1(is+2)
        vg1(ie+1) = 0.0
        vg1(ie+2) = -vg1(ie  )
      endif
c
      if (nx2z .gt. 1 .and. x2fac .lt. 0.0) then
        do 110 j=js+1,je
          vg2(j) = v2(is,j)
110     continue
        vg2(js  ) = 0.0
        vg2(js-1) = -vg2(js+1)
        vg2(js-2) = -vg2(js+2)
        vg2(je+1) = 0.0
        vg2(je+2) = -vg2(je  )
      endif
      return
      endif
c
c  igcon=2:  qa=central sound speed; vg is computed as a
c  linear function of x from x=0 to x(outer boundary).
c
      if (igcon .eq. 2) then
      qa = sqrt(gamma*(gamma-1.0)*e(is,js)/d(is,js))
      if (x1fac .ne. 0.0) then
        qb = x1fac*qa/x1a(ie)
        do 200 i=is+1,ie
          vg1(i)=-qb*x1a(i)
200     continue
        vg1(is  ) = 0.0
        vg1(is-1) = -vg1(is+1)
        vg1(is-2) = -vg1(is+2)
        vg1(ie+1) = 0.0
        vg1(ie+2) = -vg1(ie  )
      end if
c
      if (x2fac .ne. 0.0) then
        qb = x2fac*qa/x2a(je)
        do 210 j=js+1,je
          vg2(j)=-qb*x2a(j)
210     continue
        vg2(js  ) = 0.0
        vg2(js-1) = -vg2(js+1)
        vg2(js-2) = -vg2(js+2)
        vg2(je+1) = 0.0
        vg2(je+2) = -vg2(je  )
      end if
      return
      endif
c
c  igcon = 3.  Uniform translation of the grid at a velocity x1[2]fac.
c
      if (igcon .eq. 3) then
      do 300 i = is-2, ie+2
        vg1(i) = x1fac
300   continue
      do 310 j = js-2, je+2
        vg2(j) = x2fac
310   continue
      return
      endif
c
      return
      end