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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 166 - (show annotations)
Wed Jul 29 14:32:55 2015 UTC (8 years, 9 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 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