335
336
342
343# ifdef DISTRIBUTE
344
346# endif
347
348
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
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
381
382
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
401 END IF
402 IF (j2.gt.
mm(ng)+1)
THEN
403 j2=j1
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
429 r1=1.0_r8
430 r2=0.0_r8
431 ELSE IF (zbot.ge.zgrd) THEN
432 k1=1
433 k2=1
434 r1=1.0_r8
435 r2=0.0_r8
436 ELSE
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
491 ELSE
492 apos(np)=aspv
493 END IF
494 END DO
495
496
497
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
506 xoff=0.0_r8
507 ELSE
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
527 END IF
528 IF (j2.gt.
mm(ng)+1)
THEN
529 j2=j1
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
557 r1=1.0_r8
558 r2=0.0_r8
559 ELSE IF (zbot.ge.zgrd) THEN
560 k1=1
561 k2=1
562 r1=1.0_r8
563 r2=0.0_r8
564 ELSE
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
625 ELSE
626 apos(np)=aspv
627 END IF
628 END DO
629
630
631
632
633
634 ELSE IF (gtype.eq.
v3dvar)
THEN
635 IF (cgrid) THEN
638 ymin=
vymin(ng)+0.5_r8
639 ymax=
vymax(ng)+0.5_r8
640 yoff=0.0_r8
641 ELSE
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
661 END IF
662 IF (j2.gt.
mm(ng)+1)
THEN
663 j2=j1
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
691 r1=1.0_r8
692 r2=0.0_r8
693 ELSE IF (zbot.ge.zgrd) THEN
694 k1=1
695 k2=1
696 r1=1.0_r8
697 r2=0.0_r8
698 ELSE
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
759 ELSE
760 apos(np)=aspv
761 END IF
762 END DO
763
764
765
766
767
768 ELSE IF (gtype.eq.
w3dvar)
THEN
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
786 END IF
787 IF (j2.gt.
mm(ng)+1)
THEN
788 j2=j1
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
814 r1=1.0_r8
815 r2=0.0_r8
816 ELSE IF (zbot.ge.zgrd) THEN
817 k1=0
818 k2=0
819 r1=1.0_r8
820 r2=0.0_r8
821 ELSE
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
876 ELSE
877 apos(np)=aspv
878 END IF
879 END DO
880 END IF
881# ifdef DISTRIBUTE
882
883
884
885
886
888 CALL mp_collect (ng, model, npos, 0.0_r8, bounded)
889# endif
890
891
892
893
894
895 DO np=1,npos
896 IF (bounded(np).lt.1.0_r8) THEN
898 END IF
899 END DO
900 RETURN
integer, dimension(:), allocatable n
integer, parameter r3dvar
integer, parameter u3dvar
integer, parameter w3dvar
integer, parameter v3dvar