mod_DFT_DISP_WDA_d.F90

Go to the documentation of this file.

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
!
MODULE mod_dft_disp_wda_d
  IMPLICIT NONE
  PRIVATE 
  PUBLIC disp_weighted_densities_d
  PUBLIC disp_mu_d
  PUBLIC disp_dfdrho_wda_d

 CONTAINS

 !  Differentiation of disp_weighted_densities in forward (tangent) mode:
!   variations   of useful results: rhop_wd
!   with respect to varying inputs: rhop
  SUBROUTINE disp_weighted_densities_d(rhop, rhopd, rhop_wd, rhop_wdd, &
&   user)
    Use petscmanagement
    USE basic_variables, ONLY : ncomp
    USE mod_dft, ONLY : fa_disp, ab_disp, zp, dzp
    IMPLICIT NONE

#include <finclude/petscsys.h>

!passed
type(userctx)    :: user
petscscalar       :: rhop(ncomp,user%gxs:user%gxe)
petscscalar       :: rhopd(ncomp,user%gxs:user%gxe)
REAL, INTENT(OUT) :: rhop_wd(user%gxs:user%gxe,ncomp)
REAL, INTENT(OUT) :: rhop_wdd(user%gxs:user%gxe,ncomp)

    
! ! !local
! !     INTEGER :: i, j, k
! !     INTEGER :: n
! ! !Fehlermeldung einbauen, falls dim > 400!! 
! !     REAL :: x_int(400), y_int(400), y2(400)
! !     REAL :: y_intd(400), y2d(400)
! !     REAL :: zmin, dz, zz
! !     REAL :: xlo, xhi, int1
! !     REAL :: int1d
! !     INTEGER :: fa_disp_max
! !     INTRINSIC MAXVAL
! !     zmin = 1d-6
! !     fa_disp_max = MAXVAL(fa_disp(1:ncomp))
! !     rhop_wdd = 0.0
! !     y2d = 0.0
! !     DO k=1,ncomp
! !         Do i = user%xs - fa_disp_max , user%xe + fa_disp_max
! !         n = 1
! !         x_int = 0.0
! !         y_int = 0.0
! !         y_intd = 0.0
! !         DO j=i-fa_disp(k),i+fa_disp(k)
! !           IF (zp(i) - zp(j+1) .LT. ab_disp(k) .AND. zp(i) - zp(j) .GE. &
! ! &             ab_disp(k)) THEN
! ! !the position of j+1 is already within i-d/2 while j is still outside this range in this case, the integration steplength (dz) is
! ! ! just the distance, which j+1 overlaps with i-d/2 and what is integrated is the interpolated value of the integrand
! ! !here always n=1!
! !             IF (n .NE. 1) THEN
! !               GOTO 100
! !             ELSE
! ! !distance between grid points j and i
! !               zz = zp(j) - zp(i)
! ! !the part of the intervall between zp(j) and zp(j+1) which is already within i-d/2
! !               dz = zp(j+1) - (zp(i)-ab_disp(k))
! ! !if(dz < epsilon(dz)) dz = epsilon(dz) !bei unguenstiger Kombination von sig und ngrid kann dz unter Machinengenauigkeit epsilon 
! ! !liegen, dann ist x(2) = x(1) + dz = x(1) -> das fuehrt zu Abbruch in Spline Interpolation           
! ! !array containing x-values for spline integration 
! !               x_int(n) = 0.0
! !               y_intd(n) = 0.0
! !               y_int(n) = 0.0
! !             END IF
! !           ELSE IF (zp(j) .GT. zp(i) - ab_disp(k) .AND. zp(j) .LE. zp(i) &
! ! &             + ab_disp(k)) THEN
! ! !grid point j  within i+-d2 
! !             n = n + 1
! ! !first time in this If condition, dz is stil the old value from above!
! !             x_int(n) = x_int(n-1) + dz
! !             zz = zp(j) - zp(i)
! !             dz = dzp
! !             y_intd(n) = (ab_disp(k)*ab_disp(k)-zz*zz)*rhopd(k, j)
! !             y_int(n) = rhop(k, j)*(ab_disp(k)*ab_disp(k)-zz*zz)
! !             IF (zp(j) .LT. zp(i) + ab_disp(k) .AND. zp(j+1) .GE. zp(i) +&
! ! &               ab_disp(k)) THEN
! ! !zp(j) is still within zp(i)+d2 but zp(j+1) is already outside zp(i)+d2
! !               dz = zp(i) + ab_disp(k) - zp(j)
! ! !If(dz <= epsilon(dz)) exit   !wie oben, kann auch hier bei ungluecklicher Wahl von sig und ngrid dz < eps werden und somit x(n) 
! ! != x(n-1) -> Abbruch in Spline interpolation. Dann einfach ngrid aendern!
! !               zz = zp(j) - zp(i)
! !               n = n + 1
! !               x_int(n) = x_int(n-1) + dz
! !               y_intd(n) = 0.0
! !               y_int(n) = 0.0
! !             END IF
! !           END IF
! !         END DO
! !         xlo = x_int(1)
! !         xhi = x_int(n)
! !         CALL SPLINE_D(x_int(1:n), y_int(1:n), y_intd(1:n), n, 1.e30, &
! ! &               1.e30, y2(1:n), y2d(1:n))
! !         CALL SPLINT_INTEGRAL_D(x_int(1:n), y_int(1:n), y_intd(1:n), y2(1&
! ! &                        :n), y2d(1:n), n, xlo, xhi, int1, int1d)
! !         rhop_wdd(i, k) = 0.75*int1d/ab_disp(k)**3
! !         rhop_wd(i, k) = 0.75*int1/ab_disp(k)**3
! !         IF (rhop_wd(i, k) .LT. 0.0) THEN
! !           rhop_wdd(i, k) = 0.0
! !           rhop_wd(i, k) = 0.0
! !           DO j=2,n
! !             rhop_wdd(i, k) = rhop_wdd(i, k) + (x_int(j)-x_int(j-1))*(&
! ! &             y_intd(j)+y_intd(j-1))/2.0
! !             rhop_wd(i, k) = rhop_wd(i, k) + (y_int(j)+y_int(j-1))/2.0*(&
! ! &             x_int(j)-x_int(j-1))
! !           END DO
! !         END IF
! !         IF (rhop_wd(i, k) .LT. 0.0) THEN
! !           rhop_wdd(i, k) = rhopd(k, i)
! !           rhop_wd(i, k) = rhop(k, i)
! !         END IF
! !       END DO
! !     END DO
! !     GOTO 110
! !  100 STOP
! !  110 CONTINUE
  END SUBROUTINE disp_weighted_densities_d

  
  
