67
68
73
74 implicit none
75
76
77
78 integer, intent(in) :: ng, LBi, UBi, LBj, UBj, LBk, UBk
79 integer, intent(in) :: Lstr, Lend, itime, ifield, isBval, gtype
80
81 logical, intent(in) :: maskit
82 logical, intent(in) :: my_thread(Lstr:Lend)
83 logical, intent(in) :: bounded(Nfloats(ng))
84
85 real(dp), intent(in) :: Fspval
86
87 real(r8), intent(in) :: nudg(Lstr:Lend)
88
89 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
90 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
91# ifdef SOLVE3D
92 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
93# endif
94# ifdef MASKING
95 real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
96# endif
97 real(r8), intent(in) :: A(LBi:UBi,LBj:UBj,LBk:UBk)
98
99 real(r8), intent(inout) :: track(NFV(ng),0:NFT,Nfloats(ng))
100
101
102
103 logical :: Irvar, Iuvar, Jrvar, Jvvar, Krvar, Kwvar, Lmask
104 logical :: halo
105
106 integer :: Ir, Iu, Jr, Jv, Kr, Kw, l
107 integer :: i1, i2, j1, j2, k1, k2, khm, khp, vtype
108
109 real(r8) :: p1, p2, q1, q2, r1, r2
110 real(r8) :: s111, s211, s121, s221, s112, s212, s122, s222
111 real(r8) :: t111, t211, t121, t221, t112, t212, t122, t222
112
113# ifdef MASKING
114 integer :: Irn, Irnm1, Irnp1, Jrn, Jrnm1, Jrnp1
115
116 real(r8) :: cff1, cff2, cff3
117# endif
118
119
120
121
122
123
124
125 vtype=abs(gtype)
137
138
139
140# ifdef MASKING
141 lmask=maskit
142# else
143 lmask=.false.
144# endif
145
146
147
148 IF (gtype.ge.0) THEN
149 s111=1.0_r8
150 s121=1.0_r8
151 s211=1.0_r8
152 s221=1.0_r8
153 s112=1.0_r8
154 s122=1.0_r8
155 s212=1.0_r8
156 s222=1.0_r8
157 t111=1.0_r8
158 t121=1.0_r8
159 t211=1.0_r8
160 t221=1.0_r8
161 t112=1.0_r8
162 t122=1.0_r8
163 t212=1.0_r8
164 t222=1.0_r8
165 END IF
166
167
168
169
170
171 DO l=lstr,lend
172 IF (my_thread(l)) THEN
173 IF (.not.bounded(l)) THEN
174 track(ifield,itime,l)=fspval
175 ELSE
176
177
178
179 IF (krvar) THEN
180 kr=int(track(
izgrd,itime,l)+0.5_r8)
181 k1=min(max(kr ,1),
n(ng))
182 k2=min(max(kr+1,1),
n(ng))
183 r2=real(k2-k1,r8)*(track(
izgrd,itime,l)+ &
184 & 0.5_r8-real(k1,r8))
185 ELSE IF (kwvar) THEN
186 kw=int(track(
izgrd,itime,l))
187 k1=min(max(kw ,0),
n(ng))
188 k2=min(max(kw+1,0),
n(ng))
189 r2=real(k2-k1,r8)*(track(
izgrd,itime,l)-real(k1,r8))
190 ELSE
191 k1=1
192 k2=1
193 r2=0.0_r8
194 END IF
195 r1=1.0_r8-r2
196
197
198
199
200
201 IF (irvar.and.jrvar) THEN
202
203 ir=int(track(
ixgrd,itime,l))
204 jr=int(track(
iygrd,itime,l))
205
206 i1=min(max(ir ,0),
lm(ng)+1)
207 i2=min(max(ir+1,1),
lm(ng)+1)
208 j1=min(max(jr ,0),
mm(ng)+1)
209 j2=min(max(jr+1,1),
mm(ng)+1)
210
211 p2=real(i2-i1,r8)*(track(
ixgrd,itime,l)-real(i1,r8))
212 q2=real(j2-j1,r8)*(track(
iygrd,itime,l)-real(j1,r8))
213 p1=1.0_r8-p2
214 q1=1.0_r8-q2
215# ifdef SOLVE3D
216
217 IF (gtype.eq.-
w3dvar)
THEN
218 khm=min(max(k1 ,1),
n(ng))
219 khp=min(max(k1+1,1),
n(ng))
220 s111=2.0_r8*pm(i1,j1)*pn(i1,j1)/ &
221 & (hz(i1,j1,khm)+hz(i1,j1,khp))
222 s211=2.0_r8*pm(i2,j1)*pn(i2,j1)/ &
223 & (hz(i2,j1,khm)+hz(i2,j1,khp))
224 s121=2.0_r8*pm(i1,j2)*pn(i1,j2)/ &
225 & (hz(i1,j2,khm)+hz(i1,j2,khp))
226 s221=2.0_r8*pm(i2,j2)*pn(i2,j2)/ &
227 & (hz(i2,j2,khm)+hz(i2,j2,khp))
228 t111=2.0_r8/(hz(i1,j1,khm)+hz(i1,j1,khp))
229 t211=2.0_r8/(hz(i2,j1,khm)+hz(i2,j1,khp))
230 t121=2.0_r8/(hz(i1,j2,khm)+hz(i1,j2,khp))
231 t221=2.0_r8/(hz(i2,j2,khm)+hz(i2,j2,khp))
232 khm=min(max(k2 ,1),
n(ng))
233 khp=min(max(k2+1,1),
n(ng))
234 s112=2.0_r8*pm(i1,j1)*pn(i1,j1)/ &
235 & (hz(i1,j1,khm)+hz(i1,j1,khp))
236 s212=2.0_r8*pm(i2,j1)*pn(i2,j1)/ &
237 & (hz(i2,j1,khm)+hz(i2,j1,khp))
238 s122=2.0_r8*pm(i1,j2)*pn(i1,j2)/ &
239 & (hz(i1,j2,khm)+hz(i1,j2,khp))
240 s222=2.0_r8*pm(i2,j2)*pn(i2,j2)/ &
241 & (hz(i2,j2,khm)+hz(i2,j2,khp))
242 t112=2.0_r8/(hz(i1,j1,khm)+hz(i1,j1,khp))
243 t212=2.0_r8/(hz(i2,j1,khm)+hz(i2,j1,khp))
244 t122=2.0_r8/(hz(i1,j2,khm)+hz(i1,j2,khp))
245 t222=2.0_r8/(hz(i2,j2,khm)+hz(i2,j2,khp))
246 END IF
247# endif
248
249 IF (lmask) THEN
250# ifdef MASKING
251 cff1=p1*q1*r1*amask(i1,j1)+ &
252 & p2*q1*r1*amask(i2,j1)+ &
253 & p1*q2*r1*amask(i1,j2)+ &
254 & p2*q2*r1*amask(i2,j2)+ &
255 & p1*q1*r2*amask(i1,j1)+ &
256 & p2*q1*r2*amask(i2,j1)+ &
257 & p1*q2*r2*amask(i1,j2)+ &
258 & p2*q2*r2*amask(i2,j2)
259 IF (cff1.gt.0.0_r8) THEN
260 cff2=p1*q1*r1*amask(i1,j1)*s111*a(i1,j1,k1)+ &
261 & p2*q1*r1*amask(i2,j1)*s211*a(i2,j1,k1)+ &
262 & p1*q2*r1*amask(i1,j2)*s121*a(i1,j2,k1)+ &
263 & p2*q2*r1*amask(i2,j2)*s221*a(i2,j2,k1)+ &
264 & p1*q1*r2*amask(i1,j1)*s112*a(i1,j1,k2)+ &
265 & p2*q1*r2*amask(i2,j1)*s212*a(i2,j1,k2)+ &
266 & p1*q2*r2*amask(i1,j2)*s122*a(i1,j2,k2)+ &
267 & p2*q2*r2*amask(i2,j2)*s222*a(i2,j2,k2)
268 cff3=(p1*q1*r1*amask(i1,j1)*t111+ &
269 & p2*q1*r1*amask(i2,j1)*t211+ &
270 & p1*q2*r1*amask(i1,j2)*t121+ &
271 & p2*q2*r1*amask(i2,j2)*t221+ &
272 & p1*q1*r2*amask(i1,j1)*t112+ &
273 & p2*q1*r2*amask(i2,j1)*t212+ &
274 & p1*q2*r2*amask(i1,j2)*t122+ &
275 & p2*q2*r2*amask(i2,j2)*t222)*nudg(l)
276 track(ifield,itime,l)=cff2/cff1+cff3
277 ELSE
278 track(ifield,itime,l)=0.0_r8
279 END IF
280# endif
281 ELSE
282 track(ifield,itime,l)=p1*q1*r1*s111*a(i1,j1,k1)+ &
283 & p2*q1*r1*s211*a(i2,j1,k1)+ &
284 & p1*q2*r1*s121*a(i1,j2,k1)+ &
285 & p2*q2*r1*s221*a(i2,j2,k1)+ &
286 & p1*q1*r2*s112*a(i1,j1,k2)+ &
287 & p2*q1*r2*s212*a(i2,j1,k2)+ &
288 & p1*q2*r2*s122*a(i1,j2,k2)+ &
289 & p2*q2*r2*s222*a(i2,j2,k2)+ &
290 & (p1*q1*r1*t111+ &
291 & p2*q1*r1*t211+ &
292 & p1*q2*r1*t121+ &
293 & p2*q2*r1*t221+ &
294 & p1*q1*r2*t112+ &
295 & p2*q1*r2*t212+ &
296 & p1*q2*r2*t122+ &
297 & p2*q2*r2*t222)*nudg(l)
298 END IF
299
300
301
302
303
304 ELSE
305 ir=int(track(
ixgrd,itime,l))
306 jr=int(track(
iygrd,itime,l))
307 iu=int(track(
ixgrd,itime,l)+0.5_r8)
308 jv=int(track(
iygrd,itime,l)+0.5_r8)
309
310 halo=.false.
311# ifdef MASKING
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326 IF (lmask) THEN
327 irn=nint(track(
ixgrd,itime,l))
328 jrn=nint(track(
iygrd,itime,l))
330 IF (irn.ge.
lm(ng))
THEN
332 ELSE
333 irnp1=irn+1
334 END IF
335 IF (irn.le.1) THEN
337 ELSE
338 irnm1=irn-1
339 END IF
340 ELSE
341 irnm1=irn-1
342 irnp1=irn+1
343 END IF
345 IF (jrn.ge.
mm(ng))
THEN
347 ELSE
348 jrnp1=jrn+1
349 END IF
350 IF (jrn.le.1) THEN
352 ELSE
353 jrnm1=jrn-1
354 END IF
355 ELSE
356 jrnm1=jrn-1
357 jrnp1=jrn+1
358 END IF
359 IF (amask(irn,jrn).lt.0.5_r8) THEN
360 halo=.true.
361 ELSE IF ((ir.lt.irn).and. &
362 & (amask(irn-1,jrn).lt.0.5_r8)) THEN
363 halo=.true.
364 ELSE IF ((ir.eq.irn).and. &
365 & (amask(irn+1,jrn).lt.0.5_r8)) THEN
366 halo=.true.
367 ELSE IF ((jr.lt.jrn).and. &
368 & (amask(irn,jrn-1).lt.0.5_r8)) THEN
369 halo=.true.
370 ELSE IF ((jr.eq.jrn).and. &
371 & (amask(irn,jrn+1).lt.0.5_r8)) THEN
372 halo=.true.
373 ELSE IF ((ir.lt.irn).and.(jr.lt.jrn).and. &
374 & (amask(irn-1,jrn-1).lt.0.5_r8)) THEN
375 halo=.true.
376 ELSE IF ((ir.eq.irn).and.(jr.lt.jrn).and. &
377 & (amask(irn+1,jrn-1).lt.0.5_r8)) THEN
378 halo=.true.
379 ELSE IF ((ir.lt.irn).and.(jr.eq.jrn).and. &
380 & (amask(irn-1,jrn+1).lt.0.5_r8)) THEN
381 halo=.true.
382 ELSE IF ((ir.eq.irn).and.(jr.eq.jrn).and. &
383 & (amask(irn+1,jrn+1).lt.0.5_r8)) THEN
384 halo=.true.
385 END IF
386 END IF
387# endif
388
389
390
391
392
393 IF (iuvar) THEN
394 IF (halo) THEN
395# ifdef MASKING
396
397
398
399
400
401
402 i1=min(max(iu ,1),
lm(ng)+1)
403 i2=min(max(iu+1,1),
lm(ng)+1)
404 j1=jrn
405
406 p2=real(i2-i1,r8)* &
407 & (track(
ixgrd,itime,l)-real(i1,r8)+0.5_r8)
408 p1=1.0_r8-p2
409 q1=1.0_r8
410
411 IF (gtype.lt.0) THEN
412 s111=0.5_r8*(pm(i1-1,j1)+pm(i1,j1))
413 s211=0.5_r8*(pm(i2-1,j1)+pm(i2,j1))
414 s112=s111
415 s212=s112
416 END IF
417
418 track(ifield,itime,l)=p1*q1*r1*s111*a(i1,j1,k1)+ &
419 & p2*q1*r1*s211*a(i2,j1,k1)+ &
420 & p1*q1*r2*s112*a(i1,j1,k2)+ &
421 & p2*q1*r2*s212*a(i2,j1,k2)+ &
422 & nudg(l)
423# endif
424 ELSE
425
426
427
428 i1=min(max(iu ,1),
lm(ng)+1)
429 i2=min(max(iu+1,1),
lm(ng)+1)
430 j1=min(max(jr ,0),
mm(ng)+1)
431 j2=min(max(jr+1,0),
mm(ng)+1)
432
433 p2=real(i2-i1,r8)* &
434 & (track(
ixgrd,itime,l)-real(i1,r8)+0.5_r8)
435 q2=real(j2-j1,r8)* &
436 & (track(
iygrd,itime,l)-real(j1,r8))
437 p1=1.0_r8-p2
438 q1=1.0_r8-q2
439
440 IF (gtype.lt.0) THEN
441 s111=0.5_r8*(pm(i1-1,j1)+pm(i1,j1))
442 s211=0.5_r8*(pm(i2-1,j1)+pm(i2,j1))
443 s121=0.5_r8*(pm(i1-1,j2)+pm(i1,j2))
444 s221=0.5_r8*(pm(i2-1,j2)+pm(i2,j2))
445 s112=s111
446 s212=s112
447 s122=s121
448 s222=s221
449 END IF
450
451 track(ifield,itime,l)=p1*q1*r1*s111*a(i1,j1,k1)+ &
452 & p2*q1*r1*s211*a(i2,j1,k1)+ &
453 & p1*q2*r1*s121*a(i1,j2,k1)+ &
454 & p2*q2*r1*s221*a(i2,j2,k1)+ &
455 & p1*q1*r2*s112*a(i1,j1,k2)+ &
456 & p2*q1*r2*s212*a(i2,j1,k2)+ &
457 & p1*q2*r2*s122*a(i1,j2,k2)+ &
458 & p2*q2*r2*s222*a(i2,j2,k2)+ &
459 & nudg(l)
460 END IF
461
462
463
464
465
466 ELSE IF (jvvar) THEN
467 IF (halo) THEN
468# ifdef MASKING
469
470
471
472
473
474
475 i1=irn
476 j1=min(max(jv ,1),
mm(ng)+1)
477 j2=min(max(jv+1,1),
mm(ng)+1)
478
479 q2=real(j2-j1,r8)* &
480 & (track(
iygrd,itime,l)-real(j1,r8)+0.5_r8)
481 p1=1.0_r8
482 q1=1.0_r8-q2
483
484 IF (gtype.lt.0) THEN
485 s111=0.5_r8*(pn(i1,j1-1)+pn(i1,j1))
486 s121=0.5_r8*(pn(i1,j2-1)+pn(i1,j2))
487 s112=s111
488 s122=s121
489 END IF
490
491 track(ifield,itime,l)=p1*q1*r1*s111*a(i1,j1,k1)+ &
492 & p1*q2*r1*s121*a(i1,j2,k1)+ &
493 & p1*q1*r2*s112*a(i1,j1,k2)+ &
494 & p1*q2*r2*s122*a(i1,j2,k2)+ &
495 & nudg(l)
496# endif
497 ELSE
498
499
500
501 i1=min(max(ir ,0),
lm(ng)+1)
502 i2=min(max(ir+1,1),
lm(ng)+1)
503 j1=min(max(jv ,1),
mm(ng)+1)
504 j2=min(max(jv+1,1),
mm(ng)+1)
505
506 p2=real(i2-i1,r8)* &
507 & (track(
ixgrd,itime,l)-real(i1,r8))
508 q2=real(j2-j1,r8)* &
509 & (track(
iygrd,itime,l)-real(j1,r8)+0.5_r8)
510 p1=1.0_r8-p2
511 q1=1.0_r8-q2
512
513 IF (gtype.lt.0) THEN
514 s111=0.5_r8*(pn(i1,j1-1)+pn(i1,j1))
515 s211=0.5_r8*(pn(i2,j1-1)+pn(i2,j1))
516 s121=0.5_r8*(pn(i1,j2-1)+pn(i1,j2))
517 s221=0.5_r8*(pn(i2,j2-1)+pn(i2,j2))
518 s112=s111
519 s212=s112
520 s122=s121
521 s222=s221
522 END IF
523
524 track(ifield,itime,l)=p1*q1*r1*s111*a(i1,j1,k1)+ &
525 & p2*q1*r1*s211*a(i2,j1,k1)+ &
526 & p1*q2*r1*s121*a(i1,j2,k1)+ &
527 & p2*q2*r1*s221*a(i2,j2,k1)+ &
528 & p1*q1*r2*s112*a(i1,j1,k2)+ &
529 & p2*q1*r2*s212*a(i2,j1,k2)+ &
530 & p1*q2*r2*s122*a(i1,j2,k2)+ &
531 & p2*q2*r2*s222*a(i2,j2,k2)+ &
532 & nudg(l)
533 END IF
534 END IF
535 END IF
536 END IF
537 END IF
538 END DO
539
540 RETURN
integer, dimension(:), allocatable n
integer, parameter r3dvar
integer, parameter u3dvar
integer, dimension(:), allocatable lm
integer, parameter u2dvar
integer, parameter w3dvar
integer, dimension(:), allocatable mm
integer, parameter r2dvar
integer, parameter v2dvar
integer, parameter v3dvar
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic