!-----------------------------------------------------------------------
  ! Copyright (C) 2010-2015 Henry Lambert, Feliciano Giustino
  ! This file is distributed under the terms of the GNU General Public         
  ! License. See the file `LICENSE' in the root directory of the               
  ! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
  !-----------------------------------------------------------------------
SUBROUTINE construct_w(scrcoul_g, scrcoul_pade_g, w_ryd, xq_ibk)
  USE kinds,         ONLY : DP
  USE constants,     ONLY : e2, fpi, RYTOEV, tpi, eps8, pi
  USE control_gw,    ONLY : lgamma, eta, godbyneeds, padecont, modielec, trunc_2d, do_imag
  USE freq_gw,       ONLY : fpol, fiu, nfs, nfsmax, &
                            nwcoul, nwgreen, nwalloc, nwsigma, wtmp, wcoul, &
                            wgreen, wsigma, wsigmamin, wsigmamax, &
                            deltaw, wcoulmax
  USE gwsigma,       ONLY : sigma_c_st
  USE gvect,         ONLY : g, ngm, nl
  USE disp,          ONLY : nqs, nq1, nq2, nq3, wq, x_q, xk_kpoints
  USE cell_base,     ONLY : tpiba2, tpiba, omega, alat, at
  USE mp_global,     ONLY : mp_global_end
  USE mp_world,      ONLY : nproc, mpime

  IMPLICIT NONE

  COMPLEX(DP) :: scrcoul_pade_g (sigma_c_st%ngmt, sigma_c_st%ngmt)
  COMPLEX(DP) :: z(nfs), u(nfs), a(nfs)
  COMPLEX(DP)  :: scrcoul_g    (sigma_c_st%ngmt, sigma_c_st%ngmt, nfs) 

  REAL(DP) :: qg2, qg, qxy, qz
  REAL(DP) :: w_ryd
  REAL(DP) :: rcut, spal, zcut
  REAL(DP) :: xq_ibk(3), xq_ibz(3)

  INTEGER :: ig, igp, irr, icounter, ir, irp
  INTEGER  :: iwim, iw, ikq

  LOGICAL             :: pade_catch
  LOGICAL             :: found_q
  LOGICAL             :: limq, inv_q, found

  zcut = 0.50d0*sqrt(at(1,3)**2 + at(2,3)**2 + at(3,3)**2)*alat
  rcut = (float(3)/float(4)/pi*omega*float(nq1*nq2*nq3))**(float(1)/float(3))
  scrcoul_pade_g(:,:) = (0.0d0, 0.0d0)
  IF(.NOT.modielec) THEN
    DO ig = 1, sigma_c_st%ngmt
       DO igp = 1, sigma_c_st%ngmt
          DO iwim = 1, nfs
              z(iwim) = fiu(iwim)
              a(iwim) = scrcoul_g (ig,igp,iwim)
          ENDDO
          IF (padecont) THEN
              call pade_eval ( nfs, z, a, dcmplx(0.0d0, w_ryd), scrcoul_pade_g (ig,igp))
          ELSE IF (godbyneeds .and. do_imag) THEN
              scrcoul_pade_g(ig,igp) = -a(2)/(w_ryd**2+(a(1))**2)
          ELSE IF (godbyneeds .and. .not. do_imag) THEN
              scrcoul_pade_g(ig,igp) = a(2)/(dcmplx(w_ryd**2,0.0d0)-(a(1)-(0.0d0,1.0d0)*eta)**2)
          ELSE
               WRITE(6,'("No screening model chosen!")')
               STOP
               CALL mp_global_end()
          ENDIF
       ENDDO
    ENDDO
  ELSE IF (modielec) THEN
!Generate Mod_dielectric
       CALL mod_diel(scrcoul_pade_g(1,1), xq_ibk, w_ryd, 3)
!Apply COULOMB
       IF(.not.trunc_2d) THEN
          DO ig =1, sigma_c_st%ngmt 
            qg2 = (g(1,ig) + xq_ibk(1))**2 + (g(2,ig) + xq_ibk(2))**2 + (g(3,ig)+xq_ibk(3))**2
            limq = (qg2.lt.eps8) 
            IF(.not.limq) THEN
               scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx(e2*fpi/(tpiba2*qg2), 0.0d0)
            ENDIF
            qg = sqrt(qg2)
            spal = 1.0d0 - cos(rcut*sqrt(tpiba2)*qg)
!Normal case using truncated coulomb potential.
            IF(.not.limq) THEN
               scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx(spal, 0.0d0)
            ELSE
               scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx((fpi*e2*(rcut**2))/2.0d0, 0.0d0)
            ENDIF
          ENDDO
       ELSE
         DO ig = 1, sigma_c_st%ngmt
            qg2 = (g(1,ig) + xq_ibk(1))**2 + (g(2,ig) + xq_ibk(2))**2 + (g(3,ig)+xq_ibk(3))**2
            qxy  = sqrt((g(1,ig) + xq_ibk(1))**2 + (g(2,ig) + xq_ibk(2))**2)
            qz   = sqrt((g(3,ig)+xq_ibk(3))**2)
            limq = (qg2.lt.eps8) 
            !write(1000+mpime, *) scrcoul_pade_g(ig,ig)
            IF (qxy.gt.eps8) then
              spal = 1.0d0 + EXP(-tpiba*qxy*zcut)*((qz/qxy)*sin(tpiba*qz*zcut) - cos(tpiba*qz*zcut))
              scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx((e2*fpi/(tpiba2*qg2))*spal, 0.0d0)
            ELSE IF (qxy.lt.eps8.and.qz.gt.eps8) then
              spal = 1.0d0 - cos(tpiba*qz*zcut) - tpiba*qz*zcut*sin(tpiba*qz*zcut)
              scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx((e2*fpi/(tpiba2*qg2))*spal, 0.0d0)
            ELSE  
              scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx(rcut,0.0d0)
            ENDIF
            !write(1000+mpime, *) scrcoul_pade_g(ig,ig)
         ENDDO
     ENDIF
  ENDIF
END SUBROUTINE construct_w