!  Differentiation of disp_mu in forward (tangent) mode:
!   variations   of useful results: my_disp
!   with respect to varying inputs: rhop_wd
  SUBROUTINE disp_mu_d(rhop_wd, rhop_wdd, f_disp, my_disp, my_dispd, &
&   df_disp_drk, user)
    
    Use petscmanagement
    USE parameters, ONLY : pi
    USE eos_constants, ONLY : ap, bp
    USE basic_variables, ONLY : ncomp, t, parame
    USE eos_variables, Only: dhs, sig_ij, uij
    USE mod_dft, ONLY : fa_disp
    IMPLICIT NONE
#include <finclude/petscsys.h>

!passed
type(userctx)    :: user
REAL, INTENT(IN)  :: rhop_wd(user%gxs:user%gxe,ncomp)
REAL, INTENT(IN)  :: rhop_wdd(user%gxs:user%gxe,ncomp)
REAL, INTENT(OUT) :: my_disp(user%gxs:user%gxe,ncomp) !chemPot_disp / kT
REAL, INTENT(OUT) :: my_dispd(user%gxs:user%gxe,ncomp) !chemPot_disp / kT
REAL, INTENT(OUT) :: f_disp(user%gxs:user%gxe)  !F_disp / NkT
REAL, INTENT(OUT) :: df_disp_drk(user%gxs:user%gxe,ncomp) ! d(F/NkT) / drho_k = (mu/kT - F_disp / NkT)/rho                    

