c*************************************************************

	SUBROUTINE diffw_mpz(dlnro,grand_k,nu,nu_cin,n2mu,n2t,r,ro,y,
	1 deff,dh,dv,dhv)

c routine private du module mod_evol

c calcul des coefficients de diffusion turbulente pour la rotation
c Dh, Dv, Deff, suivant Matis, Palacio & Zahn A&A 2004
c en optional dérivées de Dh et Dv / Omega, U, dOmega, dU

c entrées :
c	dlnro : d ln ro / dnu
c	grand_k : K
c	nu : m^2/3
c	nu_cin : viscosité cinématique
c	y(:,0:1) : variables propres à la rotation et dérivées / nu

c sorties :
c	Deff, Dh, Dv : notations évidentes

c	Auteur: P. Morel, Département Cassiopée, O.C.A., CESAM2k

c---------------------------------------------------------------------

	USE mod_donnees, ONLY : msol, nom_rot, nrot, pi, rsol 
	USE mod_kind

	IMPLICIT NONE

	REAL (kind=dp), INTENT(in), DIMENSION(2,0:1) :: y	
	REAL (kind=dp), INTENT(in) :: dlnro, grand_k, nu, nu_cin, n2mu, n2t,
	1 r, ro
	REAL (kind=dp), INTENT(out), DIMENSION(2,2,0:1) :: dhv
	REAL (kind=dp), INTENT(out) :: deff, dh, dv
	
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: yd	
	REAL (kind=dp), DIMENSION(2,0:1) :: dgskdh
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: hve, hved
	REAL (kind=dp), PARAMETER :: betta=1.5d-6, dx=1.d-5, ric=1.d0/6.d0,
	1 unpdx=1.d0+dx
	REAL (kind=dp), SAVE :: cte1_0, cte3_0, cte5_0, cte6_0, cts1_0, cts2_0
	REAL (kind=dp) :: cte1, cte2, cte3, cte5, den, dh2, dstor, gksdh,
	1 gksdhd, nu12, stor, stor0

	INTEGER :: i, id, j
	
	LOGICAL :: deriv=.FALSE.	
	LOGICAL, SAVE :: init=.TRUE.
			
c------------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

	IF(init)THEN
	 init=.FALSE.
	 cte1_0=8.d0*pi*rsol**6*betta/9.d0/msol
	 cte3_0=4.d0*pi*rsol**6*betta/3.d0/msol	 
	 cte5_0=64.d0*pi**2*ric*rsol**6/9.d0/msol**2
	 cte6_0=rsol**2/30.d0	 
	 cts1_0=8.d0*pi*rsol**6*betta/9.d0/msol
	 cts2_0=rsol**3*betta/3.d0
	ENDIF
	
c initialisations
	dhv=0.d0
	nu12=SQRT(ABS(nu))		
	cte1=cte1_0*r**6*ro/nu12	
	cte2=(cts1_0*r**3*ro/nu12*dlnro-cts2_0)*r**3
	cte3=cte3_0*r**6*ro/nu12

c Dh
	dh2=cte1*y(1,0)*y(2,1)+cte2*y(1,0)*y(2,0)-cte3*y(1,1)*y(2,0)
	IF(SQRT(ABS(dh2)) <= nu_cin)THEN
	 dh=nu_cin ; dv=nu_cin ; deff=nu_cin ; RETURN
	ELSE
	 dhv(1,1,0)=cte1*y(2,1)+cte2*y(2,0)	!d dh2 / d Omega
	 dhv(1,1,1)=-cte3*y(2,0)			!d dh2 / d dOmega	
	 dhv(1,2,0)=cte2*y(1,0)-cte3*y(1,1)	!d dh2 / d U
	 dhv(1,2,1)=cte1*y(1,0)			!d dh2 / d dU	 	
	 IF(dh2 < 0.d0)THEN	   
	  dh2=-dh2 ; dhv(1,:,:)=-dhv(1,:,:)
	 ENDIF
	 dh=SQRT(dh2) ; dhv(1,:,:)=dhv(1,:,:)/2.d0/dh

c Dv  	  	
	 gksdh=1.d0+grand_k/dh
	 dgskdh(:,:)=-grand_k*dhv(1,:,:)/dh2	 
	 cte5=cte5_0*r**6*ro**2/nu	 
	 den=n2t+n2mu*gksdh 
	 dv=cte5*dh*gksdh*y(1,1)**2/den
	 IF(dv > nu_cin)THEN
	  dhv(2,1,0)=dv*(dhv(1,1,0)/dh
	1 +dgskdh(1,0)*(1.d0/gksdh-n2mu/den))		!d dv / d Omega
	  dhv(2,2,0)=dv*(dhv(1,2,0)/dh
	1 +dgskdh(2,0)*(1.d0/gksdh-n2mu/den))		!d dv / d U
	  dhv(2,1,1)=dv*(2.d0/y(1,1)+dhv(1,1,1)/dh+
	1 dgskdh(1,1)*(1.d0/gksdh-n2mu/den))	 	!d dv / d dOmega
	  dhv(2,2,1)=dv*(dhv(1,2,1)/dh
	1 +dgskdh(2,1)*(1.d0/gksdh-n2mu/den))		!d dv / d dU
	  dv=dv+nu_cin
	 ELSE
	  dv=nu_cin
	 ENDIF

c Deff
	 deff=cte6_0*(r*y(2,0))**2/dh+dv

	ENDIF		!dh <= nu_cin
	
c test de dérivées
c	 deriv=nu >= 6.E-01 .AND. nu <= 6.3E-01
c	 deriv=nu >= 6.2E-01 .AND. nu <= 6.25E-01
c	 deriv=dh > 1.d2
	 IF(deriv)THEN
	  PRINT*,'nu' ; WRITE(*,2000)nu
	  PRINT*,'y' ; WRITE(*,2000)y(:,0)
	  PRINT*,'dy' ; WRITE(*,2000)y(:,1)
	  PRINT*,'grand_k,gksdh,n2t,n2mu,den,Dh2,Dh,Dv'
	  WRITE(*,2000)grand_k,gksdh,n2t,n2mu,den,dh2,dh,dv
c	  PRINT*,'dv,dhv,dgskdh*,dgskdh,1.d0/gksdh,n2mu/den'
c	  WRITE(*,2000)dv,dhv(1,1,0)/dh,dgskdh(1,0)*(1.d0/gksdh-n2mu/den),
c	1 dgskdh(1,0),1.d0/gksdh,n2mu/den	  	 	  
	  ALLOCATE(hve(2),hved(2),yd(2,0:1))
	  hve(1)=dh ; hve(2)=dv	 
	  DO id=0,1
	   B1: DO i=1,2
	    yd=y
	    stor0=yd(i,id) ; stor=stor0*unpdx !; IF(ABS(stor) < dx)stor=dx
	    yd(i,id)=stor ; dstor=stor-stor0 
	    hved(1)=SQRT(ABS(cte1*yd(1,0)*yd(2,1)+cte2*yd(1,0)*yd(2,0)
	1   -cte3*yd(1,1)*yd(2,0)))
	    gksdhd=1.d0+grand_k/hved(1)
	    den=n2t+n2mu*gksdhd	   
	    hved(2)=cte5*hved(1)*gksdhd*yd(1,1)**2/den
	    IF(hved(2) > nu_cin)hved(2)=hved(2)+nu_cin
	    IF(id == 0)THEN
	     PRINT*,'diffw_mpz pour diff, dérivées ',nom_rot(i)
	    ELSE
	     PRINT*,'diffw_mpz pour diff, dérivées d',nom_rot(i)
	    ENDIF
c	    WRITE(*,*)'Dh,Dv'
c	    WRITE(*,2000)hved(1),dh,hved(2),dv
	    WRITE(*,2000)(dhv(j,i,id),(hved(j)-hve(j))/dstor,j=1,2),
	1   dgskdh(i,id),(gksdhd-gksdh)/dstor
	    PAUSE
	   ENDDO B1
	  ENDDO
	  DEALLOCATE(hved,yd) ; deriv=.FALSE.
	 ENDIF 	!deriv	

	RETURN
	
	END SUBROUTINE diffw_mpz