ROMS
Loading...
Searching...
No Matches
extract_sta.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE extract_sta_mod
3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! This routine extracts field at the requested (Xpos,Ypos,Zpos) !
12! positions. The extraction is done using linear interpolation. !
13! The (Xpos,Ypos) positions are in fractional grid coordinates. !
14! Zpos is in fractional grid coordinates (Zpos >= 0) or actual !
15! depths (Zpos < 0), if applicable. !
16! !
17! On Input: !
18! !
19! ng Nested grid number. !
20! model Calling model identifier. !
21! Cgrid Switch to interpolate at native C-grid (TRUE) or to !
22! interpolate at RHO-points (FALSE). !
23! ifield Field ID. !
24! gtype Grid type. !
25! LBi I-dimension Lower bound. !
26! UBi I-dimension Upper bound. !
27! LBj J-dimension Lower bound. !
28! UBj J-dimension Upper bound. !
29! LBk K-dimension Lower bound, if any. Otherwise, a value !
30! of one is expected. !
31! LBk K-dimension Upper bound, if any. Otherwise, a value !
32! of one is expected. !
33! UBk K-dimension Upper bound. !
34! Ascl Factor to scale field after extraction. !
35! A Tile array (2D or 3D) to process. !
36! Npos Number of values to extract. !
37! Xpos X-extraction positions (grid coordinates). !
38! Ypos Y-extraction positions (grid coordinates). !
39! Zpos Z-extraction positions (grid coordinates or depth). !
40! !
41! On Output: !
42! !
43! Apos Extracted values. !
44! !
45! Note: !
46! !
47! Starting F95 zero values can be signed (-0 or +0) following the !
48! IEEE 754 floating point standard. This can be advantageous in !
49! some computations but not here when "Ascl" is negative and "Apos" !
50! is zero. This will produce different output files in serial !
51! and distributed memory applications. Since comparing serial and !
52! parallel output is essential for tracking parallel partition !
53! bugs, "positive zero" is enforced. !
54! !
55!=======================================================================
56!
57 implicit none
58
59 CONTAINS
60!
61!***********************************************************************
62 SUBROUTINE extract_sta2d (ng, model, Cgrid, ifield, gtype, &
63 & LBi, UBi, LBj, UBj, Ascl, A, &
64 & Npos, Xpos, Ypos, Apos)
65!***********************************************************************
66!
67 USE mod_param
68 USE mod_parallel
69 USE mod_grid
70 USE mod_ncparam
71 USE mod_scalars
72
73# ifdef DISTRIBUTE
74!
75 USE distribute_mod, ONLY : mp_collect
76# endif
77!
78! Imported variable declarations.
79!
80 logical, intent(in) :: Cgrid
81!
82 integer, intent(in) :: ng, model, ifield, gtype, Npos
83 integer, intent(in) :: LBi, UBi, LBj, UBj
84 real(dp), intent(in) :: Ascl
85
86#ifdef ASSUMED_SHAPE
87 real(r8), intent(in) :: A(LBi:,LBj:)
88 real(r8), intent(in) :: Xpos(:), Ypos(:)
89 real(r8), intent(out) :: Apos(Npos)
90#else
91 real(r8), intent(in) :: A(LBi:UBi,LBj:UBj)
92 real(r8), intent(in) :: Xpos(Npos), Ypos(Npos)
93 real(r8), intent(out) :: Apos(Npos)
94#endif
95!
96! Local variable declarations.
97!
98 integer :: i1, i2, j1, j2, np
99
100 real(r8), parameter :: Aspv = 0.0_r8
101
102 real(r8) :: Xmin, Xmax, Ymin, Ymax
103 real(r8) :: Xgrd, Xoff, Ygrd, Yoff
104 real(r8) :: p1, p2, q1, q2, r1, r2, wsum
105 real(r8) :: w111, w211, w121, w221
106
107 real(r8), dimension(Npos) :: bounded
108!
109!-----------------------------------------------------------------------
110! Interpolate from 2D field at RHO-points.
111!-----------------------------------------------------------------------
112!
113 IF (gtype.eq.r2dvar) THEN
114 xmin=rxmin(ng)
115 xmax=rxmax(ng)
116 ymin=rymin(ng)
117 ymax=rymax(ng)
118 DO np=1,npos
119 xgrd=xpos(np)
120 ygrd=ypos(np)
121 bounded(np)=0.0_r8
122 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
123 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
124 i1=int(xgrd)
125 j1=int(ygrd)
126 i2=i1+1
127 j2=j1+1
128 IF (i2.gt.lm(ng)+1) THEN
129 i2=i1 ! station at the eastern boundary
130 END IF
131 IF (j2.gt.mm(ng)+1) THEN
132 j2=j1 ! station at the northern boundary
133 END IF
134 bounded(np)=1.0_r8
135 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
136 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
137 p1=1.0_r8-p2
138 q1=1.0_r8-q2
139 w111=p1*q1
140 w211=p2*q1
141 w121=p1*q2
142 w221=p2*q2
143#ifdef MASKING
144 w111=w111*grid(ng)%rmask(i1,j1)
145 w211=w211*grid(ng)%rmask(i2,j1)
146 w121=w121*grid(ng)%rmask(i1,j2)
147 w221=w221*grid(ng)%rmask(i2,j2)
148 wsum=w111+w211+w121+w221
149 IF (wsum.gt.0.0_r8) THEN
150 wsum=1.0_r8/wsum
151 w111=w111*wsum
152 w211=w211*wsum
153 w121=w121*wsum
154 w221=w221*wsum
155 ELSE
156 bounded(np)=0.0_r8
157 ENDIF
158#endif
159 apos(np)=ascl*(w111*a(i1,j1)+ &
160 & w211*a(i2,j1)+ &
161 & w121*a(i1,j2)+ &
162 & w221*a(i2,j2))
163 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
164 ELSE ! zero
165 apos(np)=aspv
166 END IF
167 END DO
168!
169!-----------------------------------------------------------------------
170! Interpolate from 2D field at U-points.
171!-----------------------------------------------------------------------
172!
173 ELSE IF (gtype.eq.u2dvar) THEN
174 IF (cgrid) THEN
175 xmin=uxmin(ng)+0.5_r8
176 xmax=uxmax(ng)+0.5_r8
177 ymin=uymin(ng)
178 ymax=uymax(ng)
179 xoff=0.0_r8
180 ELSE
181 xmin=rxmin(ng)
182 xmax=rxmax(ng)
183 ymin=rymin(ng)
184 ymax=rymax(ng)
185 xoff=0.5_r8
186 END IF
187 DO np=1,npos
188 xgrd=xpos(np)+xoff
189 ygrd=ypos(np)
190 bounded(np)=0.0_r8
191 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
192 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
193 i1=int(xgrd)
194 j1=int(ygrd)
195 i2=i1+1
196 j2=j1+1
197 IF (i2.gt.lm(ng)+1) THEN
198 i2=i1 ! station at the eastern boundary
199 END IF
200 IF (j2.gt.mm(ng)+1) THEN
201 j2=j1 ! station at the northern boundary
202 END IF
203 bounded(np)=1.0_r8
204 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
205 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
206 p1=1.0_r8-p2
207 q1=1.0_r8-q2
208 w111=p1*q1
209 w211=p2*q1
210 w121=p1*q2
211 w221=p2*q2
212#ifdef MASKING
213 w111=w111*grid(ng)%umask(i1,j1)
214 w211=w211*grid(ng)%umask(i2,j1)
215 w121=w121*grid(ng)%umask(i1,j2)
216 w221=w221*grid(ng)%umask(i2,j2)
217 wsum=w111+w211+w121+w221
218 IF (wsum.gt.0.0_r8) THEN
219 wsum=1.0_r8/wsum
220 w111=w111*wsum
221 w211=w211*wsum
222 w121=w121*wsum
223 w221=w221*wsum
224 ELSE
225 bounded(np)=0.0_r8
226 END IF
227#endif
228 apos(np)=ascl*(w111*a(i1,j1)+ &
229 & w211*a(i2,j1)+ &
230 & w121*a(i1,j2)+ &
231 & w221*a(i2,j2))
232 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
233 ELSE ! zero
234 apos(np)=aspv
235 END IF
236 END DO
237!
238!-----------------------------------------------------------------------
239! Interpolate from 2D field at V-points.
240!-----------------------------------------------------------------------
241!
242 ELSE IF (gtype.eq.v2dvar) THEN
243 IF (cgrid) THEN
244 xmin=vxmin(ng)
245 xmax=vxmax(ng)
246 ymin=vymin(ng)+0.5_r8
247 ymax=vymax(ng)+0.5_r8
248 yoff=0.0_r8
249 ELSE
250 xmin=rxmin(ng)
251 xmax=rxmax(ng)
252 ymin=rymin(ng)
253 ymax=rymax(ng)
254 yoff=0.5_r8
255 END IF
256 DO np=1,npos
257 xgrd=xpos(np)
258 ygrd=ypos(np)+yoff
259 bounded(np)=0.0_r8
260 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
261 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
262 i1=int(xgrd)
263 j1=int(ygrd)
264 i2=i1+1
265 j2=j1+1
266 IF (i2.gt.lm(ng)+1) THEN
267 i2=i1 ! station at the eastern boundary
268 END IF
269 IF (j2.gt.mm(ng)+1) THEN
270 j2=j1 ! station at the northern boundary
271 END IF
272 bounded(np)=1.0_r8
273 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
274 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
275 p1=1.0_r8-p2
276 q1=1.0_r8-q2
277 w111=p1*q1
278 w211=p2*q1
279 w121=p1*q2
280 w221=p2*q2
281#ifdef MASKING
282 w111=w111*grid(ng)%vmask(i1,j1)
283 w211=w211*grid(ng)%vmask(i2,j1)
284 w121=w121*grid(ng)%vmask(i1,j2)
285 w221=w221*grid(ng)%vmask(i2,j2)
286 wsum=w111+w211+w121+w221
287 IF (wsum.gt.0.0_r8) THEN
288 wsum=1.0_r8/wsum
289 w111=w111*wsum
290 w211=w211*wsum
291 w121=w121*wsum
292 w221=w221*wsum
293 ELSE
294 bounded(np)=0.0_r8
295 END IF
296#endif
297 apos(np)=ascl*(w111*a(i1,j1)+ &
298 & w211*a(i2,j1)+ &
299 & w121*a(i1,j2)+ &
300 & w221*a(i2,j2))
301 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
302 ELSE ! zero
303 apos(np)=aspv
304 END IF
305 END DO
306 END IF
307#ifdef DISTRIBUTE
308!
309!-----------------------------------------------------------------------
310! Collect all extracted data.
311!-----------------------------------------------------------------------
312!
313 CALL mp_collect (ng, model, npos, aspv, apos)
314 CALL mp_collect (ng, model, npos, 0.0_r8, bounded)
315#endif
316!
317!-----------------------------------------------------------------------
318! Set unbounded data to special value.
319!-----------------------------------------------------------------------
320!
321 DO np=1,npos
322 IF (bounded(np).lt.1.0_r8) THEN
323 apos(np)=spval
324 END IF
325 END DO
326 RETURN
327 END SUBROUTINE extract_sta2d
328
329#ifdef SOLVE3D
330!
331!***********************************************************************
332 SUBROUTINE extract_sta3d (ng, model, Cgrid, ifield, gtype, &
333 & LBi, UBi, LBj, UBj, LBk, UBk, Ascl, A, &
334 & Npos, Xpos, Ypos, Zpos, Apos)
335!***********************************************************************
336!
337 USE mod_param
338 USE mod_parallel
339 USE mod_grid
340 USE mod_ncparam
341 USE mod_scalars
342
343# ifdef DISTRIBUTE
344!
345 USE distribute_mod, ONLY : mp_collect
346# endif
347!
348! Imported variable declarations.
349!
350 logical, intent(in) :: Cgrid
351!
352 integer, intent(in) :: ng, model, ifield, gtype, Npos
353 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
354 real(dp), intent(in) :: Ascl
355!
356# ifdef ASSUMED_SHAPE
357 real(r8), intent(in) :: A(LBi:,LBj:,LBk:)
358 real(r8), intent(in) :: Xpos(:), Ypos(:), Zpos(:)
359 real(r8), intent(out) :: Apos(:)
360# else
361 real(r8), intent(in) :: A(LBi:UBi,LBj:UBj,LBk:UBk)
362 real(r8), intent(in) :: Xpos(Npos), Ypos(Npos), Zpos(Npos)
363 real(r8), intent(out) :: Apos(Npos)
364# endif
365!
366! Local variable declarations.
367!
368 integer :: i1, i2, j1, j2, k, k1, k2, np
369
370 real(r8), parameter :: Aspv = 0.0_r8
371
372 real(r8) :: Xmin, Xmax, Ymin, Ymax
373 real(r8) :: Xgrd, Xoff, Ygrd, Yoff, Zgrd, Zbot, Ztop
374 real(r8) :: dz, p1, p2, q1, q2, r1, r2, wsum
375 real(r8) :: w111, w211, w121, w221, w112, w212, w122, w222
376
377 real(r8), dimension(Npos) :: bounded
378!
379!-----------------------------------------------------------------------
380! Interpolate from 3D field at RHO-points.
381!-----------------------------------------------------------------------
382!
383 IF (gtype.eq.r3dvar) THEN
384 xmin=rxmin(ng)
385 xmax=rxmax(ng)
386 ymin=rymin(ng)
387 ymax=rymax(ng)
388 DO np=1,npos
389 xgrd=xpos(np)
390 ygrd=ypos(np)
391 zgrd=zpos(np)
392 bounded(np)=0.0_r8
393 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
394 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
395 i1=int(xgrd)
396 j1=int(ygrd)
397 i2=i1+1
398 j2=j1+1
399 IF (i2.gt.lm(ng)+1) THEN
400 i2=i1 ! station at the eastern boundary
401 END IF
402 IF (j2.gt.mm(ng)+1) THEN
403 j2=j1 ! station at the northern boundary
404 END IF
405 bounded(np)=1.0_r8
406 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
407 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
408 p1=1.0_r8-p2
409 q1=1.0_r8-q2
410 w111=p1*q1
411 w211=p2*q1
412 w121=p1*q2
413 w221=p2*q2
414 w112=0.0_r8
415 w212=0.0_r8
416 w122=0.0_r8
417 w222=0.0_r8
418 IF (zgrd.ge.0.0_r8) THEN
419 k1=int(zgrd)
420 k2=int(zgrd)
421 r1=1.0_r8
422 r2=0.0_r8
423 ELSE
424 ztop=grid(ng)%z_r(i1,j1,n(ng))
425 zbot=grid(ng)%z_r(i1,j1,1)
426 IF (zgrd.ge.ztop) THEN
427 k1=n(ng) ! If shallower, assign
428 k2=n(ng) ! station to surface
429 r1=1.0_r8 ! level
430 r2=0.0_r8
431 ELSE IF (zbot.ge.zgrd) THEN
432 k1=1 ! If deeper, assign
433 k2=1 ! station to bottom
434 r1=1.0_r8 ! level
435 r2=0.0_r8
436 ELSE
437 DO k=n(ng),2,-1
438 ztop=grid(ng)%z_r(i1,j1,k)
439 zbot=grid(ng)%z_r(i1,j1,k-1)
440 IF ((ztop.gt.zgrd).and.(zgrd.ge.zbot)) THEN
441 k1=k-1
442 k2=k
443 END IF
444 END DO
445 dz=grid(ng)%z_r(i1,j1,k2)-grid(ng)%z_r(i1,j1,k1)
446 r2=(zgrd-grid(ng)%z_r(i1,j1,k1))/dz
447 r1=1.0_r8-r2
448 END IF
449 END IF
450 w112=w111*r2
451 w212=w211*r2
452 w122=w121*r2
453 w222=w221*r2
454 w111=w111*r1
455 w211=w211*r1
456 w121=w121*r1
457 w221=w221*r1
458# ifdef MASKING
459 w111=w111*grid(ng)%rmask(i1,j1)
460 w211=w211*grid(ng)%rmask(i2,j1)
461 w121=w121*grid(ng)%rmask(i1,j2)
462 w221=w221*grid(ng)%rmask(i2,j2)
463 w112=w112*grid(ng)%rmask(i1,j1)
464 w212=w212*grid(ng)%rmask(i2,j1)
465 w122=w122*grid(ng)%rmask(i1,j2)
466 w222=w222*grid(ng)%rmask(i2,j2)
467 wsum=w111+w211+w121+w221+w112+w212+w122+w222
468 IF (wsum.gt.0.0_r8) THEN
469 wsum=1.0_r8/wsum
470 w111=w111*wsum
471 w211=w211*wsum
472 w121=w121*wsum
473 w221=w221*wsum
474 w112=w112*wsum
475 w212=w212*wsum
476 w122=w122*wsum
477 w222=w222*wsum
478 ELSE
479 bounded(np)=0.0_r8
480 END IF
481# endif
482 apos(np)=ascl*(w111*a(i1,j1,k1)+ &
483 & w211*a(i2,j1,k1)+ &
484 & w121*a(i1,j2,k1)+ &
485 & w221*a(i2,j2,k1)+ &
486 & w112*a(i1,j1,k2)+ &
487 & w212*a(i2,j1,k2)+ &
488 & w122*a(i1,j2,k2)+ &
489 & w222*a(i2,j2,k2))
490 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
491 ELSE ! zero
492 apos(np)=aspv
493 END IF
494 END DO
495!
496!-----------------------------------------------------------------------
497! Interpolate from 3D field at U-points.
498!-----------------------------------------------------------------------
499!
500 ELSE IF (gtype.eq.u3dvar) THEN
501 IF (cgrid) THEN
502 xmin=uxmin(ng)+0.5_r8
503 xmax=uxmax(ng)+0.5_r8
504 ymin=uymin(ng)
505 ymax=uymax(ng)
506 xoff=0.0_r8
507 ELSE
508 xmin=rxmin(ng)
509 xmax=rxmax(ng)
510 ymin=rymin(ng)
511 ymax=rymax(ng)
512 xoff=0.5_r8
513 END IF
514 DO np=1,npos
515 xgrd=xpos(np)+xoff
516 ygrd=ypos(np)
517 zgrd=zpos(np)
518 bounded(np)=0.0_r8
519 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
520 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
521 i1=int(xgrd)
522 j1=int(ygrd)
523 i2=i1+1
524 j2=j1+1
525 IF (i2.gt.lm(ng)+1) THEN
526 i2=i1 ! station at the eastern boundary
527 END IF
528 IF (j2.gt.mm(ng)+1) THEN
529 j2=j1 ! station at the northern boundary
530 END IF
531 bounded(np)=1.0_r8
532 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
533 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
534 p1=1.0_r8-p2
535 q1=1.0_r8-q2
536 w111=p1*q1
537 w211=p2*q1
538 w121=p1*q2
539 w221=p2*q2
540 w112=0.0_r8
541 w212=0.0_r8
542 w122=0.0_r8
543 w222=0.0_r8
544 IF (zgrd.ge.0.0_r8) THEN
545 k1=int(zgrd)
546 k2=int(zgrd)
547 r1=1.0_r8
548 r2=0.0_r8
549 ELSE
550 ztop=0.5_r8*(grid(ng)%z_r(i1-1,j1,n(ng))+ &
551 & grid(ng)%z_r(i1 ,j1,n(ng)))
552 zbot=0.5_r8*(grid(ng)%z_r(i1-1,j1,1)+ &
553 & grid(ng)%z_r(i1 ,j1,1))
554 IF (zgrd.ge.ztop) THEN
555 k1=n(ng) ! If shallower, assign
556 k2=n(ng) ! station to surface
557 r1=1.0_r8 ! level
558 r2=0.0_r8
559 ELSE IF (zbot.ge.zgrd) THEN
560 k1=1 ! If deeper, assign
561 k2=1 ! station to bottom
562 r1=1.0_r8 ! level
563 r2=0.0_r8
564 ELSE
565 DO k=n(ng),2,-1
566 ztop=0.5_r8*(grid(ng)%z_r(i1-1,j1,k)+ &
567 & grid(ng)%z_r(i1 ,j1,k))
568 zbot=0.5_r8*(grid(ng)%z_r(i1-1,j1,k-1)+ &
569 & grid(ng)%z_r(i1 ,j1,k-1))
570 IF ((ztop.gt.zgrd).and.(zgrd.ge.zbot)) THEN
571 k1=k-1
572 k2=k
573 END IF
574 END DO
575 dz=0.5_r8*((grid(ng)%z_r(i1-1,j1,k2)+ &
576 & grid(ng)%z_r(i1 ,j1,k2))- &
577 & (grid(ng)%z_r(i1-1,j1,k1)+ &
578 & grid(ng)%z_r(i1 ,j1,k1)))
579 r2=(zgrd-0.5_r8*(grid(ng)%z_r(i1-1,j1,k1)+ &
580 & grid(ng)%z_r(i1 ,j1,k1)))/dz
581 r1=1.0_r8-r2
582 END IF
583 END IF
584 w112=w111*r2
585 w212=w211*r2
586 w122=w121*r2
587 w222=w221*r2
588 w111=w111*r1
589 w211=w211*r1
590 w121=w121*r1
591 w221=w221*r1
592# ifdef MASKING
593 w111=w111*grid(ng)%umask(i1,j1)
594 w211=w211*grid(ng)%umask(i2,j1)
595 w121=w121*grid(ng)%umask(i1,j2)
596 w221=w221*grid(ng)%umask(i2,j2)
597 w112=w112*grid(ng)%umask(i1,j1)
598 w212=w212*grid(ng)%umask(i2,j1)
599 w122=w122*grid(ng)%umask(i1,j2)
600 w222=w222*grid(ng)%umask(i2,j2)
601 wsum=w111+w211+w121+w221+w112+w212+w122+w222
602 IF (wsum.gt.0.0_r8) THEN
603 wsum=1.0_r8/wsum
604 w111=w111*wsum
605 w211=w211*wsum
606 w121=w121*wsum
607 w221=w221*wsum
608 w112=w112*wsum
609 w212=w212*wsum
610 w122=w122*wsum
611 w222=w222*wsum
612 ELSE
613 bounded(np)=0.0_r8
614 END IF
615# endif
616 apos(np)=ascl*(w111*a(i1,j1,k1)+ &
617 & w211*a(i2,j1,k1)+ &
618 & w121*a(i1,j2,k1)+ &
619 & w221*a(i2,j2,k1)+ &
620 & w112*a(i1,j1,k2)+ &
621 & w212*a(i2,j1,k2)+ &
622 & w122*a(i1,j2,k2)+ &
623 & w222*a(i2,j2,k2))
624 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
625 ELSE ! zero
626 apos(np)=aspv
627 END IF
628 END DO
629!
630!-----------------------------------------------------------------------
631! Interpolate from 3D field at V-points.
632!-----------------------------------------------------------------------
633!
634 ELSE IF (gtype.eq.v3dvar) THEN
635 IF (cgrid) THEN
636 xmin=vxmin(ng)
637 xmax=vxmax(ng)
638 ymin=vymin(ng)+0.5_r8
639 ymax=vymax(ng)+0.5_r8
640 yoff=0.0_r8
641 ELSE
642 xmin=rxmin(ng)
643 xmax=rxmax(ng)
644 ymin=rymin(ng)
645 ymax=rymax(ng)
646 yoff=0.5_r8
647 END IF
648 DO np=1,npos
649 xgrd=xpos(np)
650 ygrd=ypos(np)+yoff
651 zgrd=zpos(np)
652 bounded(np)=0.0_r8
653 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
654 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
655 i1=int(xgrd)
656 j1=int(ygrd)
657 i2=i1+1
658 j2=j1+1
659 IF (i2.gt.lm(ng)+1) THEN
660 i2=i1 ! station at the eastern boundary
661 END IF
662 IF (j2.gt.mm(ng)+1) THEN
663 j2=j1 ! station at the northern boundary
664 END IF
665 bounded(np)=1.0_r8
666 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
667 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
668 p1=1.0_r8-p2
669 q1=1.0_r8-q2
670 w111=p1*q1
671 w211=p2*q1
672 w121=p1*q2
673 w221=p2*q2
674 w112=0.0_r8
675 w212=0.0_r8
676 w122=0.0_r8
677 w222=0.0_r8
678 IF (zgrd.ge.0.0_r8) THEN
679 k1=int(zgrd)
680 k2=int(zgrd)
681 r1=1.0_r8
682 r2=0.0_r8
683 ELSE
684 ztop=0.5_r8*(grid(ng)%z_r(i1,j1-1,n(ng))+ &
685 & grid(ng)%z_r(i1,j1, n(ng)))
686 zbot=0.5_r8*(grid(ng)%z_r(i1,j1-1,1)+ &
687 & grid(ng)%z_r(i1,j1 ,1))
688 IF (zgrd.ge.ztop) THEN
689 k1=n(ng) ! If shallower, assign
690 k2=n(ng) ! station to surface
691 r1=1.0_r8 ! level
692 r2=0.0_r8
693 ELSE IF (zbot.ge.zgrd) THEN
694 k1=1 ! If deeper, assign
695 k2=1 ! station to bottom
696 r1=1.0_r8 ! level
697 r2=0.0_r8
698 ELSE
699 DO k=n(ng),2,-1
700 ztop=0.5_r8*(grid(ng)%z_r(i1,j1-1,k)+ &
701 & grid(ng)%z_r(i1,j1 ,k))
702 zbot=0.5_r8*(grid(ng)%z_r(i1,j1-1,k-1)+ &
703 & grid(ng)%z_r(i1,j1 ,k-1))
704 IF ((ztop.gt.zgrd).and.(zgrd.ge.zbot)) THEN
705 k1=k-1
706 k2=k
707 END IF
708 END DO
709 dz=0.5_r8*((grid(ng)%z_r(i1,j1-1,k2)+ &
710 & grid(ng)%z_r(i1,j1 ,k2))- &
711 & (grid(ng)%z_r(i1,j1-1,k1)+ &
712 & grid(ng)%z_r(i1,j1 ,k1)))
713 r2=(zgrd-0.5_r8*(grid(ng)%z_r(i1,j1-1,k1)+ &
714 & grid(ng)%z_r(i1,j1 ,k1)))/dz
715 r1=1.0_r8-r2
716 END IF
717 END IF
718 w112=w111*r2
719 w212=w211*r2
720 w122=w121*r2
721 w222=w221*r2
722 w111=w111*r1
723 w211=w211*r1
724 w121=w121*r1
725 w221=w221*r1
726# ifdef MASKING
727 w111=w111*grid(ng)%vmask(i1,j1)
728 w211=w211*grid(ng)%vmask(i2,j1)
729 w121=w121*grid(ng)%vmask(i1,j2)
730 w221=w221*grid(ng)%vmask(i2,j2)
731 w112=w112*grid(ng)%vmask(i1,j1)
732 w212=w212*grid(ng)%vmask(i2,j1)
733 w122=w122*grid(ng)%vmask(i1,j2)
734 w222=w222*grid(ng)%vmask(i2,j2)
735 wsum=w111+w211+w121+w221+w112+w212+w122+w222
736 IF (wsum.gt.0.0_r8) THEN
737 wsum=1.0_r8/wsum
738 w111=w111*wsum
739 w211=w211*wsum
740 w121=w121*wsum
741 w221=w221*wsum
742 w112=w112*wsum
743 w212=w212*wsum
744 w122=w122*wsum
745 w222=w222*wsum
746 ELSE
747 bounded(np)=0.0_r8
748 END IF
749# endif
750 apos(np)=ascl*(w111*a(i1,j1,k1)+ &
751 & w211*a(i2,j1,k1)+ &
752 & w121*a(i1,j2,k1)+ &
753 & w221*a(i2,j2,k1)+ &
754 & w112*a(i1,j1,k2)+ &
755 & w212*a(i2,j1,k2)+ &
756 & w122*a(i1,j2,k2)+ &
757 & w222*a(i2,j2,k2))
758 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
759 ELSE ! zero
760 apos(np)=aspv
761 END IF
762 END DO
763!
764!-----------------------------------------------------------------------
765! Interpolate from 3D field at W-points.
766!-----------------------------------------------------------------------
767!
768 ELSE IF (gtype.eq.w3dvar) THEN
769 xmin=rxmin(ng)
770 xmax=rxmax(ng)
771 ymin=rymin(ng)
772 ymax=rymax(ng)
773 DO np=1,npos
774 xgrd=xpos(np)
775 ygrd=ypos(np)
776 zgrd=zpos(np)
777 bounded(np)=0.0_r8
778 IF (((xmin.le.xgrd).and.(xgrd.lt.xmax)).and. &
779 & ((ymin.le.ygrd).and.(ygrd.lt.ymax))) THEN
780 i1=int(xgrd)
781 j1=int(ygrd)
782 i2=i1+1
783 j2=j1+1
784 IF (i2.gt.lm(ng)+1) THEN
785 i2=i1 ! station at the eastern boundary
786 END IF
787 IF (j2.gt.mm(ng)+1) THEN
788 j2=j1 ! station at the northern boundary
789 END IF
790 bounded(np)=1.0_r8
791 p2=real(i2-i1,r8)*(xgrd-real(i1,r8))
792 q2=real(j2-j1,r8)*(ygrd-real(j1,r8))
793 p1=1.0_r8-p2
794 q1=1.0_r8-q2
795 w111=p1*q1
796 w211=p2*q1
797 w121=p1*q2
798 w221=p2*q2
799 w112=0.0_r8
800 w212=0.0_r8
801 w122=0.0_r8
802 w222=0.0_r8
803 IF (zgrd.ge.0.0_r8) THEN
804 k1=int(zgrd)
805 k2=int(zgrd)
806 r1=1.0_r8
807 r2=0.0_r8
808 ELSE
809 ztop=grid(ng)%z_w(i1,j1,n(ng))
810 zbot=grid(ng)%z_w(i1,j1,0)
811 IF (zgrd.ge.ztop) THEN
812 k1=n(ng) ! If shallower, assign
813 k2=n(ng) ! station to surface
814 r1=1.0_r8 ! level
815 r2=0.0_r8
816 ELSE IF (zbot.ge.zgrd) THEN
817 k1=0 ! If deeper, assign
818 k2=0 ! station to bottom
819 r1=1.0_r8 ! level
820 r2=0.0_r8
821 ELSE
822 DO k=n(ng),2,-1
823 ztop=grid(ng)%z_w(i1,j1,k)
824 zbot=grid(ng)%z_w(i1,j1,k-1)
825 IF ((ztop.gt.zgrd).and.(zgrd.ge.zbot)) THEN
826 k1=k-1
827 k2=k
828 END IF
829 END DO
830 dz=grid(ng)%z_w(i1,j1,k2)-grid(ng)%z_w(i1,j1,k1)
831 r2=(zgrd-grid(ng)%z_w(i1,j1,k1))/dz
832 r1=1.0_r8-r2
833 END IF
834 END IF
835 w112=w111*r2
836 w212=w211*r2
837 w122=w121*r2
838 w222=w221*r2
839 w111=w111*r1
840 w211=w211*r1
841 w121=w121*r1
842 w221=w221*r1
843# ifdef MASKING
844 w111=w111*grid(ng)%rmask(i1,j1)
845 w211=w211*grid(ng)%rmask(i2,j1)
846 w121=w121*grid(ng)%rmask(i1,j2)
847 w221=w221*grid(ng)%rmask(i2,j2)
848 w112=w112*grid(ng)%rmask(i1,j1)
849 w212=w212*grid(ng)%rmask(i2,j1)
850 w122=w122*grid(ng)%rmask(i1,j2)
851 w222=w222*grid(ng)%rmask(i2,j2)
852 wsum=w111+w211+w121+w221+w112+w212+w122+w222
853 IF (wsum.gt.0.0_r8) THEN
854 wsum=1.0_r8/wsum
855 w111=w111*wsum
856 w211=w211*wsum
857 w121=w121*wsum
858 w221=w221*wsum
859 w112=w112*wsum
860 w212=w212*wsum
861 w122=w122*wsum
862 w222=w222*wsum
863 ELSE
864 bounded(np)=0.0_r8
865 END IF
866# endif
867 apos(np)=ascl*(w111*a(i1,j1,k1)+ &
868 & w211*a(i2,j1,k1)+ &
869 & w121*a(i1,j2,k1)+ &
870 & w221*a(i2,j2,k1)+ &
871 & w112*a(i1,j1,k2)+ &
872 & w212*a(i2,j1,k2)+ &
873 & w122*a(i1,j2,k2)+ &
874 & w222*a(i2,j2,k2))
875 IF (abs(apos(np)).eq.0.0_r8) apos(np)=0.0_r8 ! positive
876 ELSE ! zero
877 apos(np)=aspv
878 END IF
879 END DO
880 END IF
881# ifdef DISTRIBUTE
882!
883!-----------------------------------------------------------------------
884! Collect all extracted data.
885!-----------------------------------------------------------------------
886!
887 CALL mp_collect (ng, model, npos, aspv, apos)
888 CALL mp_collect (ng, model, npos, 0.0_r8, bounded)
889# endif
890!
891!-----------------------------------------------------------------------
892! Set unbounded data to special value.
893!-----------------------------------------------------------------------
894!
895 DO np=1,npos
896 IF (bounded(np).lt.1.0_r8) THEN
897 apos(np)=spval
898 END IF
899 END DO
900 RETURN
901 END SUBROUTINE extract_sta3d
902#endif
903 END MODULE extract_sta_mod
904
distribute_mod::mp_collect
Definition distribute.F:108
distribute_mod
Definition distribute.F:3
extract_sta_mod
Definition extract_sta.F:2
extract_sta_mod::extract_sta2d
subroutine extract_sta2d(ng, model, cgrid, ifield, gtype, lbi, ubi, lbj, ubj, ascl, a, npos, xpos, ypos, apos)
Definition extract_sta.F:65
extract_sta_mod::extract_sta3d
subroutine extract_sta3d(ng, model, cgrid, ifield, gtype, lbi, ubi, lbj, ubj, lbk, ubk, ascl, a, npos, xpos, ypos, zpos, apos)
Definition extract_sta.F:335
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::rymin
real(r8), dimension(:), allocatable rymin
Definition mod_ncparam.F:135
mod_ncparam::vymin
real(r8), dimension(:), allocatable vymin
Definition mod_ncparam.F:161
mod_ncparam::rymax
real(r8), dimension(:), allocatable rymax
Definition mod_ncparam.F:136
mod_ncparam::uymin
real(r8), dimension(:), allocatable uymin
Definition mod_ncparam.F:148
mod_ncparam::vymax
real(r8), dimension(:), allocatable vymax
Definition mod_ncparam.F:162
mod_ncparam::uxmin
real(r8), dimension(:), allocatable uxmin
Definition mod_ncparam.F:146
mod_ncparam::uxmax
real(r8), dimension(:), allocatable uxmax
Definition mod_ncparam.F:147
mod_ncparam::rxmax
real(r8), dimension(:), allocatable rxmax
Definition mod_ncparam.F:134
mod_ncparam::uymax
real(r8), dimension(:), allocatable uymax
Definition mod_ncparam.F:149
mod_ncparam::vxmin
real(r8), dimension(:), allocatable vxmin
Definition mod_ncparam.F:159
mod_ncparam::vxmax
real(r8), dimension(:), allocatable vxmax
Definition mod_ncparam.F:160
mod_ncparam::rxmin
real(r8), dimension(:), allocatable rxmin
Definition mod_ncparam.F:133
mod_parallel
Definition mod_parallel.F:2
mod_param
Definition mod_param.F:2
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::r3dvar
integer, parameter r3dvar
Definition mod_param.F:721
mod_param::u3dvar
integer, parameter u3dvar
Definition mod_param.F:722
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::u2dvar
integer, parameter u2dvar
Definition mod_param.F:718
mod_param::w3dvar
integer, parameter w3dvar
Definition mod_param.F:724
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_param::r2dvar
integer, parameter r2dvar
Definition mod_param.F:717
mod_param::v2dvar
integer, parameter v2dvar
Definition mod_param.F:719
mod_param::v3dvar
integer, parameter v3dvar
Definition mod_param.F:723
mod_scalars
Definition mod_scalars.F:2
mod_scalars::spval
real(dp), parameter spval
Definition mod_scalars.F:821