! ! !local
! !     INTEGER :: k, ii, kk, m
! !     REAL :: m_mean, m_rk(ncomp)
! !     REAL :: m_meand, m_rkd(ncomp)
! !     REAL :: apar(0:6), bpar(0:6)
! !     REAL :: apard(0:6), bpard(0:6)
! !     REAL :: ap_rk(ncomp, 0:6), bp_rk(ncomp, 0:6)
! !     REAL :: ap_rkd(ncomp, 0:6), bp_rkd(ncomp, 0:6)
! !     REAL :: xi(ncomp), z3
! !     REAL :: xid(ncomp), z3d
! !     REAL :: i1, i2, i1_rk, i2_rk
! !     REAL :: i1d, i2d, i1_rkd, i2_rkd
! !     REAL :: order1, order2, ord1_rk, ord2_rk
! !     REAL :: order1d, order2d, ord1_rkd, ord2_rkd
! !     REAL :: c1_con, c2_con, c1_rk, rho2
! !     REAL :: c1_cond, c2_cond, c1_rkd, rho2d
! !     REAL :: rhop_wd_sum, eta_disp, eta_rk, zms
! !     REAL :: rhop_wd_sumd, eta_dispd, zmsd
! !     INTEGER :: fa_disp_max
! !     INTRINSIC MAXVAL
! !     INTRINSIC SUM
! !     INTRINSIC REAL
! !     REAL :: pwy1
! !     REAL :: pwr1
! !     REAL :: pwr1d
! !     REAL :: pwy2
! !     REAL :: pwr2
! !     REAL :: pwr2d
! !     fa_disp_max = MAXVAL(fa_disp(1:ncomp))
! !     my_dispd = 0.0
! !     xid = 0.0
! !     m_rkd = 0.0
! !     ap_rkd = 0.0
! !     bpard = 0.0
! !     apard = 0.0
! !     bp_rkd = 0.0
! !     DO k=1,ncomp
! !         Do ii = user%xs-fa_disp_max,user%xe+fa_disp_max
! !         rhop_wd_sumd = SUM(rhop_wdd(ii, 1:ncomp))
! !         rhop_wd_sum = SUM(rhop_wd(ii, 1:ncomp))
! !         m_mean = 0.0
! !         eta_disp = 0.0
! !         eta_dispd = 0.0
! !         m_meand = 0.0
! !         DO kk=1,ncomp
! !           xid(kk) = (rhop_wdd(ii, kk)*rhop_wd_sum-rhop_wd(ii, kk)*&
! ! &           rhop_wd_sumd)/rhop_wd_sum**2
! !           xi(kk) = rhop_wd(ii, kk)/rhop_wd_sum
! !           m_meand = m_meand + parame(kk,1)*xid(kk)
! !           m_mean = m_mean + xi(kk)*parame(kk,1)
! !           eta_dispd = eta_dispd + parame(kk,1)*dhs(kk)**3*rhop_wdd(ii, kk)
! !           eta_disp = eta_disp + rhop_wd(ii, kk)*parame(kk,1)*dhs(kk)**3
! !         END DO
! !         eta_dispd = pi*eta_dispd/6.0
! !         eta_disp = eta_disp*pi/6.0
! !         eta_rk = parame(k,1)*dhs(k)**3*pi/6.0
! !         m_rkd(k) = (-(m_meand*rhop_wd_sum)-(parame(k,1)-m_mean)*rhop_wd_sumd&
! ! &         )/rhop_wd_sum**2
! !         m_rk(k) = (parame(k,1)-m_mean)/rhop_wd_sum
! !         DO m=0,6
! !           apard(m) = ap(m, 2)*m_meand/m_mean**2 + ap(m, 3)*(m_meand*(1.0&
! ! &           -2.0/m_mean)/m_mean**2+(1.0-1.0/m_mean)*2.0*m_meand/m_mean**&
! ! &           2)
! !           apar(m) = ap(m, 1) + (1.0-1.0/m_mean)*ap(m, 2) + (1.0-1.0/&
! ! &           m_mean)*(1.0-2.0/m_mean)*ap(m, 3)
! !           bpard(m) = bp(m, 2)*m_meand/m_mean**2 + bp(m, 3)*(m_meand*(1.0&
! ! &           -2.0/m_mean)/m_mean**2+(1.0-1.0/m_mean)*2.0*m_meand/m_mean**&
! ! &           2)
! !           bpar(m) = bp(m, 1) + (1.0-1.0/m_mean)*bp(m, 2) + (1.0-1.0/&
! ! &           m_mean)*(1.0-2.0/m_mean)*bp(m, 3)
! ! ! --- derivatives of apar, bpar to rho_k ---------------------------
! !           ap_rkd(k, m) = (m_rkd(k)*m_mean**2-m_rk(k)*2*m_mean*m_meand)*(&
! ! &           ap(m, 2)+(3.0-4.0/m_mean)*ap(m, 3))/m_mean**4 + m_rk(k)*ap(m&
! ! &           , 3)*4.0*m_meand/m_mean**4
! !           ap_rk(k, m) = m_rk(k)/m_mean**2*(ap(m, 2)+(3.0-4.0/m_mean)*ap(&
! ! &           m, 3))
! !           bp_rkd(k, m) = (m_rkd(k)*m_mean**2-m_rk(k)*2*m_mean*m_meand)*(&
! ! &           bp(m, 2)+(3.0-4.0/m_mean)*bp(m, 3))/m_mean**4 + m_rk(k)*bp(m&
! ! &           , 3)*4.0*m_meand/m_mean**4
! !           bp_rk(k, m) = m_rk(k)/m_mean**2*(bp(m, 2)+(3.0-4.0/m_mean)*bp(&
! ! &           m, 3))
! !         END DO
! !         i1 = 0.0
! !         i2 = 0.0
! !         i1_rk = 0.0
! !         i2_rk = 0.0
! !         i1_rkd = 0.0
! !         i1d = 0.0
! !         i2d = 0.0
! !         i2_rkd = 0.0
! !         DO m=0,6
! !           pwy1 = REAL(m)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy1 .EQ. &
! ! &             INT(pwy1))) THEN
! !             pwr1d = pwy1*eta_disp**(pwy1-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy1 .EQ. 1.0) THEN
! !             pwr1d = eta_dispd
! !           ELSE
! !             pwr1d = 0.0
! !           END IF
! !           pwr1 = eta_disp**pwy1
! !           i1d = i1d + apard(m)*pwr1 + apar(m)*pwr1d
! !           i1 = i1 + apar(m)*pwr1
! !           pwy1 = REAL(m)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy1 .EQ. &
! ! &             INT(pwy1))) THEN
! !             pwr1d = pwy1*eta_disp**(pwy1-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy1 .EQ. 1.0) THEN
! !             pwr1d = eta_dispd
! !           ELSE
! !             pwr1d = 0.0
! !           END IF
! !           pwr1 = eta_disp**pwy1
! !           i2d = i2d + bpard(m)*pwr1 + bpar(m)*pwr1d
! !           i2 = i2 + bpar(m)*pwr1
! !           pwy1 = REAL(m - 1)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy1 .EQ. &
! ! &             INT(pwy1))) THEN
! !             pwr1d = pwy1*eta_disp**(pwy1-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy1 .EQ. 1.0) THEN
! !             pwr1d = eta_dispd
! !           ELSE
! !             pwr1d = 0.0
! !           END IF
! !           pwr1 = eta_disp**pwy1
! !           pwy2 = REAL(m)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy2 .EQ. &
! ! &             INT(pwy2))) THEN
! !             pwr2d = pwy2*eta_disp**(pwy2-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy2 .EQ. 1.0) THEN
! !             pwr2d = eta_dispd
! !           ELSE
! !             pwr2d = 0.0
! !           END IF
! !           pwr2 = eta_disp**pwy2
! !           i1_rkd = i1_rkd + REAL(m)*eta_rk*(apard(m)*pwr1+apar(m)*pwr1d)&
! ! &           + ap_rkd(k, m)*pwr2 + ap_rk(k, m)*pwr2d
! !           i1_rk = i1_rk + apar(m)*REAL(m)*pwr1*eta_rk + ap_rk(k, m)*pwr2
! !           pwy1 = REAL(m - 1)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy1 .EQ. &
! ! &             INT(pwy1))) THEN
! !             pwr1d = pwy1*eta_disp**(pwy1-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy1 .EQ. 1.0) THEN
! !             pwr1d = eta_dispd
! !           ELSE
! !             pwr1d = 0.0
! !           END IF
! !           pwr1 = eta_disp**pwy1
! !           pwy2 = REAL(m)
! !           IF (eta_disp .GT. 0.0 .OR. (eta_disp .LT. 0.0 .AND. pwy2 .EQ. &
! ! &             INT(pwy2))) THEN
! !             pwr2d = pwy2*eta_disp**(pwy2-1)*eta_dispd
! !           ELSE IF (eta_disp .EQ. 0.0 .AND. pwy2 .EQ. 1.0) THEN
! !             pwr2d = eta_dispd
! !           ELSE
! !             pwr2d = 0.0
! !           END IF
! !           pwr2 = eta_disp**pwy2
! !           i2_rkd = i2_rkd + REAL(m)*eta_rk*(bpard(m)*pwr1+bpar(m)*pwr1d)&
! ! &           + bp_rkd(k, m)*pwr2 + bp_rk(k, m)*pwr2d
! !           i2_rk = i2_rk + bpar(m)*REAL(m)*pwr1*eta_rk + bp_rk(k, m)*pwr2
! !         END DO
! !         ord1_rk = 0.0
! !         ord2_rk = 0.0
! !         order1 = 0.0
! !         order2 = 0.0
! !         order1d = 0.0
! !         order2d = 0.0
! !         ord2_rkd = 0.0
! !         ord1_rkd = 0.0
! !         DO kk=1,ncomp
! ! !sig_ij(kk,k) = 0.5 * ( dhs(kk) + dhs(k) )
! ! !uij(kk,k)    = (1.0 - kij(kk,k)) * SQRT( eps(kk) * eps(k) )
! !           order1d = order1d + parame(kk,1)*parame(k,1)*sig_ij(kk, k)**3*uij(kk, &
! ! &           k)*(xid(kk)*xi(k)+xi(kk)*xid(k))/t
! !           order1 = order1 + xi(kk)*xi(k)*parame(kk,1)*parame(k,1)*sig_ij(kk, k)&
! ! &           **3*uij(kk, k)/t
! !           order2d = order2d + parame(kk,1)*parame(k,1)*sig_ij(kk, k)**3*uij(kk, &
! ! &           k)**2*(xid(kk)*xi(k)+xi(kk)*xid(k))/t**2
! !           order2 = order2 + xi(kk)*xi(k)*parame(kk,1)*parame(k,1)*sig_ij(kk, k)&
! ! &           **3*(uij(kk, k)/t)**2
! !           ord1_rkd = ord1_rkd + 2.0*parame(k,1)*parame(kk,1)*sig_ij(kk, k)**3*&
! ! &           uij(kk, k)*(rhop_wd_sumd*xi(kk)+rhop_wd_sum*xid(kk))/t
! !           ord1_rk = ord1_rk + 2.0*parame(k,1)*rhop_wd_sum*xi(kk)*parame(kk,1)*&
! ! &           sig_ij(kk, k)**3*uij(kk, k)/t
! !           ord2_rkd = ord2_rkd + 2.0*parame(k,1)*parame(kk,1)*sig_ij(kk, k)**3*&
! ! &           uij(kk, k)**2*(rhop_wd_sumd*xi(kk)+rhop_wd_sum*xid(kk))/t**2
! !           ord2_rk = ord2_rk + 2.0*parame(k,1)*rhop_wd_sum*xi(kk)*parame(kk,1)*&
! ! &           sig_ij(kk, k)**3*(uij(kk, k)/t)**2
! !         END DO
! !         z3d = eta_dispd
! !         z3 = eta_disp
! !         zmsd = -z3d
! !         zms = 1.0 - z3
! !         c1_cond = -((((m_meand*(8.0*z3-2.0*z3*z3)+m_mean*(8.0*z3d-2.0*(&
! ! &         z3d*z3+z3*z3d)))*zms**4-m_mean*(8.0*z3-2.0*z3*z3)*4*zms**3*&
! ! &         zmsd)/(zms**4)**2+(((1.0-m_mean)*(20.0*z3d-27.0*(z3d*z3+z3*z3d&
! ! &         )+12.0*3*z3**2*z3d-2.0*4*z3**3*z3d)-m_meand*(20.0*z3-27.0*z3*&
! ! &         z3+12.0*z3**3-2.0*z3**4))*zms**2*(2.0-z3)**2-(1.0-m_mean)*(&
! ! &         20.0*z3-27.0*z3*z3+12.0*z3**3-2.0*z3**4)*2*zms*(2.0-z3)*(zmsd*&
! ! &         (2.0-z3)-zms*z3d))/((zms*(2.0-z3))**2)**2)/(1.0+m_mean*(8.0*z3&
! ! &         -2.0*z3*z3)/zms**4+(1.0-m_mean)*(20.0*z3-27.0*z3*z3+12.0*z3**3&
! ! &         -2.0*z3**4)/(zms*(2.0-z3))**2)**2)
! !         c1_con = 1.0/(1.0+m_mean*(8.0*z3-2.0*z3*z3)/zms**4+(1.0-m_mean)*&
! ! &         (20.0*z3-27.0*z3*z3+12.0*z3**3-2.0*z3**4)/(zms*(2.0-z3))**2)
! !         c2_cond = -((c1_cond*c1_con+c1_con*c1_cond)*(m_mean*(-(4.0*z3*z3&
! ! &         )+20.0*z3+8.0)/zms**5+(1.0-m_mean)*(2.0*z3**3+12.0*z3*z3-48.0*&
! ! &         z3+40.0)/(zms*(2.0-z3))**3)+c1_con**2*(((m_meand*(-(4.0*z3*z3)&
! ! &         +20.0*z3+8.0)+m_mean*(20.0*z3d-4.0*(z3d*z3+z3*z3d)))*zms**5-&
! ! &         m_mean*(-(4.0*z3*z3)+20.0*z3+8.0)*5*zms**4*zmsd)/(zms**5)**2+(&
! ! &         ((1.0-m_mean)*(2.0*3*z3**2*z3d+12.0*(z3d*z3+z3*z3d)-48.0*z3d)-&
! ! &         m_meand*(2.0*z3**3+12.0*z3*z3-48.0*z3+40.0))*zms**3*(2.0-z3)**&
! ! &         3-(1.0-m_mean)*(2.0*z3**3+12.0*z3*z3-48.0*z3+40.0)*3*zms**2*(&
! ! &         2.0-z3)**2*(zmsd*(2.0-z3)-zms*z3d))/((zms*(2.0-z3))**3)**2))
! !         c2_con = -(c1_con*c1_con*(m_mean*(-(4.0*z3*z3)+20.0*z3+8.0)/zms&
! ! &         **5+(1.0-m_mean)*(2.0*z3**3+12.0*z3*z3-48.0*z3+40.0)/(zms*(2.0&
! ! &         -z3))**3))
! !         c1_rkd = eta_rk*c2_cond - ((c1_cond*c1_con+c1_con*c1_cond)*m_rk(&
! ! &         k)+c1_con**2*m_rkd(k))*((8.0*z3-2.0*z3*z3)/zms**4-(-(2.0*z3**4&
! ! &         )+12.0*z3**3-27.0*z3*z3+20.0*z3)/(zms*(2.0-z3))**2) - c1_con**&
! ! &         2*m_rk(k)*(((8.0*z3d-2.0*(z3d*z3+z3*z3d))*zms**4-(8.0*z3-2.0*&
! ! &         z3*z3)*4*zms**3*zmsd)/(zms**4)**2-((12.0*3*z3**2*z3d-2.0*4*z3&
! ! &         **3*z3d-27.0*(z3d*z3+z3*z3d)+20.0*z3d)*zms**2*(2.0-z3)**2-(-(&
! ! &         2.0*z3**4)+12.0*z3**3-27.0*z3*z3+20.0*z3)*2*zms*(2.0-z3)*(zmsd&
! ! &         *(2.0-z3)-zms*z3d))/((zms*(2.0-z3))**2)**2)
! !         c1_rk = c2_con*eta_rk - c1_con*c1_con*m_rk(k)*((8.0*z3-2.0*z3*z3&
! ! &         )/zms**4-(-(2.0*z3**4)+12.0*z3**3-27.0*z3*z3+20.0*z3)/(zms*(&
! ! &         2.0-z3))**2)
! !         rho2d = rhop_wd_sumd*rhop_wd_sum + rhop_wd_sum*rhop_wd_sumd
! !         rho2 = rhop_wd_sum*rhop_wd_sum
! !         my_dispd(ii, k) = -(2.0*pi*((order1d*rho2+order1*rho2d)*i1_rk+&
! ! &         order1*rho2*i1_rkd+ord1_rkd*i1+ord1_rk*i1d)) - pi*((c1_cond*&
! ! &         m_mean+c1_con*m_meand)*(order2*rho2*i2_rk+ord2_rk*i2)+c1_con*&
! ! &         m_mean*((order2d*rho2+order2*rho2d)*i2_rk+order2*rho2*i2_rkd+&
! ! &         ord2_rkd*i2+ord2_rk*i2d)) - pi*((c1_cond*m_rk(k)+c1_con*m_rkd(&
! ! &         k)+c1_rkd*m_mean+c1_rk*m_meand)*order2*rho2*i2+(c1_con*m_rk(k)&
! ! &         +c1_rk*m_mean)*((order2d*rho2+order2*rho2d)*i2+order2*rho2*i2d&
! ! &         ))
! !         my_disp(ii, k) = -(2.0*pi*(order1*rho2*i1_rk+ord1_rk*i1)) - pi*&
! ! &         c1_con*m_mean*(order2*rho2*i2_rk+ord2_rk*i2) - pi*(c1_con*m_rk&
! ! &         (k)+c1_rk*m_mean)*order2*rho2*i2
! !         f_disp(ii) = -(2.0*pi*rhop_wd_sum*i1*order1) - pi*rhop_wd_sum*&
! ! &         c1_con*m_mean*i2*order2
! !         df_disp_drk(ii, k) = (my_disp(ii, k)-f_disp(ii))/rhop_wd_sum
! !       END DO
! !     END DO
  END SUBROUTINE disp_mu_d

  
  
