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

Contents of /trunk/Sources/dyn3d/PPM3d/fzppm.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: 5248 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 fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, &
2 wk1, wk2, wz2, delp, kord)
3 PARAMETER (kmax=150)
4 PARAMETER (r23=2./3., r3=1./3.)
5 REAL wz(imr, jnp, nlay), p(imr, jnp, nlay), dc(imr, jnp, nlay), &
6 wk1(imr, *), delp(imr, jnp, nlay), dq(imr, jnp, nlay), &
7 dqdt(imr, jnp, nlay)
8 ! Assuming JNP >= NLAY
9 REAL ar(imr, *), al(imr, *), a6(imr, *), flux(imr, *), wk2(imr, *), &
10 wz2(imr, *)
11
12 jmr = jnp - 1
13 imjm = imr*jnp
14 nlaym1 = nlay - 1
15
16 lmt = kord - 3
17
18 ! ****6***0*********0*********0*********0*********0*********0**********72
19 ! Compute DC for PPM
20 ! ****6***0*********0*********0*********0*********0*********0**********72
21
22 DO k = 1, nlaym1
23 DO i = 1, imjm
24 dqdt(i, 1, k) = p(i, 1, k+1) - p(i, 1, k)
25 END DO
26 END DO
27
28 DO k = 2, nlaym1
29 DO i = 1, imjm
30 c0 = delp(i, 1, k)/(delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
31 c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))
32 c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))
33 tmp = c0*(c1*dqdt(i,1,k)+c2*dqdt(i,1,k-1))
34 qmax = max(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - p(i, 1, k)
35 qmin = p(i, 1, k) - min(p(i,1,k-1), p(i,1,k), p(i,1,k+1))
36 dc(i, 1, k) = sign(min(abs(tmp),qmax,qmin), tmp)
37 END DO
38 END DO
39
40
41 ! ****6***0*********0*********0*********0*********0*********0**********72
42 ! Loop over latitudes (to save memory)
43 ! ****6***0*********0*********0*********0*********0*********0**********72
44
45 DO j = 1, jnp
46 IF ((j==2 .OR. j==jmr) .AND. j1/=2) GO TO 2000
47
48 DO k = 1, nlay
49 DO i = 1, imr
50 wz2(i, k) = wz(i, j, k)
51 wk1(i, k) = p(i, j, k)
52 wk2(i, k) = delp(i, j, k)
53 flux(i, k) = dc(i, j, k) !this flux is actually the monotone slope
54 END DO
55 END DO
56
57 ! ****6***0*********0*********0*********0*********0*********0**********72
58 ! Compute first guesses at cell interfaces
59 ! First guesses are required to be continuous.
60 ! ****6***0*********0*********0*********0*********0*********0**********72
61
62 ! three-cell parabolic subgrid distribution at model top
63 ! two-cell parabolic with zero gradient subgrid distribution
64 ! at the surface.
65
66 ! First guess top edge value
67 DO i = 1, imr
68 ! three-cell PPM
69 ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
70 a = 3.*(dqdt(i,j,2)-dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/(wk2(i,1)+wk2(i, &
71 2)))/((wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)))
72 b = 2.*dqdt(i, j, 1)/(wk2(i,1)+wk2(i,2)) - r23*a*(2.*wk2(i,1)+wk2(i,2))
73 al(i, 1) = wk1(i, 1) - wk2(i, 1)*(r3*a*wk2(i,1)+0.5*b)
74 al(i, 2) = wk2(i, 1)*(a*wk2(i,1)+b) + al(i, 1)
75
76 ! Check if change sign
77 IF (wk1(i,1)*al(i,1)<=0.) THEN
78 al(i, 1) = 0.
79 flux(i, 1) = 0.
80 ELSE
81 flux(i, 1) = wk1(i, 1) - al(i, 1)
82 END IF
83 END DO
84
85 ! Bottom
86 DO i = 1, imr
87 ! 2-cell PPM with zero gradient right at the surface
88
89 fct = dqdt(i, j, nlaym1)*wk2(i, nlay)**2/((wk2(i,nlay)+wk2(i, &
90 nlaym1))*(2.*wk2(i,nlay)+wk2(i,nlaym1)))
91 ar(i, nlay) = wk1(i, nlay) + fct
92 al(i, nlay) = wk1(i, nlay) - (fct+fct)
93 IF (wk1(i,nlay)*ar(i,nlay)<=0.) ar(i, nlay) = 0.
94 flux(i, nlay) = ar(i, nlay) - wk1(i, nlay)
95 END DO
96
97
98 ! ****6***0*********0*********0*********0*********0*********0**********72
99 ! 4th order interpolation in the interior.
100 ! ****6***0*********0*********0*********0*********0*********0**********72
101
102 DO k = 3, nlaym1
103 DO i = 1, imr
104 c1 = dqdt(i, j, k-1)*wk2(i, k-1)/(wk2(i,k-1)+wk2(i,k))
105 c2 = 2./(wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
106 a1 = (wk2(i,k-2)+wk2(i,k-1))/(2.*wk2(i,k-1)+wk2(i,k))
107 a2 = (wk2(i,k)+wk2(i,k+1))/(2.*wk2(i,k)+wk2(i,k-1))
108 al(i, k) = wk1(i, k-1) + c1 + c2*(wk2(i,k)*(c1*(a1-a2)+a2*flux(i, &
109 k-1))-wk2(i,k-1)*a1*flux(i,k))
110 END DO
111 END DO
112
113 DO i = 1, imr*nlaym1
114 ar(i, 1) = al(i, 2)
115 END DO
116
117 DO i = 1, imr*nlay
118 a6(i, 1) = 3.*(wk1(i,1)+wk1(i,1)-(al(i,1)+ar(i,1)))
119 END DO
120
121 ! ****6***0*********0*********0*********0*********0*********0**********72
122 ! Top & Bot always monotonic
123 CALL lmtppm(flux(1,1), a6(1,1), ar(1,1), al(1,1), wk1(1,1), imr, 0)
124 CALL lmtppm(flux(1,nlay), a6(1,nlay), ar(1,nlay), al(1,nlay), &
125 wk1(1,nlay), imr, 0)
126
127 ! Interior depending on KORD
128 IF (lmt<=2) CALL lmtppm(flux(1,2), a6(1,2), ar(1,2), al(1,2), wk1(1,2), &
129 imr*(nlay-2), lmt)
130
131 ! ****6***0*********0*********0*********0*********0*********0**********72
132
133 DO i = 1, imr*nlaym1
134 IF (wz2(i,1)>0.) THEN
135 cm = wz2(i, 1)/wk2(i, 1)
136 flux(i, 2) = ar(i, 1) + 0.5*cm*(al(i,1)-ar(i,1)+a6(i,1)*(1.-r23*cm))
137 ELSE
138 cp = wz2(i, 1)/wk2(i, 2)
139 flux(i, 2) = al(i, 2) + 0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp))
140 END IF
141 END DO
142
143 DO i = 1, imr*nlaym1
144 flux(i, 2) = wz2(i, 1)*flux(i, 2)
145 END DO
146
147 DO i = 1, imr
148 dq(i, j, 1) = dq(i, j, 1) - flux(i, 2)
149 dq(i, j, nlay) = dq(i, j, nlay) + flux(i, nlay)
150 END DO
151
152 DO k = 2, nlaym1
153 DO i = 1, imr
154 dq(i, j, k) = dq(i, j, k) + flux(i, k) - flux(i, k+1)
155 END DO
156 END DO
157 2000 END DO
158 RETURN
159 END SUBROUTINE fzppm
160
  ViewVC Help
Powered by ViewVC 1.1.21