c---------------------------------------------------------------------- SUBROUTINE grad2_bp(grav,lum,ray,t,nel,ychim,ioni,g_rad,dg_rad) ! Version 2.2 (16/02/2007) ! Cette SUBROUTINE calcule l'acceleration radiative et sa derivee ! par rapport a l'abondance. Elle DOit etre appelee pour chaque couche, ! et a chaque fois que la dIFfusion est calculee. gr_popion3 ! Cette SUBROUTINE fait partie du module grcesam2. ! Version 2.1 (14/06/2004) ! Auteur: ! Georges ALECIAN ! LUTH - UMR 8102, Observatoire de MeuDOn ! F-92195 MEUDON CEDEX, FRANCE ! Tel: 01 45 07 74 20, + 33 1 45 07 74 20 c Adaptation a CESAM2k par B. Pichon & P. Morel, OCA c--------------------------------------------------------------- USE mod_donnees, ONLY: nchim, nucleo USE mod_kind IMPLICIT NONE REAL(DP), DIMENSION(0:,:), INTENT(in) :: ioni REAL(DP), DIMENSION(:), INTENT(in) :: ychim REAL(DP), INTENT(in) :: grav, lum, nel, ray, t REAL(DP), DIMENSION(:), INTENT(out) :: g_rad REAL(DP), DIMENSION(:,:), INTENT(out) :: dg_rad ! pour le g_cont et gr_popion3 REAL(DP), DIMENSION(0:pzi) :: fparti, theta REAL(DP), DIMENSION(30,0:pzi) :: ff REAL(DP), DIMENSION(0:pzi,pnchim) :: popi REAL(DP), DIMENSION(nchim) :: khi, N_chim REAL(DP), DIMENSION(0:pzi,nchim) :: C_ion, pond_ion ! mes variables locales REAL(DP), PARAMETER :: pondrare = 1.5d0 REAL(DP) :: b, bco, Ck_s, CX, dgc_kj, dgr_kj, f, fact_temp,gc_kj,gr_kj, 1 g_grav, pondt, q, x INTEGER :: iso, j, k !************************** 2000 FORMAT(8es10.3) ! securites popi=0.d0 !tableau g_rad = 0.d0 dg_rad = 0.d0 IF(ychim(1)*t*ray < 1.d-50) RETURN khi = 1.d0 !tableau gc_kj = 0.d0 dgc_kj = 0.d0 ! module du vecteur gravite g_grav = grav CCCC CCCC ! On considere que les isotopes ont leurs transitions atomiques ! presque aux memes frequences (bf comme bb). ! Cela implique que l'effet de saturation DOit etre presque le meme ! pour tous les isotopes d'un meme element. Cette saturation est reglee par ! la concentration C_ion(k,j) (en pratique: abondance en nb par rapport a H). ! On veut DOnc que C_ion(k,j) soit le meme pour tous les isotopes d'un element ! et soit DOnne par la concentration totale pour l'element. ! Infos: Somme[ychim(j)*nucleo(j)]=1 DO j=1,nchim N_chim(j)=0. DO iso=1,nchim N_chim(j)=N_chim(j)+ychim(iso)*isotops(j,iso) ENDDO ENDDO DO j=1,nchim DO k=0,NINT(zi(j)) C_ion(k,j)=ioni(k,j)*N_chim(j)/N_chim(1) ENDDO khi(j)=(N_chim(j)/N_chim(1)) / C_sol(j) ENDDO ! coefficients de ponderation pour le g_rad total DO j=1,nchim pondt=0.d0 DO k=0,NINT(zi(j)) pondt = pondt + ioni(k,j)*(1.d0+rar_flag(k,j)*(pondrare-1.d0)) ENDDO DO k=0,NINT(zi(j)) pond_ion(k,j) = ioni(k,j)*(1.d0+rar_flag(k,j)*(pondrare-1.d0))/pondt ENDDO ENDDO DO j=1,nchim IF(id_ppxi(j) == 1) THEN call gr_parti3(t,nel,fparti,ff,j) call gr_popion3(t,nel,fparti,ff,popi,j) DO k=0,NINT(zi(j)) C_ion(k,j)=popi(k,j)*N_chim(j)/N_chim(1) ENDDO pondt=0.d0 DO k=0,NINT(zi(j)) pondt = pondt + popi(k,j)*(1. + rar_flag(k,j)*(pondrare-1.) ) ENDDO DO k=0,NINT(zi(j)) pond_ion(k,j) = popi(k,j)*(1. + rar_flag(k,j)*(pondrare-1.) )/pondt ENDDO ENDIF ENDDO ! Fin substitution !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! constantes: ! pi**2 * kbol**3 * echarg**2 ! 5.57E-05 = --------------------------------- ! 2.* clight**4 * hpl**2 * me * amu ! ! ! me * mp* clight ! 9.83E-23 = ----------------- ! 2. * pi ! ! CX : fraction de masse de l'hydrogene !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF(tgrad == 0) THEN DO j=1,nchim IF(id_ppxi(j) == 1) THEN g_rad(j) = 0.d0 ELSE IF(igrad == 0) THEN g_rad(j) = 0.d0 ELSEIF(igrad == 1) THEN g_rad(j) = g_grav *(1-pond_ion(NINT(zi(j)),j)) ! correction noyau nu, ! le noyau nu ne peut pas avoir g_rad = -g, on suppose que la photoionisation est 0 ENDIF ENDIF ENDDO ELSEIF(tgrad == 1) THEN DO j=1,nchim IF(id_ppxi(j) == 1) THEN g_rad(j) = g_grav ELSE IF(igrad == 0) THEN g_rad(j) = 0.d0 ELSEIF(igrad == 1) THEN g_rad(j) = g_grav*(1-pond_ion(NINT(zi(j)),j)) !corr. noyau nu, ! le noyau nu ne peut pas avoir g_rad = -g, on suppose photoionisation = 0 ENDIF ENDIF ENDDO ELSEIF(tgrad == 2) THEN DO j=1,nchim IF(id_ppxi(j) == 1) THEN call gr_parti3(t,nel,fparti,ff,j) call gr_theta(t,fparti,ff,theta,j) ENDIF fact_temp = lum / ( 16.d0 * ATAN(1.0d0) * sigma ) fact_temp = fact_temp * ( 1.d0 + 2.30834657d-4 ) ! prov fact_temp = fact_temp / t / ray**2 CX = N_chim(1) * nucleo(1) q = 5.57d-05 * fact_temp / nucleo(j) b = 9.83d-23 * nel * (1./SQRT(t)) / CX DO k=0,NINT(zi(j)) IF(id_ppxi(j)== 1) THEN IF(niv_z(k,j) > 0) THEN bco=7.16d-26*fact_temp*nel/SQRT(t)/(nucleo(j)* 1 niv_z(k,j)*niv_z(k,j)) gc_kj = bco * theta(k) ELSE gc_kj = 0.d0 ENDIF IF(psistar(k,j) > 1.d-30) THEN Ck_s = b * psistar(k,j)**2 gr_kj=q*phistar(k,j)*(1.+xistar(k,j)*C_ion(k,j))* 1 (1.+C_ion(k,j)/Ck_s)**alphstar(k,j)+ acstar(k,j) * gc_kj 2 * (khi(j)/(khi(j)+1.))**betstar(k,j) dgr_kj=popi(k,j)*gr_kj*(xistar(k,j)/(1.d0+ xistar(k,j)*C_ion(k,j)) 1 + alphstar(k,j)/(C_ion(k,j) + Ck_s) ) dgc_kj = betstar(k,j)*gc_kj*(khi(j)**betstar(k,j)) 1 /((khi(j)+1.)**(betstar(k,j)+1)) ELSE gr_kj = acstar(k,j) * gc_kj* (khi(j)/(khi(j)+1.))**betstar(k,j) dgr_kj = 0.d0 dgc_kj = betstar(k,j)*gc_kj*(khi(j)**betstar(k,j)) 1 /((khi(j)+1.d0)**(betstar(k,j)+1)) ENDIF g_rad(j) = g_rad(j) + pond_ion(k,j)*gr_kj dg_rad(j,j) = dg_rad(j,j) + pond_ion(k,j)* (dgr_kj/N_chim(1) + 1 dgc_kj*N_chim(j) ) ELSE !pour les elements non-indentifies IF(igrad == 0) THEN g_rad(j) = 0.d0 ELSEIF(igrad == 1) THEN g_rad(j) = -g_grav ENDIF dg_rad(j,j) = 0.d0 ENDIF ENDDO ! On fait tendre l'acceleration radiative vers zero a partir de T > 1.d7K, ! pour atteindre 0 quand T=1.4d7. x = 1.4d+07/t - 1.d0 IF(x >= 0.4d0) THEN f = 1. ELSEIF(x < 0.d0) THEN f = 0.d0 ELSE f = ((x/0.4d0)**2)*(-2.d0*(x/0.4d0)+3.d0) ENDIF g_rad(j) = f * g_rad(j) dg_rad(j,j) = f * dg_rad(j,j) ENDDO ENDIF RETURN END SUBROUTINE grad2_bp