c****************************************************************** SUBROUTINE diffw_p03(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 Palacios & al A&A 2003 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 :: ch=1.d0, dx=1.d-5, ric=2.d0/5.d0, 1 unpdx=1.d0+dx REAL (kind=dp), SAVE :: cte1_0, cte2_0, cte3_0, cte5_0, cte6_0, 1 cts1_0, cts2_0 REAL (kind=dp) :: cte1, cte2, cte3, cte5, cts1, cts2, den, dstor, 1 gksdh, 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. cts1_0=8.d0*pi*rsol**4/9.d0/msol cts2_0=rsol/3.d0 cte1_0=1.d0/ch cte2_0=8.d0*pi*rsol**4/9.d0/msol cte3_0=4.d0*pi*rsol**4/3.d0/msol cte5_0=64.d0*pi**2*ric*rsol**6/9.d0/msol**2 cte6_0=rsol**2/30.d0 ENDIF c initialisations dhv=0.d0 nu12=SQRT(ABS(nu)) cts1=cts1_0*r**4/nu12 cts2=cts2_0*r cte1=(cts1-cts2)*r**4*ro*dlnro/nu12*cte1_0 cte2=cte2_0*r**4*ro/nu12*cte1_0 cte3=cte3_0*r**4*ro/nu12*cte1_0 c Dh dh=cte1*y(2,0)+cte2*y(2,1)-cte3*y(2,0)*y(1,1)/y(1,0) dh=ABS(dh) IF(dh <= nu_cin)THEN dh=nu_cin ; dv=nu_cin ; deff=nu_cin ; RETURN ELSE dhv(1,1,0)=+cte3*y(2,0)*y(1,1)/y(1,0)**2 !d dh / d Omega dhv(1,1,1)=-cte3*y(2,0)/y(1,0) !d dh / d dOmega dhv(1,2,0)=cte1-cte3*y(1,1)/y(1,0) !d dh / d U dhv(1,2,1)=cte2 !d dh / d dU c Dv gksdh=1.d0+grand_k/dh dgskdh(:,:)=-grand_k*dhv(1,:,:)/dh**2 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 (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,Dh,Dv' WRITE(*,2000)grand_k,gksdh,n2t,n2mu,den,dh,dv ALLOCATE(hve(2),hved(2),yd(2,0:1)) yd=y ; hve(1)=dh ; hve(2)=dv DO id=0,1 B1: DO i=1,2 stor0=yd(i,id) ; stor=stor0*unpdx !; IF(ABS(stor) < dx)stor=dx yd(i,id)=stor ; dstor=stor-stor0 hved(1)=cte1*yd(2,0)+cte2*yd(2,1)-cte3*yd(2,0)*yd(1,1)/yd(1,0) gksdhd=1.d0+grand_k/hved(1) den=ABS(n2t)+n2mu*gksdhd hved(2)=cte5*hved(1)*gksdhd*yd(1,1)**2/den 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(*,2000)hved(1),hve(1),yd(1:2,0),yd(1:2,1) WRITE(*,2000)(dhv(j,i,id),(hved(j)-hve(j))/dstor,j=1,2), 1 dgskdh(i,id),(gksdhd-gksdh)/dstor yd(i,id)=stor0 ; PAUSE ENDDO B1 ENDDO DEALLOCATE(hved,yd) ; deriv=.FALSE. ENDIF !deriv RETURN END SUBROUTINE diffw_p03