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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 176 - (show annotations)
Tue Feb 23 17:00:39 2016 UTC (8 years, 2 months ago) by guez
File size: 2988 byte(s)
1e-60 was underflowing in simple precision.
1 SUBROUTINE qckxyz(q, qtmp, imr, jnp, nlay, j1, j2, cosp, acosp, cross, ic, &
2 nstep)
3
4 implicit none
5
6 real, PARAMETER:: my_tiny=tiny(0.)
7 integer imr, jnp, nlay, j1, j2, ic, nstep
8 real q(imr, jnp, nlay), qtmp(imr, jnp), cosp(*), acosp(*)
9 LOGICAL cross
10 real dup, qly, qup, sum
11 integer i, icr, ip, ipx, ipy, l, len, nlaym1
12
13 nlaym1 = nlay - 1
14 len = imr*(j2-j1+1)
15 ip = 0
16
17 ! Top layer
18 l = 1
19 icr = 1
20 CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, my_tiny)
21 IF (ipy/=0) then
22 CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, my_tiny)
23 IF (ipx/=0) then
24
25 IF (cross) THEN
26 CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, my_tiny)
27 END IF
28 IF (icr/=0) then
29
30 ! Vertical filling...
31 DO i = 1, len
32 IF (q(i,j1,1)<0.) THEN
33 ip = ip + 1
34 q(i, j1, 2) = q(i, j1, 2) + q(i, j1, 1)
35 q(i, j1, 1) = 0.
36 END IF
37 END DO
38 end IF
39 end IF
40 end IF
41
42 DO l = 2, nlaym1
43 icr = 1
44
45 CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, my_tiny)
46 IF (ipy/=0) then
47 CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, my_tiny)
48 IF (ipx/=0) then
49 IF (cross) THEN
50 CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, my_tiny)
51 END IF
52 IF (icr/=0) then
53
54 DO i = 1, len
55 IF (q(i,j1,l)<0.) THEN
56
57 ip = ip + 1
58 ! From above
59 qup = q(i, j1, l-1)
60 qly = -q(i, j1, l)
61 dup = min(qly, qup)
62 q(i, j1, l-1) = qup - dup
63 q(i, j1, l) = dup - qly
64 ! Below
65 q(i, j1, l+1) = q(i, j1, l+1) + q(i, j1, l)
66 q(i, j1, l) = 0.
67 END IF
68 END DO
69 end IF
70 end IF
71 end IF
72 END DO
73
74 ! BOTTOM LAYER
75 sum = 0.
76 l = nlay
77
78 CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, my_tiny)
79 IF (ipy/=0) then
80 CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, my_tiny)
81 IF (ipx/=0) then
82
83 CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, my_tiny)
84 IF (icr/=0) then
85
86 DO i = 1, len
87 IF (q(i,j1,l)<0.) THEN
88 ip = ip + 1
89
90 ! From above
91
92 qup = q(i, j1, nlaym1)
93 qly = -q(i, j1, l)
94 dup = min(qly, qup)
95 q(i, j1, nlaym1) = qup - dup
96 ! From "below" the surface.
97 sum = sum + qly - dup
98 q(i, j1, l) = 0.
99 END IF
100 END DO
101 end IF
102 end IF
103 end IF
104
105 IF (ip>imr) THEN
106 WRITE (6, *) 'IC=', ic, ' STEP=', nstep, ' Vertical filling pts=', ip
107 END IF
108
109 IF (sum>1.E-25) THEN
110 WRITE (6, *) ic, nstep, ' Mass source from the ground=', sum
111 END IF
112
113 END SUBROUTINE qckxyz

  ViewVC Help
Powered by ViewVC 1.1.21