ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 2361 by gezelter, Wed Oct 12 21:00:50 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 518 | Line 518 | contains  
518      drdxi = -d(1) / rij
519      drdyi = -d(2) / rij
520      drdzi = -d(3) / rij
521 <    drduxi = 0.0d0
522 <    drduyi = 0.0d0
523 <    drduzi = 0.0d0
521 >    drduxi = 0.0_dp
522 >    drduyi = 0.0_dp
523 >    drduzi = 0.0_dp
524  
525      drdxj = d(1) / rij
526      drdyj = d(2) / rij
527      drdzj = d(3) / rij
528 <    drduxj = 0.0d0
529 <    drduyj = 0.0d0
530 <    drduzj = 0.0d0
528 >    drduxj = 0.0_dp
529 >    drduyj = 0.0_dp
530 >    drduzj = 0.0_dp
531  
532      ! find the atom type id (atid) for each atom:
533   #ifdef IS_MPI
# Line 549 | Line 549 | contains  
549         sigma_i = ShapeMap%Shapes(st1)%sigma
550         s_i = ShapeMap%Shapes(st1)%sigma
551         eps_i = ShapeMap%Shapes(st1)%epsilon
552 <       dsigmaidx = 0.0d0
553 <       dsigmaidy = 0.0d0
554 <       dsigmaidz = 0.0d0
555 <       dsigmaidux = 0.0d0
556 <       dsigmaiduy = 0.0d0
557 <       dsigmaiduz = 0.0d0
558 <       dsidx = 0.0d0
559 <       dsidy = 0.0d0
560 <       dsidz = 0.0d0
561 <       dsidux = 0.0d0
562 <       dsiduy = 0.0d0
563 <       dsiduz = 0.0d0
564 <       depsidx = 0.0d0
565 <       depsidy = 0.0d0
566 <       depsidz = 0.0d0
567 <       depsidux = 0.0d0
568 <       depsiduy = 0.0d0
569 <       depsiduz = 0.0d0
552 >       dsigmaidx = 0.0_dp
553 >       dsigmaidy = 0.0_dp
554 >       dsigmaidz = 0.0_dp
555 >       dsigmaidux = 0.0_dp
556 >       dsigmaiduy = 0.0_dp
557 >       dsigmaiduz = 0.0_dp
558 >       dsidx = 0.0_dp
559 >       dsidy = 0.0_dp
560 >       dsidz = 0.0_dp
561 >       dsidux = 0.0_dp
562 >       dsiduy = 0.0_dp
563 >       dsiduz = 0.0_dp
564 >       depsidx = 0.0_dp
565 >       depsidy = 0.0_dp
566 >       depsidz = 0.0_dp
567 >       depsidux = 0.0_dp
568 >       depsiduy = 0.0_dp
569 >       depsiduz = 0.0_dp
570      else
571  
572   #ifdef IS_MPI
# Line 599 | Line 599 | contains  
599  
600         dctidx = - zi * xi / r3
601         dctidy = - zi * yi / r3
602 <       dctidz = 1.0d0 / rij - zi2 / r3
602 >       dctidz = 1.0_dp / rij - zi2 / r3
603         dctidux = yi / rij ! - (zi * xi2) / r3
604         dctiduy = -xi / rij !- (zi * yi2) / r3
605 <       dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
605 >       dctiduz = 0.0_dp !zi / rij - (zi2 * zi) / r3
606  
607         ! this is an attempt to try to truncate the singularity when
608         ! sin(theta) is near 0.0:
609  
610         sti2 = 1.0_dp - cti*cti
611 <       if (dabs(sti2) .lt. 1.0d-12) then
612 <          proji = sqrt(rij * 1.0d-12)
613 <          dcpidx = 1.0d0 / proji
614 <          dcpidy = 0.0d0
611 >       if (abs(sti2) .lt. 1.0e-12_dp) then
612 >          proji = sqrt(rij * 1.0e-12_dp)
613 >          dcpidx = 1.0_dp / proji
614 >          dcpidy = 0.0_dp
615            dcpidux = xi / proji
616 <          dcpiduy = 0.0d0
617 <          dspidx = 0.0d0
618 <          dspidy = 1.0d0 / proji
619 <          dspidux = 0.0d0
616 >          dcpiduy = 0.0_dp
617 >          dspidx = 0.0_dp
618 >          dspidy = 1.0_dp / proji
619 >          dspidux = 0.0_dp
620            dspiduy = yi / proji
621         else
622            proji = sqrt(xi2 + yi2)
623            proji3 = proji*proji*proji
624 <          dcpidx = 1.0d0 / proji - xi2 / proji3
624 >          dcpidx = 1.0_dp / proji - xi2 / proji3
625            dcpidy = - xi * yi / proji3
626            dcpidux = xi / proji - (xi2 * xi) / proji3
627            dcpiduy = - (xi * yi2) / proji3
628            dspidx = - xi * yi / proji3
629 <          dspidy = 1.0d0 / proji - yi2 / proji3
629 >          dspidy = 1.0_dp / proji - yi2 / proji3
630            dspidux = - (yi * xi2) / proji3
631            dspiduy = yi / proji - (yi2 * yi) / proji3
632         endif
633  
634         cpi = xi / proji
635 <       dcpidz = 0.0d0
636 <       dcpiduz = 0.0d0
635 >       dcpidz = 0.0_dp
636 >       dcpiduz = 0.0_dp
637  
638         spi = yi / proji
639 <       dspidz = 0.0d0
640 <       dspiduz = 0.0d0
639 >       dspidz = 0.0_dp
640 >       dspiduz = 0.0_dp
641  
642         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
643              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 648 | Line 648 | contains  
648         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
649              CHEBYSHEV_UN, um_i, dum_i)
650  
651 <       sigma_i = 0.0d0
652 <       s_i = 0.0d0
653 <       eps_i = 0.0d0
654 <       dsigmaidx = 0.0d0
655 <       dsigmaidy = 0.0d0
656 <       dsigmaidz = 0.0d0
657 <       dsigmaidux = 0.0d0
658 <       dsigmaiduy = 0.0d0
659 <       dsigmaiduz = 0.0d0
660 <       dsidx = 0.0d0
661 <       dsidy = 0.0d0
662 <       dsidz = 0.0d0
663 <       dsidux = 0.0d0
664 <       dsiduy = 0.0d0
665 <       dsiduz = 0.0d0
666 <       depsidx = 0.0d0
667 <       depsidy = 0.0d0
668 <       depsidz = 0.0d0
669 <       depsidux = 0.0d0
670 <       depsiduy = 0.0d0
671 <       depsiduz = 0.0d0
651 >       sigma_i = 0.0_dp
652 >       s_i = 0.0_dp
653 >       eps_i = 0.0_dp
654 >       dsigmaidx = 0.0_dp
655 >       dsigmaidy = 0.0_dp
656 >       dsigmaidz = 0.0_dp
657 >       dsigmaidux = 0.0_dp
658 >       dsigmaiduy = 0.0_dp
659 >       dsigmaiduz = 0.0_dp
660 >       dsidx = 0.0_dp
661 >       dsidy = 0.0_dp
662 >       dsidz = 0.0_dp
663 >       dsidux = 0.0_dp
664 >       dsiduy = 0.0_dp
665 >       dsiduz = 0.0_dp
666 >       depsidx = 0.0_dp
667 >       depsidy = 0.0_dp
668 >       depsidz = 0.0_dp
669 >       depsidux = 0.0_dp
670 >       depsiduy = 0.0_dp
671 >       depsiduz = 0.0_dp
672  
673         do lm = 1, ShapeMap%Shapes(st1)%nContactFuncs
674            l = ShapeMap%Shapes(st1)%ContactFuncLValue(lm)
# Line 806 | Line 806 | contains  
806         sigma_j = ShapeMap%Shapes(st2)%sigma
807         s_j = ShapeMap%Shapes(st2)%sigma
808         eps_j = ShapeMap%Shapes(st2)%epsilon
809 <       dsigmajdx = 0.0d0
810 <       dsigmajdy = 0.0d0
811 <       dsigmajdz = 0.0d0
812 <       dsigmajdux = 0.0d0
813 <       dsigmajduy = 0.0d0
814 <       dsigmajduz = 0.0d0
815 <       dsjdx = 0.0d0
816 <       dsjdy = 0.0d0
817 <       dsjdz = 0.0d0
818 <       dsjdux = 0.0d0
819 <       dsjduy = 0.0d0
820 <       dsjduz = 0.0d0
821 <       depsjdx = 0.0d0
822 <       depsjdy = 0.0d0
823 <       depsjdz = 0.0d0
824 <       depsjdux = 0.0d0
825 <       depsjduy = 0.0d0
826 <       depsjduz = 0.0d0
809 >       dsigmajdx = 0.0_dp
810 >       dsigmajdy = 0.0_dp
811 >       dsigmajdz = 0.0_dp
812 >       dsigmajdux = 0.0_dp
813 >       dsigmajduy = 0.0_dp
814 >       dsigmajduz = 0.0_dp
815 >       dsjdx = 0.0_dp
816 >       dsjdy = 0.0_dp
817 >       dsjdz = 0.0_dp
818 >       dsjdux = 0.0_dp
819 >       dsjduy = 0.0_dp
820 >       dsjduz = 0.0_dp
821 >       depsjdx = 0.0_dp
822 >       depsjdy = 0.0_dp
823 >       depsjdz = 0.0_dp
824 >       depsjdux = 0.0_dp
825 >       depsjduy = 0.0_dp
826 >       depsjduz = 0.0_dp
827      else
828  
829   #ifdef IS_MPI
# Line 857 | Line 857 | contains  
857  
858         dctjdx = - zj * xj / r3
859         dctjdy = - zj * yj / r3
860 <       dctjdz = 1.0d0 / rij - zj2 / r3
860 >       dctjdz = 1.0_dp / rij - zj2 / r3
861         dctjdux = yj / rij !- (zi * xj2) / r3
862         dctjduy = -xj / rij !- (zj * yj2) / r3
863 <       dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
863 >       dctjduz = 0.0_dp !zj / rij - (zj2 * zj) / r3
864  
865         ! this is an attempt to try to truncate the singularity when
866         ! sin(theta) is near 0.0:
867  
868         stj2 = 1.0_dp - ctj*ctj
869 <       if (dabs(stj2) .lt. 1.0d-12) then
870 <          projj = sqrt(rij * 1.0d-12)
871 <          dcpjdx = 1.0d0 / projj
872 <          dcpjdy = 0.0d0
869 >       if (abs(stj2) .lt. 1.0e-12_dp) then
870 >          projj = sqrt(rij * 1.0e-12_dp)
871 >          dcpjdx = 1.0_dp / projj
872 >          dcpjdy = 0.0_dp
873            dcpjdux = xj / projj
874 <          dcpjduy = 0.0d0
875 <          dspjdx = 0.0d0
876 <          dspjdy = 1.0d0 / projj
877 <          dspjdux = 0.0d0
874 >          dcpjduy = 0.0_dp
875 >          dspjdx = 0.0_dp
876 >          dspjdy = 1.0_dp / projj
877 >          dspjdux = 0.0_dp
878            dspjduy = yj / projj
879         else
880            projj = sqrt(xj2 + yj2)
881            projj3 = projj*projj*projj
882 <          dcpjdx = 1.0d0 / projj - xj2 / projj3
882 >          dcpjdx = 1.0_dp / projj - xj2 / projj3
883            dcpjdy = - xj * yj / projj3
884            dcpjdux = xj / projj - (xj2 * xj) / projj3
885            dcpjduy = - (xj * yj2) / projj3
886            dspjdx = - xj * yj / projj3
887 <          dspjdy = 1.0d0 / projj - yj2 / projj3
887 >          dspjdy = 1.0_dp / projj - yj2 / projj3
888            dspjdux = - (yj * xj2) / projj3
889            dspjduy = yj / projj - (yj2 * yj) / projj3
890         endif
891  
892         cpj = xj / projj
893 <       dcpjdz = 0.0d0
894 <       dcpjduz = 0.0d0
893 >       dcpjdz = 0.0_dp
894 >       dcpjduz = 0.0_dp
895  
896         spj = yj / projj
897 <       dspjdz = 0.0d0
898 <       dspjduz = 0.0d0
897 >       dspjdz = 0.0_dp
898 >       dspjduz = 0.0_dp
899  
900  
901   !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
# Line 909 | Line 909 | contains  
909         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
910              CHEBYSHEV_UN, um_j, dum_j)
911  
912 <       sigma_j = 0.0d0
913 <       s_j = 0.0d0
914 <       eps_j = 0.0d0
915 <       dsigmajdx = 0.0d0
916 <       dsigmajdy = 0.0d0
917 <       dsigmajdz = 0.0d0
918 <       dsigmajdux = 0.0d0
919 <       dsigmajduy = 0.0d0
920 <       dsigmajduz = 0.0d0
921 <       dsjdx = 0.0d0
922 <       dsjdy = 0.0d0
923 <       dsjdz = 0.0d0
924 <       dsjdux = 0.0d0
925 <       dsjduy = 0.0d0
926 <       dsjduz = 0.0d0
927 <       depsjdx = 0.0d0
928 <       depsjdy = 0.0d0
929 <       depsjdz = 0.0d0
930 <       depsjdux = 0.0d0
931 <       depsjduy = 0.0d0
932 <       depsjduz = 0.0d0
912 >       sigma_j = 0.0_dp
913 >       s_j = 0.0_dp
914 >       eps_j = 0.0_dp
915 >       dsigmajdx = 0.0_dp
916 >       dsigmajdy = 0.0_dp
917 >       dsigmajdz = 0.0_dp
918 >       dsigmajdux = 0.0_dp
919 >       dsigmajduy = 0.0_dp
920 >       dsigmajduz = 0.0_dp
921 >       dsjdx = 0.0_dp
922 >       dsjdy = 0.0_dp
923 >       dsjdz = 0.0_dp
924 >       dsjdux = 0.0_dp
925 >       dsjduy = 0.0_dp
926 >       dsjduz = 0.0_dp
927 >       depsjdx = 0.0_dp
928 >       depsjdy = 0.0_dp
929 >       depsjdz = 0.0_dp
930 >       depsjdux = 0.0_dp
931 >       depsjduy = 0.0_dp
932 >       depsjduz = 0.0_dp
933  
934         do lm = 1, ShapeMap%Shapes(st2)%nContactFuncs
935            l = ShapeMap%Shapes(st2)%ContactFuncLValue(lm)
# Line 1099 | Line 1099 | contains  
1099   !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1100   !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1101   !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1102 <    depsdxi = eps_j * depsidx / (2.0d0 * eps)
1103 <    depsdyi = eps_j * depsidy / (2.0d0 * eps)
1104 <    depsdzi = eps_j * depsidz / (2.0d0 * eps)
1105 <    depsduxi = eps_j * depsidux / (2.0d0 * eps)
1106 <    depsduyi = eps_j * depsiduy / (2.0d0 * eps)
1107 <    depsduzi = eps_j * depsiduz / (2.0d0 * eps)
1102 >    depsdxi = eps_j * depsidx / (2.0_dp * eps)
1103 >    depsdyi = eps_j * depsidy / (2.0_dp * eps)
1104 >    depsdzi = eps_j * depsidz / (2.0_dp * eps)
1105 >    depsduxi = eps_j * depsidux / (2.0_dp * eps)
1106 >    depsduyi = eps_j * depsiduy / (2.0_dp * eps)
1107 >    depsduzi = eps_j * depsiduz / (2.0_dp * eps)
1108  
1109 <    depsdxj = eps_i * depsjdx / (2.0d0 * eps)
1110 <    depsdyj = eps_i * depsjdy / (2.0d0 * eps)
1111 <    depsdzj = eps_i * depsjdz / (2.0d0 * eps)
1112 <    depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1113 <    depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1114 <    depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1109 >    depsdxj = eps_i * depsjdx / (2.0_dp * eps)
1110 >    depsdyj = eps_i * depsjdy / (2.0_dp * eps)
1111 >    depsdzj = eps_i * depsjdz / (2.0_dp * eps)
1112 >    depsduxj = eps_i * depsjdux / (2.0_dp * eps)
1113 >    depsduyj = eps_i * depsjduy / (2.0_dp * eps)
1114 >    depsduzj = eps_i * depsjduz / (2.0_dp * eps)
1115  
1116   !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1117  
# Line 1147 | Line 1147 | contains  
1147      rt12 = rt6*rt6
1148      rt126 = rt12 - rt6
1149  
1150 <    pot_temp = 4.0d0 * eps * rt126
1150 >    pot_temp = 4.0_dp * eps * rt126
1151  
1152      vpair = vpair + pot_temp
1153      if (do_pot) then
1154   #ifdef IS_MPI
1155 <       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5d0*pot_temp*sw
1156 <       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5d0*pot_temp*sw
1155 >       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5_dp*pot_temp*sw
1156 >       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5_dp*pot_temp*sw
1157   #else
1158         pot = pot + pot_temp*sw
1159   #endif
# Line 1161 | Line 1161 | contains  
1161  
1162   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1163  
1164 <    dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1165 <    dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1166 <    dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
1167 <    dvduxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxi + 4.0d0*depsduxi*rt126
1168 <    dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1169 <    dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1164 >    dvdxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxi + 4.0_dp*depsdxi*rt126
1165 >    dvdyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyi + 4.0_dp*depsdyi*rt126
1166 >    dvdzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzi + 4.0_dp*depsdzi*rt126
1167 >    dvduxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxi + 4.0_dp*depsduxi*rt126
1168 >    dvduyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyi + 4.0_dp*depsduyi*rt126
1169 >    dvduzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzi + 4.0_dp*depsduzi*rt126
1170  
1171 <    dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1172 <    dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1173 <    dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
1174 <    dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1175 <    dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1176 <    dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1171 >    dvdxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxj + 4.0_dp*depsdxj*rt126
1172 >    dvdyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyj + 4.0_dp*depsdyj*rt126
1173 >    dvdzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzj + 4.0_dp*depsdzj*rt126
1174 >    dvduxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxj + 4.0_dp*depsduxj*rt126
1175 >    dvduyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyj + 4.0_dp*depsduyj*rt126
1176 >    dvduzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzj + 4.0_dp*depsduzj*rt126
1177   !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1178      ! do the torques first since they are easy:
1179      ! remember that these are still in the body fixed axes
1180  
1181 <    txi = 0.0d0
1182 <    tyi = 0.0d0
1183 <    tzi = 0.0d0
1181 >    txi = 0.0_dp
1182 >    tyi = 0.0_dp
1183 >    tzi = 0.0_dp
1184  
1185 <    txj = 0.0d0
1186 <    tyj = 0.0d0
1187 <    tzj = 0.0d0
1185 >    txj = 0.0_dp
1186 >    tyj = 0.0_dp
1187 >    tzj = 0.0_dp
1188  
1189      txi = (dvduyi - dvduzi) * sw
1190      tyi = (dvduzi - dvduxi) * sw
# Line 1343 | Line 1343 | contains  
1343      end DO
1344  
1345      ! start with 0,0:
1346 <    PLM(0,0) = 1.0D0
1346 >    PLM(0,0) = 1.0_DP
1347  
1348      ! x = +/- 1 functions are easy:
1349 <    IF (abs(X).EQ.1.0D0) THEN
1349 >    IF (abs(X).EQ.1.0_DP) THEN
1350         DO I = 1, m
1351            PLM(0, I) = X**I
1352 <          DLM(0, I) = 0.5D0*I*(I+1.0D0)*X**(I+1)
1352 >          DLM(0, I) = 0.5_DP*I*(I+1.0_DP)*X**(I+1)
1353         end DO
1354         DO J = 1, m
1355            DO I = 1, l
1356               IF (I.EQ.1) THEN
1357                  DLM(I, J) = 1.0D+300
1358               ELSE IF (I.EQ.2) THEN
1359 <                DLM(I, J) = -0.25D0*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1359 >                DLM(I, J) = -0.25_DP*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1360               ENDIF
1361            end DO
1362         end DO
# Line 1364 | Line 1364 | contains  
1364      ENDIF
1365  
1366      LS = 1
1367 <    IF (abs(X).GT.1.0D0) LS = -1
1368 <    XQ = sqrt(LS*(1.0D0-X*X))
1369 <    XS = LS*(1.0D0-X*X)
1367 >    IF (abs(X).GT.1.0_DP) LS = -1
1368 >    XQ = sqrt(LS*(1.0_DP-X*X))
1369 >    XS = LS*(1.0_DP-X*X)
1370  
1371      DO I = 1, l
1372 <       PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1372 >       PLM(I, I) = -LS*(2.0_DP*I-1.0_DP)*XQ*PLM(I-1, I-1)
1373      enddo
1374  
1375      DO I = 0, l
1376 <       PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1376 >       PLM(I, I+1)=(2.0_DP*I+1.0_DP)*X*PLM(I, I)
1377      enddo
1378  
1379      DO I = 0, l
1380         DO J = I+2, m
1381 <          PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1382 <               (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1381 >          PLM(I, J)=((2.0_DP*J-1.0_DP)*X*PLM(I,J-1) - &
1382 >               (I+J-1.0_DP)*PLM(I,J-2))/(J-I)
1383         end DO
1384      end DO
1385  
1386 <    DLM(0, 0)=0.0D0
1386 >    DLM(0, 0)=0.0_DP
1387      DO J = 1, m
1388         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1389      end DO
1390  
1391      DO I = 1, l
1392         DO J = I, m
1393 <          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1393 >          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0_DP)/XQ*PLM(I-1, J)
1394         end DO
1395      end DO
1396  
# Line 1427 | Line 1427 | contains  
1427      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1428      integer :: k
1429  
1430 <    A = 2.0D0
1431 <    B = 0.0D0
1432 <    C = 1.0D0
1433 <    Y0 = 1.0D0
1434 <    Y1 = 2.0D0*X
1435 <    DY0 = 0.0D0
1436 <    DY1 = 2.0D0
1437 <    PL(0) = 1.0D0
1438 <    PL(1) = 2.0D0*X
1439 <    DPL(0) = 0.0D0
1440 <    DPL(1) = 2.0D0
1430 >    A = 2.0_DP
1431 >    B = 0.0_DP
1432 >    C = 1.0_DP
1433 >    Y0 = 1.0_DP
1434 >    Y1 = 2.0_DP*X
1435 >    DY0 = 0.0_DP
1436 >    DY1 = 2.0_DP
1437 >    PL(0) = 1.0_DP
1438 >    PL(1) = 2.0_DP*X
1439 >    DPL(0) = 0.0_DP
1440 >    DPL(1) = 2.0_DP
1441      IF (function_type.EQ.CHEBYSHEV_TN) THEN
1442         Y1 = X
1443 <       DY1 = 1.0D0
1443 >       DY1 = 1.0_DP
1444         PL(1) = X
1445 <       DPL(1) = 1.0D0
1445 >       DPL(1) = 1.0_DP
1446      ELSE IF (function_type.EQ.LAGUERRE) THEN
1447 <       Y1 = 1.0D0-X
1448 <       DY1 = -1.0D0
1449 <       PL(1) = 1.0D0-X
1450 <       DPL(1) = -1.0D0
1447 >       Y1 = 1.0_DP-X
1448 >       DY1 = -1.0_DP
1449 >       PL(1) = 1.0_DP-X
1450 >       DPL(1) = -1.0_DP
1451      ENDIF
1452      DO K = 2, m
1453         IF (function_type.EQ.LAGUERRE) THEN
1454 <          A = -1.0D0/K
1455 <          B = 2.0D0+A
1456 <          C = 1.0D0+A
1454 >          A = -1.0_DP/K
1455 >          B = 2.0_DP+A
1456 >          C = 1.0_DP+A
1457         ELSE IF (function_type.EQ.HERMITE) THEN
1458 <          C = 2.0D0*(K-1.0D0)
1458 >          C = 2.0_DP*(K-1.0_DP)
1459         ENDIF
1460         YN = (A*X+B)*Y1-C*Y0
1461         DYN = A*Y1+(A*X+B)*DY1-C*DY0

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines