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

Annotation of /trunk/Sources/dyn3d/PPM3d/xtp.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: 2516 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 xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq, q, uc, fx1, xmass, &
2     iord)
3     DIMENSION uc(imr, *), dc(-iml:imr+iml+1), xmass(imr, jnp), fx1(imr+1), &
4     dq(imr, jnp), qtmp(-iml:imr+1+iml)
5     DIMENSION pu(imr, jnp), q(imr, jnp), isave(imr)
6    
7     imp = imr + 1
8    
9     ! van Leer at high latitudes
10     jvan = max(1, jnp/18)
11     j1vl = j1 + jvan
12     j2vl = j2 - jvan
13    
14     DO j = j1, j2
15    
16     DO i = 1, imr
17     qtmp(i) = q(i, j)
18     END DO
19    
20     IF (j>=jn .OR. j<=js) GO TO 2222
21     ! ************* Eulerian **********
22    
23     qtmp(0) = q(imr, j)
24     qtmp(-1) = q(imr-1, j)
25     qtmp(imp) = q(1, j)
26     qtmp(imp+1) = q(2, j)
27    
28     IF (iord==1 .OR. j==j1 .OR. j==j2) THEN
29     DO i = 1, imr
30     iu = float(i) - uc(i, j)
31     fx1(i) = qtmp(iu)
32     END DO
33     ELSE
34     CALL xmist(imr, iml, qtmp, dc)
35     dc(0) = dc(imr)
36    
37     IF (iord==2 .OR. j<=j1vl .OR. j>=j2vl) THEN
38     DO i = 1, imr
39     iu = float(i) - uc(i, j)
40     fx1(i) = qtmp(iu) + dc(iu)*(sign(1.,uc(i,j))-uc(i,j))
41     END DO
42     ELSE
43     CALL fxppm(imr, iml, uc(1,j), qtmp, dc, fx1, iord)
44     END IF
45    
46     END IF
47    
48     DO i = 1, imr
49     fx1(i) = fx1(i)*xmass(i, j)
50     END DO
51    
52     GO TO 1309
53    
54     ! ***** Conservative (flux-form) Semi-Lagrangian transport *****
55    
56     2222 CONTINUE
57    
58     DO i = -iml, 0
59     qtmp(i) = q(imr+i, j)
60     qtmp(imp-i) = q(1-i, j)
61     END DO
62    
63     IF (iord==1 .OR. j==j1 .OR. j==j2) THEN
64     DO i = 1, imr
65     itmp = int(uc(i,j))
66     isave(i) = i - itmp
67     iu = i - uc(i, j)
68     fx1(i) = (uc(i,j)-itmp)*qtmp(iu)
69     END DO
70     ELSE
71     CALL xmist(imr, iml, qtmp, dc)
72    
73     DO i = -iml, 0
74     dc(i) = dc(imr+i)
75     dc(imp-i) = dc(1-i)
76     END DO
77    
78     DO i = 1, imr
79     itmp = int(uc(i,j))
80     rut = uc(i, j) - itmp
81     isave(i) = i - itmp
82     iu = i - uc(i, j)
83     fx1(i) = rut*(qtmp(iu)+dc(iu)*(sign(1.,rut)-rut))
84     END DO
85     END IF
86    
87     DO i = 1, imr
88     IF (uc(i,j)>1.) THEN
89     ! DIR$ NOVECTOR
90     DO ist = isave(i), i - 1
91     fx1(i) = fx1(i) + qtmp(ist)
92     END DO
93     ELSE IF (uc(i,j)<-1.) THEN
94     DO ist = i, isave(i) - 1
95     fx1(i) = fx1(i) - qtmp(ist)
96     END DO
97     ! DIR$ VECTOR
98     END IF
99     END DO
100     DO i = 1, imr
101     fx1(i) = pu(i, j)*fx1(i)
102     END DO
103    
104     ! ***************************************
105    
106     1309 fx1(imp) = fx1(1)
107     DO i = 1, imr
108     dq(i, j) = dq(i, j) + fx1(i) - fx1(i+1)
109     END DO
110    
111     ! ***************************************
112    
113     END DO
114     RETURN
115     END SUBROUTINE xtp
116    
  ViewVC Help
Powered by ViewVC 1.1.21