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

Annotation of /trunk/Sources/dyn3d/PPM3d/fzppm.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 166 - (hide 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 guez 166 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