c************************************************************ PROGRAM test_coef_diff c test pour les coefficients de diffusion c CESAM95 c Auteur: P. Morel, Departement J.D. Cassini, O.C.A. c--------------------------------------------------------------------- USE mod_chim, ONLY : diffm, difft USE mod_donnees, ONLY : abe7, ac12, ac13, ah, ah2, ahe3, ahe4, 1 ali7, an14, an15, ao16, ao17, aradia, clight, g, ini_ctes_phys, lit_nl, 2 lsol, msol, nom_diffm, nom_difft, nom_fich2, nom_osc, 3 pi, print_ctes, rsol USE mod_exploit, ONLY : read_osc USE mod_kind USE mod_numerique, ONLY : pause USE mod_sta_chim, ONLY : ihe4, iw, nbelem, nchim, nom_elem, 1 nucleo, zi IMPLICIT NONE INTEGER, PARAMETER :: pnosc=800, pnelem=15 REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: dd REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: var, y, d, d0, dv REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: r, p, t, ro, l, m, 1 xchim, v, v0, glob, kap, gradrad, gradad, gamma1, cp, delta REAL (kind=dp), DIMENSION(3) :: m_zc=0.d0, r_zc=0.d0 REAL (kind=dp), parameter :: dx=1.d-5, unpdx=1.d0+dx REAL (kind=dp) :: mstar, stor, stor0, dstor, cte8, 1 drox=0.d0, dkapx=0.d0, dgradradx=0.d0, dgradadx=0.d0 REAL (kind=sp), ALLOCATABLE, DIMENSION(:) :: xdes, dh, dhe4, dli7, 1 vh, vhe4, vli7 REAL (kind=sp) :: xmax, xmin, ymax, ymin INTEGER :: iconst, ivar, i, itot, j, k, ll, ides, ili7=-1, lim=0 LOGICAL :: ok, dessin, ps, melange, deriv CHARACTER (len=4), ALLOCATABLE, DIMENSION(:) :: nom_elem_lu CHARACTER (len=9) :: device, titre='diffusion' CHARACTER (len=80) :: modele, nom_fich c--------------------------------------------------------------------- c itot: nombre de points c iconst: nombre de constantes global c ivar: nombre de variables du tableau var sans les elements chimiques c modele: nom du modele c nom_elem: nom des elements utilises c nbelem: nombre d'elements utilises c nchim: nombre d'elements chimiques utilises c glob: variables globales c glob(1)=mstar*msol c glob(2)=rstar*rsol c glob(3)=ltot*lsol c glob(4)=z c glob(5)=x0 c glob(6)=alpha c glob(7)=9./4. c glob(8)=1./162. c glob(9)=1. c glob(10)=1. c glob(11)=d2p c glob(12)=d2ro c glob(13)=age c glob(14)=vsal c glob(15)=0. c var: variables c var(1,i)=r*rsol c var(2,i)=log(m) -1.d38 au centre c var(3,i)=t c var(4,i)=p c var(5,i)=ro c var(6,i)=gradient c var(7,i)=l c var(8,i)=kap c var(9,i)=energie thermo+gravifique c var(10,i)=grand Gamma1 c var(11,i)=gradient adiabatique c var(12,i)=delta c var(13,i)=cp c var(14,i)=1 / (mu elec.) c var(15,i)=vaissala, 0 au centre c var(16,i)=w(i) c var(17,i)=d ln kappa / d ln T c var(18,i)=d ln kappa / d ln ro c var(19,i)=d epsilon(nuc) / d ln T c var(20,i)=d epsilon(nuc) / d ln ro c var(21,i)=Ptot/Pgas c var(22,i)=gradrad c var(22+j,i)=xchim(j), j=1,nbelem 2000 FORMAT(8es10.3) 2001 FORMAT(10es8.1) c le test utilise le fichier de donnees nom_fich2 nom_fich2='test' !sans extension .don omise c qui a créé le fichier le fichier d'oscillations nom_fich=nom_fich2 WRITE(*,1)nom_fich 1 FORMAT('on lit le fichier d''oscillation:',a80) CALL lit_nl ; CALL ini_ctes_phys CALL PRINT_ctes(6) cte8=3.d0/16.d0/pi/aradia/clight/g*lsol/msol ALLOCATE(glob(15),var(22+pnelem,pnosc),nom_elem_lu(pnelem)) CALL read_osc(nom_fich,itot,iconst,ivar,modele,glob,var, 1 nbelem,nchim,iw,nom_elem_lu,ok) mstar=glob(1)/msol ALLOCATE(r(itot),p(itot),t(itot),ro(itot),l(itot),m(itot), 1 kap(itot),y(nbelem,itot),gradrad(itot),gradad(itot), 2 xdes(itot),dh(itot),dhe4(itot),dli7(itot),vh(itot), 3 vhe4(itot),vli7(itot),gamma1(itot),cp(itot),delta(itot), 4 nom_elem(nbelem),nucleo(nbelem),zi(nbelem)) m=0.d0 IF(.NOT.ok)stop DO i=1,itot IF(var(2,i) > -1.d30)m(i)=exp(var(2,i)) p(i)=var(4,i) t(i)=var(3,i) ro(i)=var(5,i) l(i)=var(7,i)/lsol r(i)=var(1,i)/rsol kap(i)=var(8,i) gradrad(i)=var(22,i) gradad(i)=var(11,i) gamma1(i)=var(10,i) cp(i)=var(13,i) delta(i)=var(12,i) DO j=1,nbelem y(j,i)=max(var(ivar+j,i),1.d-30) ENDDO ENDDO nom_elem=nom_elem_lu(1:nbelem) DEALLOCATE(var,glob,nom_elem_lu) DO i=1,nbelem IF(nom_elem(i) == ' H1 ')THEN nucleo(i)=ah zi(i)=1 ELSEIF(nom_elem(i) == ' H2 ')THEN nucleo(i)=ah2 zi(i)=1 ELSEIF(nom_elem(i) == 'He3 ')THEN nucleo(i)=ahe3 zi(i)=2 ELSEIF(nom_elem(i) == 'He4 ')THEN nucleo(i)=ahe4 zi(i)=2 ihe4=i ELSEIF(nom_elem(i) == 'Li7 ')THEN nucleo(i)=ali7 zi(i)=3 ili7=i ELSEIF(nom_elem(i) == 'Be7 ')THEN nucleo(i)=abe7 zi(i)=4 ELSEIF(nom_elem(i) == 'C12 ')THEN nucleo(i)=ac12 zi(i)=6 ELSEIF(nom_elem(i) == 'C13 ')THEN nucleo(i)=ac13 zi(i)=6 ELSEIF(nom_elem(i) == 'N14 ')THEN nucleo(i)=an14 zi(i)=7 ELSEIF(nom_elem(i) == 'N15 ')THEN nucleo(i)=an15 zi(i)=7 ELSEIF(nom_elem(i) == 'O16 ')THEN nucleo(i)=ao16 zi(i)=8 ELSEIF(nom_elem(i) == 'O17 ')THEN nucleo(i)=ao17 zi(i)=8 ELSEIF(nom_elem(i) == 'Si28')THEN nucleo(i)=28 zi(i)=14 ELSEIF(nom_elem(i) == ' Ex ')THEN nucleo(i)=28 zi(i)=13 ELSE PRINT*,'erreur de nom d''element: ',nom_elem(i) ; STOP ENDIF ENDDO ALLOCATE(xchim(nbelem),v0(nbelem),dv(nbelem,nbelem),d0(nbelem,nbelem), 1 dd(nbelem,nbelem,nbelem),v(nbelem),d(nbelem,nbelem)) ides=0 c le test b1: DO ll=1,itot c IF(t(ll) < 3.872E+06)CYCLE b1 c IF(t(ll) > 3.933E+06)CYCLE b1 c IF(t(ll) < 1.d5)CYCLE b1 c IF(t(ll) < 1.d6)CYCLE b1 melange=gradrad(ll) > gradad(ll) IF(melange)CYCLE b1 IF(m(ll)*r(ll) <= 0.d0)CYCLE b1 xchim(:)=y(:,ll)/nucleo(:) !comp.chim. par mole WRITE(*,2000)t(ll),r(ll) ; WRITE(*,2001)xchim v0=0.d0 ; d0=0.d0 ; dv=0.d0 ; dd=0.d0 CALL diffm(p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,kap(ll),dkapx, 1 gradrad(ll),dgradradx,gradad(ll),dgradadx,mstar,xchim,d0,dd,v0,dv) PRINT*,'V ',nom_diffm WRITE(*,2001)(v0(i),i=1,nbelem) PRINT*,'D ',nom_diffm DO i=1,nbelem WRITE(*,2001)(d0(i,j),j=1,nbelem) ENDDO v0=0.d0 ; d0=0.d0 ; dv=0.d0 ; dd=0.d0 CALL difft(melange,p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox, 1 kap(ll),dkapx,gradrad(ll),dgradradx,gradad(ll), 2 dgradadx,gamma1(ll),cp(ll),delta(ll),xchim,d0,dd,v0,dv) PRINT*,'melange=',melange ; PRINT*,'V ',nom_difft WRITE(*,2001)(v0(i),i=1,nbelem) ; PRINT*,'D ',nom_difft DO i=1,nbelem WRITE(*,2001)(d0(i,j),j=1,nbelem) ENDDO CALL pause('dans coeff_diff') c DO k=1,nbelem c PRINT*,'dD / Xk, k=',k c DO i=1,nbelem c WRITE(*,2000)(dd(i,j,k),j=1,nbelem) c ENDDO c ENDDO c pause deriv=.FALSE. c deriv=.TRUE. IF(deriv)THEN WRITE(*,2000)r(ll),t(ll) DO i=1,nchim v=0.d0 ; d=0.d0 ; dv=0.d0 ; dd=0.d0 stor0=xchim(i) ; stor=stor0*unpdx ; dstor=stor-stor0 xchim(i)=stor c CALL diffm(p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,kap(ll),dkapx, c 1 gradrad(ll),dgradradx,gradad(ll),dgradadx,mstar,xchim,d,dd,v,dv) CALL difft(melange,p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox, 1 kap(ll),dkapx,gradrad(ll),dgradradx,gradad(ll), 2 dgradadx,gamma1(ll),cp(ll),delta(ll),xchim,d,dd,v,dv) xchim(i)=stor0 WRITE(*,"('derivee / ',a4)")nom_elem(i) c WRITE(*,2000)v0 ; WRITE(*,2000)v WRITE(*,2000)((v(j)-v0(j))/dstor,dv(j,i),j=1,nchim) ENDDO CALL pause('ici') ENDIF ides=ides+1 ; xdes(ides)=r(ll) dh(ides)=d0(1,1) ; dhe4(ides)=d0(ihe4,ihe4) IF(ili7 > 0)dli7(ides)=d0(ili7,ili7) vh(ides)=v0(1) ; vhe4(ides)=abs(v0(ihe4)) IF(ili7 > 0)vli7(ides)=abs(v0(ili7)) ENDDO b1 dessin=.false. ; dessin=.true. ; ps=.false. IF(dessin)THEN xmax=maxval(xdes)*1.05 ; xmin=0. IF(ili7 > 0)THEN ymax=MAX(MAXVAL(dh),MAXVAL(dhe4),MAXVAL(dli7))*1.05 ELSE ymax=MAX(MAXVAL(dh),MAXVAL(dhe4))*1.05 ENDIF ymin=-ymax/10. IF(.NOT.ps)THEN device='/xw' ELSE device='/ps' ENDIF c WRITE(*,*)'device ? /xw, /PS, /CPS' c read(5,'(a)')device CALL pgbegin(0,device,1,2) CALL pgenv(xmin,xmax,ymin,ymax,0,0) CALL pglabel('R/Rsol','Diff',titre) CALL pgsls(1) CALL pgline(ides,xdes,dh) CALL pgsls(2) CALL pgline(ides,xdes,dhe4) CALL pgsls(3) IF(ili7 > 0)CALL pgline(ides,xdes,dli7) CALL pgsls(1) IF(ili7 > 0)THEN ymax=max(maxval(vh),maxval(vhe4),maxval(vli7))*1.05 ELSE ymax=max(maxval(vh),maxval(vhe4))*1.05 ENDIF ymin=-ymax/10. CALL pgenv(xmin,xmax,ymin,ymax,0,0) CALL pglabel('R/Rsol','Vdiff',titre) CALL pgsls(1) CALL pgline(ides,xdes,vh) CALL pgsls(2) CALL pgline(ides,xdes,vhe4) CALL pgsls(3) IF(ili7 > 0)CALL pgline(ides,xdes,vli7) CALL pgend ENDIF STOP END PROGRAM test_coef_diff