!  Differentiation of disp_dfdrho_wda in forward (tangent) mode:
!   variations   of useful results: df_drho_disp
!   with respect to varying inputs: my_disp df_drho_disp
  SUBROUTINE disp_dfdrho_wda_d(ii, rhop, rhop_wd, my_disp, my_dispd, &
&   f_disp, df_disp_drk, df_drho_disp, df_drho_dispd, user)

    Use petscmanagement
    USE basic_variables, ONLY : ncomp
    USE mod_dft, ONLY : zp, dzp, fa_disp, ab_disp
    IMPLICIT NONE
#include <finclude/petscsys.h>

!passed
INTEGER, INTENT(IN) :: ii
type(userctx)      :: user
petscscalar         :: rhop(ncomp,user%gxs:user%gxe)
REAL, INTENT(IN)    :: rhop_wd(user%gxs:user%gxe,ncomp)
REAL, INTENT(IN)    :: my_disp(user%gxs:user%gxe,ncomp) !chemPot_disp / kT
REAL, INTENT(IN)    :: my_dispd(user%gxs:user%gxe,ncomp) !chemPot_disp / kT
REAL, INTENT(IN)    :: f_disp(user%gxs:user%gxe)  !F_disp / NkT
REAL, INTENT(IN)    :: df_disp_drk(user%gxs:user%gxe,ncomp) ! d(F/NkT) / drho_k = (mu/kT - F_disp / NkT)/rho                   
REAL, INTENT(OUT)   :: df_drho_disp(ncomp)
REAL, INTENT(OUT)   :: df_drho_dispd(ncomp)

