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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 166 - (hide 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 guez 166 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