50 |
|
!! @author Matthew Meineke |
51 |
|
!! @author Christopher Fennell |
52 |
|
!! @author J. Daniel Gezelter |
53 |
< |
!! @version $Id: sticky.F90,v 1.9 2005-05-12 19:43:48 chrisfen Exp $, $Date: 2005-05-12 19:43:48 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $ |
53 |
> |
!! @version $Id: sticky.F90,v 1.10 2005-05-17 22:35:01 chrisfen Exp $, $Date: 2005-05-17 22:35:01 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $ |
54 |
|
|
55 |
|
module sticky |
56 |
|
|
550 |
|
integer :: me1, me2 |
551 |
|
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig |
552 |
|
real (kind=dp) :: zi3, zi4, zi5, zj3, zj4, zj5 |
553 |
< |
real (kind=dp) :: oSw1, oSw2, prodVal |
553 |
> |
real (kind=dp) :: frac1, frac2, prodVal |
554 |
|
real (kind=dp) :: prei1, prei2, prei, prej1, prej2, prej |
555 |
|
real (kind=dp) :: walt, walti, waltj, dwaltidx, dwaltidy, dwaltidz |
556 |
|
real (kind=dp) :: dwaltjdx, dwaltjdy, dwaltjdz |
608 |
|
rI4 = rI2*rI2 |
609 |
|
rI5 = rI3*rI2 |
610 |
|
rI6 = rI3*rI3 |
611 |
< |
rI7 = rI5*rI2 |
611 |
> |
rI7 = rI4*rI3 |
612 |
|
|
613 |
|
drdx = d(1) * rI |
614 |
|
drdy = d(2) * rI |
664 |
|
|
665 |
|
call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr) |
666 |
|
|
667 |
< |
wi = 2.0d0*(xi2-yi2)*zi * rI3 |
668 |
< |
wj = 2.0d0*(xj2-yj2)*zj * rI3 |
667 |
> |
frac1 = 0.5d0 |
668 |
> |
frac2 = 0.5d0 |
669 |
|
|
670 |
+ |
wi = 2.0d0*(xi2-yi2)*zi*rI3 |
671 |
+ |
wj = 2.0d0*(xj2-yj2)*zj*rI3 |
672 |
+ |
|
673 |
|
! prodVal = zihat*zjhat |
674 |
|
! if (prodVal .ge. 0.0d0) then |
675 |
|
! wi = 0.0d0 |
679 |
|
wi2 = wi*wi |
680 |
|
wj2 = wj*wj |
681 |
|
|
682 |
< |
w = wi*wi2+wj*wj2 |
682 |
> |
w = frac1*wi*wi2 + frac2*wi + wj*wj2 |
683 |
|
|
684 |
|
zif = zihat - 0.6d0 |
685 |
|
zis = zihat + 0.8d0 |
709 |
|
#endif |
710 |
|
endif |
711 |
|
|
712 |
< |
dwidx = 3.0d0*wi2*( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 ) |
713 |
< |
dwidy = 3.0d0*wi2*( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 ) |
714 |
< |
dwidz = 3.0d0*wi2*( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 ) |
712 |
> |
dwidx = ( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 ) |
713 |
> |
dwidy = ( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 ) |
714 |
> |
dwidz = ( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 ) |
715 |
> |
|
716 |
> |
dwidx = frac1*3.0d0*wi2*dwidx + frac2*dwidx |
717 |
> |
dwidy = frac1*3.0d0*wi2*dwidy + frac2*dwidy |
718 |
> |
dwidz = frac1*3.0d0*wi2*dwidz + frac2*dwidz |
719 |
|
|
720 |
< |
dwjdx = 3.0d0*wj2*( 4.0d0*xj*zj*rI3 - 6.0d0*xj*zj*(xj2-yj2)*rI5 ) |
721 |
< |
dwjdy = 3.0d0*wj2*( -4.0d0*yj*zj*rI3 - 6.0d0*yj*zj*(xj2-yj2)*rI5 ) |
722 |
< |
dwjdz = 3.0d0*wj2*( 2.0d0*(xj2-yj2)*rI3 - 6.0d0*zj2*(xj2-yj2)*rI5 ) |
720 |
> |
dwjdx = ( 4.0d0*xj*zj*rI3 - 6.0d0*xj*zj*(xj2-yj2)*rI5 ) |
721 |
> |
dwjdy = ( -4.0d0*yj*zj*rI3 - 6.0d0*yj*zj*(xj2-yj2)*rI5 ) |
722 |
> |
dwjdz = ( 2.0d0*(xj2-yj2)*rI3 - 6.0d0*zj2*(xj2-yj2)*rI5 ) |
723 |
|
|
724 |
+ |
dwjdx = frac1*3.0d0*wj2*dwjdx + frac2*dwjdx |
725 |
+ |
dwjdy = frac1*3.0d0*wj2*dwjdy + frac2*dwjdy |
726 |
+ |
dwjdz = frac1*3.0d0*wj2*dwjdz + frac2*dwjdz |
727 |
+ |
|
728 |
|
uglyi = zif*zif*zis + zif*zis*zis |
729 |
|
uglyj = zjf*zjf*zjs + zjf*zjs*zjs |
730 |
|
|
752 |
|
!dwjpdy = 0.0d0 |
753 |
|
!dwjpdz = 0.0d0 |
754 |
|
|
755 |
< |
dwidux = 3.0d0*wi2*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 ) |
756 |
< |
dwiduy = 3.0d0*wi2*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 ) |
757 |
< |
dwiduz = 3.0d0*wi2*( -8.0d0*xi*yi*zi*rI3 ) |
755 |
> |
dwidux = ( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 ) |
756 |
> |
dwiduy = ( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 ) |
757 |
> |
dwiduz = ( -8.0d0*xi*yi*zi*rI3 ) |
758 |
|
|
759 |
< |
dwjdux = 3.0d0*wj2*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 ) |
760 |
< |
dwjduy = 3.0d0*wj2*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 ) |
761 |
< |
dwjduz = 3.0d0*wj2*( -8.0d0*xj*yj*zj*rI3 ) |
759 |
> |
dwidux = frac1*3.0d0*wi2*dwidux + frac2*dwidux |
760 |
> |
dwiduy = frac1*3.0d0*wi2*dwiduy + frac2*dwiduy |
761 |
> |
dwiduz = frac1*3.0d0*wi2*dwiduz + frac2*dwiduz |
762 |
|
|
763 |
+ |
dwjdux = ( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 ) |
764 |
+ |
dwjduy = ( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 ) |
765 |
+ |
dwjduz = ( -8.0d0*xj*yj*zj*rI3 ) |
766 |
+ |
|
767 |
+ |
dwjdux = frac1*3.0d0*wj2*dwjdux + frac2*dwjdux |
768 |
+ |
dwjduy = frac1*3.0d0*wj2*dwjduy + frac2*dwjduy |
769 |
+ |
dwjduz = frac1*3.0d0*wj2*dwjduz + frac2*dwjduz |
770 |
+ |
|
771 |
|
dwipdux = 2.0d0*yi*uglyi*rI |
772 |
|
dwipduy = -2.0d0*xi*uglyi*rI |
773 |
|
dwipduz = 0.0d0 |