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

Contents of /trunk/Sources/dyn3d/PPM3d/lmtppm.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: 1876 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 lmtppm(dc, a6, ar, al, p, im, lmt)
2
3 ! A6 = CURVATURE OF THE TEST PARABOLA
4 ! AR = RIGHT EDGE VALUE OF THE TEST PARABOLA
5 ! AL = LEFT EDGE VALUE OF THE TEST PARABOLA
6 ! DC = 0.5 * MISMATCH
7 ! P = CELL-AVERAGED VALUE
8 ! IM = VECTOR LENGTH
9
10 ! OPTIONS:
11
12 ! LMT = 0: FULL MONOTONICITY
13 ! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
14 ! LMT = 2: POSITIVE-DEFINITE CONSTRAINT
15
16 PARAMETER (r12=1./12.)
17 DIMENSION a6(im), ar(im), al(im), p(im), dc(im)
18
19 IF (lmt==0) THEN
20 ! Full constraint
21 DO i = 1, im
22 IF (dc(i)==0.) THEN
23 ar(i) = p(i)
24 al(i) = p(i)
25 a6(i) = 0.
26 ELSE
27 da1 = ar(i) - al(i)
28 da2 = da1**2
29 a6da = a6(i)*da1
30 IF (a6da<-da2) THEN
31 a6(i) = 3.*(al(i)-p(i))
32 ar(i) = al(i) - a6(i)
33 ELSE IF (a6da>da2) THEN
34 a6(i) = 3.*(ar(i)-p(i))
35 al(i) = ar(i) - a6(i)
36 END IF
37 END IF
38 END DO
39 ELSE IF (lmt==1) THEN
40 ! Semi-monotonic constraint
41 DO i = 1, im
42 IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 150
43 IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN
44 ar(i) = p(i)
45 al(i) = p(i)
46 a6(i) = 0.
47 ELSE IF (ar(i)>al(i)) THEN
48 a6(i) = 3.*(al(i)-p(i))
49 ar(i) = al(i) - a6(i)
50 ELSE
51 a6(i) = 3.*(ar(i)-p(i))
52 al(i) = ar(i) - a6(i)
53 END IF
54 150 END DO
55 ELSE IF (lmt==2) THEN
56 DO i = 1, im
57 IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 250
58 fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12
59 IF (fmin>=0.) GO TO 250
60 IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN
61 ar(i) = p(i)
62 al(i) = p(i)
63 a6(i) = 0.
64 ELSE IF (ar(i)>al(i)) THEN
65 a6(i) = 3.*(al(i)-p(i))
66 ar(i) = al(i) - a6(i)
67 ELSE
68 a6(i) = 3.*(ar(i)-p(i))
69 al(i) = ar(i) - a6(i)
70 END IF
71 250 END DO
72 END IF
73 RETURN
74 END SUBROUTINE lmtppm
75
  ViewVC Help
Powered by ViewVC 1.1.21