! ! !local
! !     INTEGER :: n, icomp, jj
! !     REAL :: zmin, zz, dz
! !     REAL :: x_int(400), y_int(400), y2(400)
! !     REAL :: y_intd(400), y2d(400)
! !     REAL :: xhi, xlo, int2
! !     REAL :: int2d
! !     REAL :: rhop_sum
! !     zmin = 1d-6
! !     y2d = 0.0
! !     DO icomp=1,ncomp
! !       n = 1
! !       x_int = 0.0
! !       y_int = 0.0
! !       y_intd = 0.0
! !       DO jj=ii-fa_disp(icomp),ii+fa_disp(icomp)
! ! !IF ( zp(jj+1) > (zl+zmin) ) THEN
! !         IF (zp(ii) - zp(jj+1) .LT. ab_disp(icomp) .AND. zp(ii) - zp(jj) &
! ! &           .GE. ab_disp(icomp)) THEN
! ! !the position of j+1 is already within i-d/2 while j is still outside this range in this case, the integration steplength (dz) is
! ! ! just the distance, which j+1 overlaps with i-d/2 and what is integrated is the interpolated value of the integrand
! ! !here always n=1!
! !           IF (n .NE. 1) THEN
! !             GOTO 100
! !           ELSE
! ! !distance between grid points j and i
! !             zz = zp(jj) - zp(ii)
! ! !the part of the intervall between zp(j) and zp(j+1) which is already within i-d/2
! !             dz = zp(jj+1) - (zp(ii)-ab_disp(icomp))
! ! !if(dz < epsilon(dz)) dz = epsilon(dz) !bei unguenstiger Kombination von sig und ngrid kann dz unter Machinengenauigkeit epsilon 
! ! !liegen, dann ist x(2) = x(1) + dz = x(1) -> das fuehrt zu Abbruch in Spline Interpolation           
! ! !array containing x-values for spline integration 
! !             x_int(n) = 0.0
! !             y_intd(n) = 0.0
! !             y_int(n) = 0.0
! !           END IF
! !         ELSE IF (zp(jj) .GT. zp(ii) - ab_disp(icomp) .AND. zp(jj) .LE. &
! ! &           zp(ii) + ab_disp(icomp)) THEN
! ! !grid point j  within i+-d2 
! !           n = n + 1
! ! !first time in this If condition, dz is stil the old value from above!
! !           x_int(n) = x_int(n-1) + dz
! !           zz = zp(jj) - zp(ii)
! !           dz = dzp
! !           y_intd(n) = (ab_disp(icomp)*ab_disp(icomp)-zz*zz)*my_dispd(jj&
! ! &           , icomp)
! !           y_int(n) = my_disp(jj, icomp)*(ab_disp(icomp)*ab_disp(icomp)-&
! ! &           zz*zz)
! ! !rhop_sum = sum(rhop(1:ncomp,jj))
! ! !y_int(n) = rhop_sum * df_disp_drk(jj,icomp) * (ab_disp(icomp)*ab_disp(icomp) - zz*zz) 
! !           IF (zp(jj) .LT. zp(ii) + ab_disp(icomp) .AND. zp(jj+1) .GE. zp&
! ! &             (ii) + ab_disp(icomp)) THEN
! ! !zp(j) is still within zp(i)+d2 but zp(j+1) is already outside zp(i)+d2
! !             dz = zp(ii) + ab_disp(icomp) - zp(jj)
! ! !If(dz <= epsilon(dz)) exit   !wie oben, kann auch hier bei ungluecklicher Wahl von sig und ngrid dz < eps werden und somit x(n) 
! ! != x(n-1) -> Abbruch in Spline interpolation. Dann einfach ngrid aendern!
! !             zz = zp(jj) - zp(ii)
! !             n = n + 1
! !             x_int(n) = x_int(n-1) + dz
! !             y_intd(n) = 0.0
! !             y_int(n) = 0.0
! !           END IF
! !         END IF
! !       END DO
! !       xlo = x_int(1)
! !       xhi = x_int(n)
! !       CALL SPLINE_D(x_int(1:n), y_int(1:n), y_intd(1:n), n, 1.e30, 1.e30&
! ! &             , y2(1:n), y2d(1:n))
! !       CALL SPLINT_INTEGRAL_D(x_int(1:n), y_int(1:n), y_intd(1:n), y2(1:n&
! ! &                      ), y2d(1:n), n, xlo, xhi, int2, int2d)
! !       df_drho_dispd(icomp) = 0.75*int2d/ab_disp(icomp)**3
! !       df_drho_disp(icomp) = int2*0.75/ab_disp(icomp)**3
! !     END DO
! !     GOTO 110
! !  100 STOP
! !  110 CONTINUE
  END SUBROUTINE disp_dfdrho_wda_d

END MODULE mod_dft_disp_wda_d