sint.F
References to this file elsewhere.
1
2 SUBROUTINE SINT(XF, &
3 ims, ime, jms, jme, icmask , &
4 its, ite, jts, jte, nf, xstag, ystag )
5 IMPLICIT NONE
6 INTEGER ims, ime, jms, jme, &
7 its, ite, jts, jte
8
9 LOGICAL icmask( ims:ime, jms:jme )
10 LOGICAL xstag, ystag
11
12 INTEGER nf, ior
13 REAL one12, one24, ep
14 PARAMETER(one12=1./12.,one24=1./24.)
15 PARAMETER(ior=2)
16 !
17 REAL XF(ims:ime,jms:jme,NF)
18 !
19 REAL Y(ims:ime,jms:jme,-IOR:IOR), &
20 Z(ims:ime,jms:jme,-IOR:IOR), &
21 F(ims:ime,jms:jme,0:1)
22 !
23 INTEGER I,J,II,JJ,IIM
24 INTEGER N2STAR, N2END, N1STAR, N1END
25 !
26 DATA EP/ 1.E-10/
27 !
28 ! PARAMETER(NONOS=1)
29 ! PARAMETER(N1OS=N1*NONOS+1-NONOS,N2OS=N2*NONOS+1-NONOS)
30 !
31 REAL W(ims:ime,jms:jme),OV(ims:ime,jms:jme),UN(ims:ime,jms:jme)
32 REAL MXM(ims:ime,jms:jme),MN(ims:ime,jms:jme)
33 REAL FL(ims:ime,jms:jme,0:1)
34 REAL XIG(81), XJG(81) ! won't use but nine of these fellers.
35 INTEGER IFRST
36 integer rr
37 COMMON /DEPAR2/ XIG,XJG,IFRST
38 DATA IFRST /1/
39
40 REAL rioff, rjoff
41 !
42 REAL donor, y1, y2, a
43 DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
44 REAL tr4, ym1, y0, yp1, yp2
45 TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1)) &
46 -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0) &
47 -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
48 REAL pp, pn, x
49 PP(X)=AMAX1(0.,X)
50 PN(X)=AMIN1(0.,X)
51 !! XIG(I) = 1./3.-FLOAT(I-1)*1./3
52 !! XJG(J) = 1./3.-FLOAT(J-1)*1./3
53
54 rr = nint(sqrt(float(nf)))
55 !! write(6,*) ' nf, rr are ',nf,rr
56 !! IF ( IFRST .EQ. 1 ) THEN
57
58 rioff = 0
59 rjoff = 0
60 if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
61 if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
62
63 DO I=1,rr
64 DO J=1,rr
65 XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
66 XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
67 ENDDO
68 ENDDO
69 IFRST = 0
70
71 !! ENDIF
72
73 ! IF ( IFRST .EQ. 1 ) THEN
74 ! DO I=1,3
75 ! DO J=1,3
76 ! XIG(J+(I-1)*3)=1./3.-FLOAT(J-1)*1./3.
77 ! XJG(J+(I-1)*3)=1./3.-FLOAT(I-1)*1./3.
78 ! ENDDO
79 ! ENDDO
80 ! IFRST = 0
81 ! ENDIF
82 !
83 N2STAR = jts
84 N2END = jte
85 N1STAR = its
86 N1END = ite
87
88 DO 2000 IIM=1,NF
89 !
90 ! HERE STARTS RESIDUAL ADVECTION
91 !
92 DO 9000 JJ=N2STAR,N2END
93 DO 50 J=-IOR,IOR
94
95 DO 51 I=-IOR,IOR
96 DO 511 II=N1STAR,N1END
97 IF ( icmask(II,JJ) ) Y(II,JJ,I)=XF(II+I,JJ+J,IIM)
98 511 CONTINUE
99 51 CONTINUE
100
101 DO 811 II=N1STAR,N1END
102 IF ( icmask(II,JJ) ) THEN
103 FL(II,JJ,0)=DONOR(Y(II,JJ,-1),Y(II,JJ,0),XIG(IIM))
104 FL(II,JJ,1)=DONOR(Y(II,JJ,0),Y(II,JJ,1),XIG(IIM))
105 ENDIF
106 811 CONTINUE
107 DO 812 II=N1STAR,N1END
108 IF ( icmask(II,JJ) ) W(II,JJ)=Y(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
109 812 CONTINUE
110 DO 813 II=N1STAR,N1END
111 IF ( icmask(II,JJ) ) THEN
112 MXM(II,JJ)= &
113 AMAX1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1), &
114 W(II,JJ))
115 MN(II,JJ)=AMIN1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),W(II,JJ))
116 ENDIF
117 813 CONTINUE
118 DO 312 II=N1STAR,N1END
119 IF ( icmask(II,JJ) ) THEN
120 F(II,JJ,0)= &
121 TR4(Y(II,JJ,-2),Y(II,JJ,-1),Y(II,JJ,0), &
122 Y(II,JJ,1),XIG(IIM))
123 F(II,JJ,1)= &
124 TR4(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),Y(II,JJ,2),&
125 XIG(IIM))
126 ENDIF
127 312 CONTINUE
128 DO 822 II=N1STAR,N1END
129 IF ( icmask(II,JJ) ) THEN
130 F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
131 F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
132 ENDIF
133 822 CONTINUE
134 DO 823 II=N1STAR,N1END
135 IF ( icmask(II,JJ) ) THEN
136 OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
137 PP(F(II,JJ,0))+EP)
138 UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))- &
139 PN(F(II,JJ,0))+EP)
140 ENDIF
141 823 CONTINUE
142 DO 824 II=N1STAR,N1END
143 IF ( icmask(II,JJ) ) THEN
144 F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+ &
145 PN(F(II,JJ,0))*AMIN1(1.,UN(II,JJ))
146 F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+ &
147 PN(F(II,JJ,1))*AMIN1(1.,OV(II,JJ))
148 ENDIF
149 824 CONTINUE
150 DO 825 II=N1STAR,N1END
151 IF ( icmask(II,JJ) ) THEN
152 Y(II,JJ,0)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
153 ENDIF
154 825 CONTINUE
155 DO 361 II=N1STAR,N1END
156 IF ( icmask(II,JJ) ) Z(II,JJ,J)=Y(II,JJ,0)
157 361 CONTINUE
158 !
159 ! END IF FIRST J LOOP
160 !
161 8000 CONTINUE
162 50 CONTINUE
163
164 DO 911 II=N1STAR,N1END
165 IF ( icmask(II,JJ) ) THEN
166 FL(II,JJ,0)=DONOR(Z(II,JJ,-1),Z(II,JJ,0),XJG(IIM))
167 FL(II,JJ,1)=DONOR(Z(II,JJ,0),Z(II,JJ,1),XJG(IIM))
168 ENDIF
169 911 CONTINUE
170 DO 912 II=N1STAR,N1END
171 IF ( icmask(II,JJ) ) W(II,JJ)=Z(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
172 912 CONTINUE
173 DO 913 II=N1STAR,N1END
174 IF ( icmask(II,JJ) ) THEN
175 MXM(II,JJ)=AMAX1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
176 MN(II,JJ)=AMIN1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
177 ENDIF
178 913 CONTINUE
179 DO 412 II=N1STAR,N1END
180 IF ( icmask(II,JJ) ) THEN
181 F(II,JJ,0)= &
182 TR4(Z(II,JJ,-2),Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1)&
183 ,XJG(IIM))
184 F(II,JJ,1)= &
185 TR4(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),Z(II,JJ,2), &
186 XJG(IIM))
187 ENDIF
188 412 CONTINUE
189 DO 922 II=N1STAR,N1END
190 IF ( icmask(II,JJ) ) THEN
191 F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
192 F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
193 ENDIF
194 922 CONTINUE
195 DO 923 II=N1STAR,N1END
196 IF ( icmask(II,JJ) ) THEN
197 OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
198 PP(F(II,JJ,0))+EP)
199 UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-PN(F(II,JJ,0))+ &
200 EP)
201 ENDIF
202 923 CONTINUE
203 DO 924 II=N1STAR,N1END
204 IF ( icmask(II,JJ) ) THEN
205 F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+PN(F(II,JJ,0)) &
206 *AMIN1(1.,UN(II,JJ))
207 F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+PN(F(II,JJ,1)) &
208 *AMIN1(1.,OV(II,JJ))
209 ENDIF
210 924 CONTINUE
211 9000 CONTINUE
212 DO 925 JJ=N2STAR,N2END
213 DO 925 II=N1STAR,N1END
214 IF ( icmask(II,JJ) ) XF(II,JJ,IIM)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
215 925 CONTINUE
216
217 !
218 2000 CONTINUE
219 RETURN
220 END
221
222 ! Version of sint that replaces mask with detailed ranges for avoiding boundaries
223 ! may help performance by getting the conditionals out of innner loops
224
225 SUBROUTINE SINTB(XF1, XF , &
226 ims, ime, jms, jme, icmask , &
227 its, ite, jts, jte, nf, xstag, ystag )
228 IMPLICIT NONE
229 INTEGER ims, ime, jms, jme, &
230 its, ite, jts, jte
231
232 LOGICAL icmask( ims:ime, jms:jme )
233 LOGICAL xstag, ystag
234
235 INTEGER nf, ior
236 REAL one12, one24, ep
237 PARAMETER(one12=1./12.,one24=1./24.)
238 PARAMETER(ior=2)
239 !
240 REAL XF(ims:ime,jms:jme,NF)
241 REAL XF1(ims:ime,jms:jme,NF)
242 !
243 REAL Y(ims:ime,jms:jme,-IOR:IOR), &
244 Z(ims:ime,jms:jme,-IOR:IOR), &
245 F(ims:ime,jms:jme,0:1)
246 !
247 INTEGER I,J,II,JJ,IIM
248 INTEGER N2STAR, N2END, N1STAR, N1END
249 !
250 DATA EP/ 1.E-10/
251 !
252 ! PARAMETER(NONOS=1)
253 ! PARAMETER(N1OS=N1*NONOS+1-NONOS,N2OS=N2*NONOS+1-NONOS)
254 !
255 REAL W(ims:ime,jms:jme),OV(ims:ime,jms:jme),UN(ims:ime,jms:jme)
256 REAL MXM(ims:ime,jms:jme),MN(ims:ime,jms:jme)
257 REAL FL(ims:ime,jms:jme,0:1)
258 REAL XIG(81), XJG(81) ! won't use but nine of these fellers.
259 INTEGER IFRST
260 integer rr
261 COMMON /DEPAR2B/ XIG,XJG,IFRST
262 DATA IFRST /1/
263
264 REAL rioff, rjoff
265 !
266 REAL donor, y1, y2, a
267 DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
268 REAL tr4, ym1, y0, yp1, yp2
269 TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1)) &
270 -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0) &
271 -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
272 REAL pp, pn, x
273 PP(X)=AMAX1(0.,X)
274 PN(X)=AMIN1(0.,X)
275
276 rr = nint(sqrt(float(nf)))
277
278 rioff = 0
279 rjoff = 0
280 if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
281 if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
282
283 DO I=1,rr
284 DO J=1,rr
285 XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
286 XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
287 ENDDO
288 ENDDO
289 IFRST = 0
290
291 !! ENDIF
292
293 ! IF ( IFRST .EQ. 1 ) THEN
294 ! DO I=1,3
295 ! DO J=1,3
296 ! XIG(J+(I-1)*3)=1./3.-FLOAT(J-1)*1./3.
297 ! XJG(J+(I-1)*3)=1./3.-FLOAT(I-1)*1./3.
298 ! ENDDO
299 ! ENDDO
300 ! IFRST = 0
301 ! ENDIF
302 !
303 N2STAR = jts
304 N2END = jte
305 N1STAR = its
306 N1END = ite
307
308 DO 2000 IIM=1,NF
309 !
310 ! HERE STARTS RESIDUAL ADVECTION
311 !
312 DO 9000 JJ=N2STAR,N2END
313 !cdir unroll=5
314 DO 50 J=-IOR,IOR
315
316 !cdir unroll=5
317 DO 51 I=-IOR,IOR
318 DO 511 II=N1STAR,N1END
319 Y(II,JJ,I)=XF1(II+I,JJ+J,IIM)
320 511 CONTINUE
321 51 CONTINUE
322
323 DO 811 II=N1STAR,N1END
324 FL(II,JJ,0)=DONOR(Y(II,JJ,-1),Y(II,JJ,0),XIG(IIM))
325 FL(II,JJ,1)=DONOR(Y(II,JJ,0),Y(II,JJ,1),XIG(IIM))
326 811 CONTINUE
327 DO 812 II=N1STAR,N1END
328 W(II,JJ)=Y(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
329 812 CONTINUE
330 DO 813 II=N1STAR,N1END
331 MXM(II,JJ)= &
332 AMAX1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1), &
333 W(II,JJ))
334 MN(II,JJ)=AMIN1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),W(II,JJ))
335 813 CONTINUE
336 DO 312 II=N1STAR,N1END
337 F(II,JJ,0)= &
338 TR4(Y(II,JJ,-2),Y(II,JJ,-1),Y(II,JJ,0), &
339 Y(II,JJ,1),XIG(IIM))
340 F(II,JJ,1)= &
341 TR4(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),Y(II,JJ,2),&
342 XIG(IIM))
343 312 CONTINUE
344 DO 822 II=N1STAR,N1END
345 F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
346 F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
347 822 CONTINUE
348 DO 823 II=N1STAR,N1END
349 OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
350 PP(F(II,JJ,0))+EP)
351 UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))- &
352 PN(F(II,JJ,0))+EP)
353 823 CONTINUE
354 DO 824 II=N1STAR,N1END
355 F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+ &
356 PN(F(II,JJ,0))*AMIN1(1.,UN(II,JJ))
357 F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+ &
358 PN(F(II,JJ,1))*AMIN1(1.,OV(II,JJ))
359 824 CONTINUE
360 DO 825 II=N1STAR,N1END
361 Y(II,JJ,0)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
362 825 CONTINUE
363 DO 361 II=N1STAR,N1END
364 Z(II,JJ,J)=Y(II,JJ,0)
365 361 CONTINUE
366 !
367 ! END IF FIRST J LOOP
368 !
369 8000 CONTINUE
370 50 CONTINUE
371
372 DO 911 II=N1STAR,N1END
373 FL(II,JJ,0)=DONOR(Z(II,JJ,-1),Z(II,JJ,0),XJG(IIM))
374 FL(II,JJ,1)=DONOR(Z(II,JJ,0),Z(II,JJ,1),XJG(IIM))
375 911 CONTINUE
376 DO 912 II=N1STAR,N1END
377 W(II,JJ)=Z(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
378 912 CONTINUE
379 DO 913 II=N1STAR,N1END
380 MXM(II,JJ)=AMAX1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
381 MN(II,JJ)=AMIN1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
382 913 CONTINUE
383 DO 412 II=N1STAR,N1END
384 F(II,JJ,0)= &
385 TR4(Z(II,JJ,-2),Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1)&
386 ,XJG(IIM))
387 F(II,JJ,1)= &
388 TR4(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),Z(II,JJ,2), &
389 XJG(IIM))
390 412 CONTINUE
391 DO 922 II=N1STAR,N1END
392 F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
393 F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
394 922 CONTINUE
395 DO 923 II=N1STAR,N1END
396 OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
397 PP(F(II,JJ,0))+EP)
398 UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-PN(F(II,JJ,0))+ &
399 EP)
400 923 CONTINUE
401 DO 924 II=N1STAR,N1END
402 F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+PN(F(II,JJ,0)) &
403 *AMIN1(1.,UN(II,JJ))
404 F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+PN(F(II,JJ,1)) &
405 *AMIN1(1.,OV(II,JJ))
406 924 CONTINUE
407 9000 CONTINUE
408 DO 925 JJ=N2STAR,N2END
409 DO 925 II=N1STAR,N1END
410 XF(II,JJ,IIM)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
411 925 CONTINUE
412
413 !
414 2000 CONTINUE
415 RETURN
416 END
417
418