621
624
625
626
627#include "implicit_f.inc"
628
629
630
631 TYPE(TTABLE) TABLE
632 INTEGER ,INTENT(IN) :: NEL
633 INTEGER ,VALUE ,INTENT(IN) :: DIMX
634 INTEGER ,DIMENSION(DIMX,TABLE%NDIM) :: IPOS
635 my_real ,
DIMENSION(DIMX,TABLE%NDIM) :: xx
636 my_real ,
DIMENSION(NEL) :: yy, dydx1
637
638
639
640 LOGICAL, DIMENSION(NEL) :: NEED_TO_COMPUTE
641 INTEGER NDIM, K, (4), I, IB(NEL,2,2,2,2),
642 . IP,IN,IM,IL,P,N,M,L,N1,N12,N123
643 my_real :: dx2,r(nel,4),unr(nel,4),dx2_0(nel)
644 INTEGER :: J
645 INTEGER :: NINDX_1,M_INDX1,NINDX_2,M_INDX2
646 INTEGER, DIMENSION(NEL) :: INDX_1,INDX_2
647
648 r(1:nel,1:4) = zero
649 ndim=table%NDIM
650 IF( SIZE(xx,2) < table%NDIM )THEN
651 CALL ancmsg(msgid=36,anmode=aninfo,
652 . c1='TABLE INTERPOLATION')
654 END IF
655
656 DO k=1,ndim
657 nxk(k)=SIZE(table%X(k)%VALUES)
658 ENDDO
659
660 DO k=1,ndim
661 ipos(1:nel,k)=
max(ipos(1:nel,k),1)
662 nindx_1 = 0
663 m_indx1 = 0
664 nindx_2 = 0
665 m_indx2 = nxk(k) + 1
666#include "vectorize.inc"
667 DO i=1,nel
668 m = ipos(i,k)
669 dx2_0(i) = table%X(k)%VALUES(m) - xx(i,k)
670 IF(dx2_0(i) >= zero)THEN
671 nindx_1 = nindx_1 + 1
672 indx_1(nindx_1) = i
673 m_indx1 =
max(m_indx1,m)
674 ELSE
675 nindx_2 = nindx_2 + 1
676 indx_2(nindx_2) = i
677 m_indx2 =
min(m_indx2,m)
678 ENDIF
679 ENDDO
680
681 need_to_compute(1:nindx_1) = .true.
682 DO n=m_indx1,1,-1
683#include "vectorize.inc"
684 DO j=1,nindx_1
685 IF(need_to_compute(j)) THEN
686 i = indx_1(j)
687 m = ipos(i,k)
688 dx2 = table%X(k)%VALUES(n) - xx(i,k)
689 IF(dx2<zero.OR.n <=1)THEN
691 need_to_compute(j) = .false.
692 ENDIF
693 ENDIF
694 ENDDO
695 ENDDO
696 need_to_compute(1:nindx_2) = .true.
697 DO n=m_indx2,nxk(k)
698#include "vectorize.inc"
699 DO j=1,nindx_2
700 IF(need_to_compute(j)) THEN
701 i = indx_2(j)
702 m = ipos(i,k)
703 dx2 = table%X(k)%VALUES(n) - xx(i,k)
704 IF(dx2>=zero.OR.n==nxk(k))THEN
705 ipos(i,k)=n-1
706 need_to_compute(j) = .false.
707 ENDIF
708 ENDIF
709 ENDDO
710 ENDDO
711 ENDDO
712
713 DO k=1,ndim
714#include "vectorize.inc"
715 DO i=1,nel
716 n = ipos(i,k)
717 r(i,k) =(table%X(k)%VALUES(n+1)-xx(i,k))/
718 . (table%X(k)%VALUES(n+1)-table%X(k)%VALUES(n))
719 END DO
720 END DO
721
722 SELECT CASE(ndim)
723
724 CASE(4)
725
726 n1 =nxk(1)
727 n12 =nxk(1)*nxk(2)
728 n123=n12 *nxk(3)
729 DO i=1,nel
730 DO p=0,1
731 ip=n123*(ipos(i,4)-1+p)
732 DO n=0,1
733 in=n12*(ipos(i,3)-1+n)
734 DO m=0,1
735 im=n1*(ipos(i,2)-1+m)
736 DO l=0,1
737 il=ipos(i,1)+l
738 ib(i,l+1,m+1,n+1,p+1)=ip+in+im+il
739 END DO
740 END DO
741 END DO
742 END DO
743 END DO
744
745 unr(1:nel,1:4)=one-r(1:nel,1:4)
746#include "vectorize.inc"
747 DO i=1,nel
748
749 yy(i)=
750 . r(i,4)*(r(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,1,1))
751 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,1,1)))
752 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,1,1))
753 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,1,1))))
754 . +unr(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,2,1))
755 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,2,1)))
756 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,2,1))
757 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,2,1)))))
758 . +unr(i,4)*(r(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,1,1))
759 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,1,1)))
760 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,1,1))
761 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,1,1))))
762 . +unr(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,2,1))
763 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,2,1)))
764 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,2,1))
765 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,2,1)))))
766
767 dydx1(i)=
768 . (r(i,4)*(r(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,1,1))
769 . -table%Y%VALUES(ib(i,1,1,1,1)))
770 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,1,1))
771 . -table%Y%VALUES(ib(i,1,2,1,1))))
772 . +unr(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,2,1))
773 . -table%Y%VALUES(ib(i,1,1,2,1)))
774 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,2,1))
775 . -table%Y%VALUES(ib(i,1,2,2,1)))))
776 . +unr(i,4)*(r(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,1,1))
777 . -table%Y%VALUES(ib(i,1,1,1,1)))
778 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,1,1))
779 . -table%Y%VALUES(ib(i,1,2,1,1))))
780 . +unr(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,2,1))
781 . -table%Y%VALUES(ib(i,1,1,2,1)))
782 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,2,1))
783 . -table%Y%VALUES(ib(i,1,2,2,1))))))/
784 . (table%X(1)%VALUES(ipos(i,1)+1)-table%X(1)%VALUES(ipos(i,1)))
785
786 END DO
787
788 CASE(3)
789
790 n1 =nxk(1)
791 n12 =nxk(1)*nxk(2)
792 DO i=1,nel
793 DO n=0,1
794 in=n12*(ipos(i,3)-1+n)
795 DO m=0,1
796 im=n1*(ipos(i,2)-1+m)
797 DO l=0,1
798 il=ipos(i,1)+l
799 ib(i,l+1,m+1,n+1,1)=in+im+il
800 END DO
801 END DO
802 END DO
803 END DO
804
805 unr(1:nel,1:3)=one-r(1:nel,1:3)
806#include "vectorize.inc"
807 DO i=1,nel
808
809 yy(i)=
810 . (r(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,1,1))
811 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,1,1)))
812 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,1,1))
813 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,1,1))))
814 . +unr(i,3)*(r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,2,1))
815 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,2,1)))
816 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,2,1))
817 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,2,1)))))
818
819 dydx1(i)=
820 . (r(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,1,1))
821 . -table%Y%VALUES(ib(i,1,1,1,1)))
822 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,1,1))
823 . -table%Y%VALUES(ib(i,1,2,1,1))))
824 . +unr(i,3)*(r(i,2)*( table%Y%VALUES(ib(i,2,1,2,1))
825 . -table%Y%VALUES(ib(i,1,1,2,1)))
826 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,2,1))
827 . -table%Y%VALUES(ib(i,1,2,2,1)))))/
828 . (table%X(1)%VALUES(ipos(i,1)+1)-table%X(1)%VALUES(ipos(i,1)))
829
830 END DO
831
832 CASE(2)
833
834 n1 =nxk(1)
835 DO i=1,nel
836 DO m=0,1
837 im=n1*(ipos(i,2)-1+m)
838 DO l=0,1
839 il=ipos(i,1)+l
840 ib(i,l+1,m+1,1,1)=im+il
841 END DO
842 END DO
843 END DO
844
845 unr(1:nel,1:2)=one-r(1:nel,1:2)
846#include "vectorize.inc"
847 DO i=1,nel
848
849 yy(i)=
850 . (r(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,1,1,1))
851 . +unr(i,1)*table%Y%VALUES(ib(i,2,1,1,1)))
852 . +unr(i,2)*(r(i,1)*table%Y%VALUES(ib(i,1,2,1,1))
853 . +unr(i,1)*table%Y%VALUES(ib(i,2,2,1,1))))
854 dydx1(i)=
855 . (r(i,2)*( table%Y%VALUES(ib(i,2,1,1,1))
856 . -table%Y%VALUES(ib(i,1,1,1,1)))
857 . +unr(i,2)*( table%Y%VALUES(ib(i,2,2,1,1))
858 . -table%Y%VALUES(ib(i,1,2,1,1))))/
859 . (table%X(1)%VALUES(ipos(i,1)+1)-table%X(1)%VALUES(ipos(i,1)))
860 END DO
861
862 CASE(1)
863
864 unr(1:nel,1:1)=one-r(1:nel,1:1)
865#include "vectorize.inc"
866 DO i=1,nel
867
868 yy(i)= r(i,1)*table%Y%VALUES(ipos(i,1))
869 . +unr(i,1)*table%Y%VALUES(ipos(i,1)+1)
870 dydx1(i)=(table%Y%VALUES(ipos(i,1)+1)-table%Y%VALUES(ipos(i,1)))/
871 . (table%X(1)%VALUES(ipos(i,1)+1)-table%X(1)%VALUES(ipos(i,1)))
872 END DO
873
874 END SELECT
875 RETURN