/[lmdze]/trunk/Sources/dyn3d/PPM3d/xadv.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/PPM3d/xadv.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 166 - (show annotations)
Wed Jul 29 14:32:55 2015 UTC (8 years, 10 months ago) by guez
File size: 1922 byte(s) 
Split ppm3d.f into files containing a single procedure.

Factorized computations of filtering matrices into a procedure
inifilr_hemisph. Had then to change the matrices from allocatable to
pointer and from customized lower bound to lower bound 1. The change
in lower bounds does not matter because the matrices are only used as
a whole as actual arguments.

Also, in infilr, instead of finding jfilt[ns][uv] from approximately
jjm /2, start at index j1 that corresponds to the equator. This is not
the same if there is a zoom in latitude.

Also, the test (rlamda(modfrst[ns][uv](j)) * cos(rlat[uv](j)) < 1) in
the loops on filtered latitudes is not useful now that we start from
j1: it is necessarily true. See notes.

Just encapsulated lwvn into a module and removed unused argument ktraer.
1 SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad)
2 REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp)
3
4 jmr = jnp - 1
5 DO j = j1, j2
6 IF (j>js .AND. j<jn) GO TO 1309
7
8 DO i = 1, imr
9 qtmp(i) = p(i, j)
10 END DO
11
12 DO i = -iml, 0
13 qtmp(i) = p(imr+i, j)
14 qtmp(imr+1-i) = p(1-i, j)
15 END DO
16
17 IF (iad==2) THEN
18 DO i = 1, imr
19 ip = nint(ua(i,j))
20 ru = ip - ua(i, j)
21 ip = i - ip
22 a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
23 b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
24 adx(i, j) = qtmp(ip) + ru*(a1*ru+b1)
25 END DO
26 ELSE IF (iad==1) THEN
27 DO i = 1, imr
28 iu = ua(i, j)
29 ru = ua(i, j) - iu
30 iiu = i - iu
31 IF (ua(i,j)>=0.) THEN
32 adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
33 ELSE
34 adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
35 END IF
36 END DO
37 END IF
38
39 DO i = 1, imr
40 adx(i, j) = adx(i, j) - p(i, j)
41 END DO
42 1309 END DO
43
44 ! Eulerian upwind
45
46 DO j = js + 1, jn - 1
47
48 DO i = 1, imr
49 qtmp(i) = p(i, j)
50 END DO
51
52 qtmp(0) = p(imr, j)
53 qtmp(imr+1) = p(1, j)
54
55 IF (iad==2) THEN
56 qtmp(-1) = p(imr-1, j)
57 qtmp(imr+2) = p(2, j)
58 DO i = 1, imr
59 ip = nint(ua(i,j))
60 ru = ip - ua(i, j)
61 ip = i - ip
62 a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
63 b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
64 adx(i, j) = qtmp(ip) - p(i, j) + ru*(a1*ru+b1)
65 END DO
66 ELSE IF (iad==1) THEN
67 ! 1st order
68 DO i = 1, imr
69 ip = i - ua(i, j)
70 adx(i, j) = ua(i, j)*(qtmp(ip)-qtmp(ip+1))
71 END DO
72 END IF
73 END DO
74
75 IF (j1/=2) THEN
76 DO i = 1, imr
77 adx(i, 2) = 0.
78 adx(i, jmr) = 0.
79 END DO
80 END IF
81 ! set cross term due to x-adv at the poles to zero.
82 DO i = 1, imr
83 adx(i, 1) = 0.
84 adx(i, jnp) = 0.
85 END DO
86 RETURN
87 END SUBROUTINE xadv
88
  ViewVC Help
Powered by ViewVC 1.1.21