program lire_cep ! Version: 29/11/99 ! Auteur: T.Guillot ! opa_f : opacite ! convp_jmj : convection ! nuc_0 : reactions thermonucleaires ! lim_0 : conditions limites externes atmosphere ! hopf : loi t(tau) utilisee implicit none integer irep character*30 type(0:5) data type/'ADIABATIQUE','NON-ADIABATIQUE','AD-PROTOPLANETE', & & 'NA-PROTOPLANETE','A-POLYTROPE N=1','R-POLYTROPE N=3'/ external etat_eos,etat_hhe,etat_gph, & & opaque,opa_f,opa_liss, & & convp_jmj,convp_ad,nuc_0, & & lim_0,lim_gg,lim_proto,lim_poly,hopf write(*,*)'Choisir: 0->modele adiabatique' write(*,*)' 1->modele non-adiabatique' write(*,*)' 2->modele de protoplanete adiabatique' write(*,*)' 3->modele de protoplanete non-adiabatique' write(*,*)' 4->polytrope d''indice n=1' write(*,*)' 5->test d''integration du noyau' read(*,*)irep if (irep.eq.0) then call re_cepam(type(0),0,etat_hhe,opaque,convp_jmj,nuc_0, & & lim_0,hopf) ! write(*,*)'Attention: Lim_ext= lim_gg!!!!' ! call re_cepam(type(0),0,etat_hhe,opaque,convp_jmj,nuc_0, ! & lim_gg,hopf) elseif (irep.eq.1) then call re_cepam(type(1),0,etat_hhe,opa_f,convp_jmj,nuc_0, & & lim_0,hopf) ! write(*,*)'Attention: OPACITES LISSEES' ! write(*,*)'Attention: Lim_ext= lim_gg!!!!' ! call re_cepam(type(1),0,etat_hhe,opa_liss,convp_jmj,nuc_0, ! & lim_gg,hopf) elseif (irep.eq.2) then call re_cepam(type(2),0,etat_hhe,opaque,convp_jmj,nuc_0, & & lim_proto,hopf) elseif (irep.eq.3) then ! call re_cepam(type(3),0,etat_hhe,opa_f,convp_jmj,nuc_0, ! & lim_proto,hopf) write(*,*)'Attention: OPACITES LISSEES' write(*,*)'Attention: Lim_ext= lim_gg!!!!' call re_cepam(type(3),0,etat_hhe,opa_liss,convp_jmj,nuc_0, & & lim_gg,hopf) elseif (irep.eq.4) then call re_cepam(type(4),0,etat_gph,opaque,convp_ad,nuc_0, & & lim_poly,hopf) elseif (irep.eq.5) then call test_noyau endif stop end !*********************************************************************** subroutine test_noyau implicit none integer i,j,n_noy,nbm,nbp real*8 p_rac,r_rac,drdp,p_noy(1),t_noy(1),ro_noy(1), & & r_noy(1),m_noy(1),drop_noy(1) include './cepam.inc' include './Communs/ctephys.inc' include './Communs/planetes.inc' include './Communs/modelp.inc' character*31 nom_fich1,nom_fich2,nom_fichn common/in_out/nom_fich1,nom_fich2 common/fich_noy/nom_fichn !======================================================================= call constantes('./Donnees/ctes_jup.don') nom_fich1='bid' nom_fich2='bid' nom_fichn='bid' nbm=10 nbp=20 mplanet=1.d0 p_ice=0.5 omega=1.d0 alpha=0.666666667 precix=1.d-6 10 continue write(*,*)'mnoyau,p_rac?' read(*,*)mnoyau,p_rac call lim_noyau(.false.,p_rac,r_rac,drdp,p_noy,t_noy,ro_noy, & & r_noy,m_noy,drop_noy,n_noy) write(*,'(1p3d13.5)')p_rac,r_rac, & & mnoyau*mjup/(4*pi/3.*(r_rac*rjup)**3) goto 10 do i=1,nbm mnoyau=10**(-4+3*(i-1.)/(nbm-1)) write(*,*)mnoyau do j=1,nbp p_rac=10**(12+2*(j-1.)/(nbm-1)) call lim_noyau(.false.,p_rac,r_rac,drdp,p_noy,t_noy,ro_noy, & & r_noy,m_noy,drop_noy,n_noy) write(*,*)p_rac,r_rac enddo enddo return end !****************************************************************** subroutine list_nuages(type, & & bp_S,q,qt_S,knot_S,n, & & m_tc,mt_tc,mr_tc,n_tc,nr_tc,knot_tc,xchims_t,xchims_tc, & & etat,opa,conv,nuc) ! Version: 05/03/95 ! Interpolation d'un modele CEPAM. ! Creation d'un fichier "nuage.des" pour calculs de condensation ! Auteur: T. Guillot, Departement Cassini, O.C.A., Observatoire de Nice ! type: chaine de caracteres caracterisant le modele implicit none include './cepam.inc' !---- Variables globales ---- integer knot_S,n,n_tc,nr_tc,knot_tc real*8 bp_S(1),q(1),qt_S(1),m_tc(1),mt_tc(1),mr_tc(1), & & xchims_t(1),xchims_tc(1) character*30 type !---- Variables locales ---- integer plnp_S,pds,pn_o parameter (plnp_S=(pn-1)*pm_S+pr_S,pds=(pn_noy-1)*pm_S+pr_S, & & pn_o=1000+pn_noy) integer i,j,k,n_o,imol real*8 p(pn_o),t(pn_o), & & gradad(pn_o),xsat(pn_o),gradwet(pn_o),gradconv(pn_o), & & cp(pn_o),ro(pn_o),grav(pn_o),vconv(pn_o),stab(pn_o) real*8 xchim(1),granr,pmax,fmax real*8 drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2,bid !.... Variables pour interpolation de quantites relatives a l'enveloppe .... integer pn_Sm,ps_Sm parameter (pn_Sm=7,ps_Sm=4) real*8 x_Sm(pn),xt_Sm(pn+2*ps_Sm),xr_Sm(pn),z_Sm(pn_Sm*pn), & & s_Sm(pn_Sm*ps_Sm*pn),dfdx(pn_Sm),fx(pn_Sm) integer n_Sm,nr_Sm,m_Sm,knot_Sm,lx !.... Variables pour interpolation de quantites relatives a l'"atmosphere" .... integer pna,pn_Sa,ps_Sa parameter (pna=100,pn_Sa=3,ps_Sa=8) real*8 x_Sa(pna),xt_Sa(pna+2*ps_Sa),xr_Sa(pna),z_Sa(pn_Sa*pna), & & s_Sa(pn_Sa*ps_Sa*pna) integer n_Sa,nr_Sa,m_Sa,knot_Sa !.... Variables specifiques au calcul de condensation .... integer nbc,i0 parameter (nbc=301) real*8 size(nbc),tcond(nbc),tfall(nbc),tcoag(nbc),tcoal(nbc) real*8 f(pn_o),f_p(pn_o),fc_p(pn_o),t_p(pn_o),tv(pn_o), & & v_p(pn_o),xsat_p(pn_o),lat(pn_o),lat_p(pn_o),tv_p(pn_o), & & cpd_p(pn_o),cpv_p(pn_o),cpc_p(pn_o),cpv(pn_o),phi(pn_o), & & vsed(pn_o),cape(pn_o) real*8 mumoy,mux,eps,humidite,p0,err,v_old,t_old,f_old, & & dc,r_p,ro_p,c1,c2 logical phase character*8 car8 character*16 c_pla,c_mol,met_grad,met_cond integer*4 long include './Communs/ctephys.inc' include './Communs/planetes.inc' include './Communs/modelp.inc' character*2 ci,cj character*4 ck character*80 titre character*31 nom_fich1,nom_fich2 common/in_out/nom_fich1,nom_fich2 external etat,opa,conv,nuc !---------------------------------------------------------------------- granr=kbol/amu !cste des gaz parfaits call idate(i,j,k) ! write(6,*)'idate',i,j,k ! pause'idate' encode(2,'(i2)',ci)i encode(2,'(i2)',cj)j encode(4,'(i4)',ck)k titre=type(:long(type))//' du '//ci//'/'//cj//'/'//ck & & //' T. Guillot, LPL' write(48,*)titre(:long(titre)) write(48,*)'fichier relatif a la condensation' 90 format(i3,13(1x,a3)) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture du fichier nuages.don c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc open(unit=13,file='nuages.don',status='old') read(13,*) read(13,*)car8,c_pla read(13,*)car8,c_mol read(13,*)car8,met_grad read(13,*)car8,met_cond read(13,*)car8,pmax close(13) write(*,*)'pmax=',pmax write(*,*)'molecule: ',c_mol if (c_mol(:3).eq.'CH4') then imol=3 mux=16. elseif (c_mol(:3).eq.'NH3') then imol=4 mux=17. elseif (c_mol(:3).eq.'H2O') then imol=5 mux=18. elseif (c_mol(:2).eq.'Fe') then imol=6 mux=56. elseif (c_mol(:6).eq.'MgSiO3') then imol=7 mux=100. endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Interpolation de l'enveloppe c ! A partir d'un modele calcule par CEPAM, on spline c ! m,r,t,ro,gradconv-gradad,vconv,v_ss en fonction de p c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call psplcep(bp_S,q,qt_S,knot_S,n, & & m_tc,mt_tc,mr_tc,n_tc,nr_tc,knot_tc,xchims_t,xchims_tc, & & x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,knot_Sm,z_Sm,s_Sm, & & etat,opa,conv,nuc) ! Interpolation de la composition chimique en 1-M/Mplanet (supposee cte) !------------------------------------------------------------ ! call pp1dn(nbelem,m_tc,mt_tc,mr_tc,n_tc,nr_tc,mc_S,0.d0, ! & knot_tc,k,dxchim,xchims_t,xchims_tc,xchim,.true.) ! On suppose que la composition chimique est Y=0.18 ! PS: cela n'est pas tout a fait coherent avec le modele obtenu par CEPAM ! PPS: on neglige la presence de NH3 et CH4 xchim(1)=1.-0.18 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Interpolation de l'atmosphere d'apres les observations c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call psplatm(xchim,x_Sa,xt_Sa,xr_Sa,n_Sa,nr_Sa,m_Sa,knot_Sa,z_Sa, & & s_Sa,etat) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Calcul de T(P) supposant Cp=cte c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc n_o=600. !nombre de couches do i=1,n_o p(i)=1d6*(1+pmax*(i-1.)/n_o) if (p(i).le.1d6) then call pp1dn(pn_Sa,x_Sa,xt_Sa,xr_Sa,n_Sa,nr_Sa,m_Sa,log(p(i)), & & knot_Sa,lx,dfdx,z_Sa,s_Sa,fx,.true.) t(i)=exp(fx(1)) call etat(p(i),t(i),xchim,.false., & & ro(i),drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2, & & gradad(i)) gradconv(i)=gradad(i)+fx(3) grav(i)=g*mjup*mplanet/rjup**2 vconv(i)=0.d0 else t(i)=t(i-1)*(p(i)/p(i-1))**gradad(i) err=1d10 do while (err.gt.1d-8) t_old=t(i) call etat(p(i),t(i),xchim,.false., & & ro(i),drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2, & & gradad(i)) t(i)=t(i-1)*(p(i)/p(i-1))**((gradad(i)+gradad(i-1))/2) err=abs(t_old-t(i))/t(i) enddo call pp1dn(pn_Sm,x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,log(p(i)), & & knot_Sm,lx,dfdx,z_Sm,s_Sm,fx,.true.) gradconv(i)=gradad(i)+fx(5) grav(i)=g*fx(1)/fx(2)**2 vconv(i)=fx(6) endif cp(i)=t(i)*dst !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Vitesse de condensation+coagulation+sedimentation c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call rossow78(c_mol,p(i),t(i),ro(i),grav(i),xsat(i),lat(i), & & nbc,size,tcond,tfall,tcoag,tcoal) i0=nbc do j=nbc,1,-1 if (min(tcond(j),tcoag(j),tcoal(j)).gt.tfall(j)) i0=j ! if (tcond(j).gt.tfall(j)) i0=j enddo vsed(i)=p(i)/ro(i)/grav(i)/tfall(i0) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Stabilite c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc eps=mux/(ro(i)*granr*t(i)/p(i)) f(i)=xsat(i)/(1-xsat(i))*eps stab(i)=1-(1-1/eps)/(f(i)+1)*f(i)*mux*lat(i)/granr/t(i) if (10*(i/10)+1.eq.i) & & write(*,'(i4,1p6d10.3)'),i,p(i),t(i),ro(i), & & vsed(i),f(i)/3.6d-5/mux,stab(i) enddo i=2 do while ((stab(i).gt.0).and.(i.lt.n_o)) i=i+1 enddo if (i.lt.n_o) then fmax=(f(i-1)*stab(i)-f(i)*stab(i-1))/(stab(i)-stab(i-1)) write(*,546)fmax,fmax/3.6d-5/mux,p(i-1),p(i) 546 format(' F=',1pd12.5,'; F/Fsol=',1pd12.5,'; P1=',1pd10.3, & & '; P2=',1pd10.3) else write(*,*)'Fmax n''existe pas' endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! POINT DE RETOUR c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 10 continue write(*,*)'Humidite relative (0<=h<=1), dc, pression P0 (bar)?', & & 'Rayon des plumes (km)?' read(*,*)humidite,dc,p0,r_p p0=p0*1d6 r_p=r_p*1d5 mumoy=ro(i)*granr*t(i)/p(i) write(*,*)'Mu=',real(mumoy),'; MuX=',real(mux) write(*,*)'Methodes: ',met_grad,' & ',met_cond write(*,*) eps=mux/mumoy do i=1,n_o f(i)=eps*xsat(i)/(1.-xsat(i))*humidite tv(i)=t(i)*(1+f(i)/eps)/(1+f(i)) enddo !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Recherche de la pression P0, supposee etre la base du nuage c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc i0=1 do while (p(i0).lt.p0) i0=i0+1 enddo write(*,*)'Base du nuage: P0=',real(p(i0)),' ;i0=',i0 call pp1dn(pn_Sm,x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,log(p(i0)), & & knot_Sm,lx,dfdx,z_Sm,s_Sm,fx,.true.) v_p(i0)=fx(7)**2 write(*,*)'Vitesse (Shaviv & Salpeter 1973):',real(sqrt(v_p(i0))) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Plume ascendante: _p c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(*,*)'Dans la plume ascendante: i,p,t_ext,t_p,f_p,fc_p,v2' f_p(i0)=eps*(xsat(i0)/(1.-xsat(i0))) !fraction massique de vapeur fc_p(i0)=0.d0 !fraction massique de condense gradwet(i0)=gradad(i0) t_p(i0)=t(i0) tv_p(i0)=t(i0)*(1+f_p(i0)/eps)/(1+f_p(i0)) tv_p(i0)=tv_p(i0)-tv(i0)*fc_p(i0) phi(i0)=0. i=i0 write(*,'(i4,f8.4,2f8.3,f5.1,1p4d10.3)'),i,p(i)/1d6,tv(i), & & tv_p(i),tv_p(i)-tv(i),f_p(i),fc_p(i),v_p(i), & & sqrt(max(0.d0,v_p(i))) do i=i0-1,1,-1 err=1.d0 v_old=1.d10 t_old=1.d10 f_old=1.d10 t_p(i)=t_p(i+1)*(1+gradwet(i+1)*log(p(i)/p(i+1))) tv_p(i)=t_p(i)*(1+f_p(i+1)/eps)/(1+f_p(i+1)) tv_p(i)=tv_p(i)-tv(i)*fc_p(i) phi(i)=phi(i+1) v_p(i)=v_p(i+1) do while (err.gt.1d-6) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Calculs de condensation c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call clapeyron(c_mol,p(i),t_p(i),xsat_p(i),lat_p(i),phase) f_p(i)=eps*xsat_p(i)/(1-xsat_p(i)) if (met_cond(:7).eq.'yanai73') then ! c1=fc_p(i+1)/(1+f_p(i+1)+fc_p(i+1)) ! c2=fc_p(i)/(1+f_p(i)+fc_p(i)) c1=fc_p(i+1) c2=fc_p(i) c2=c1-(f_p(i)*(1-c2)-f(i)-(f_p(i+1)*(1-c1)-f(i+1)))-phi(i)*( & & f_p(i+1)+(1-f_p(i+1))*c1-2*f(i+1)+ & & f_p(i)+(1-f_p(i))*c2-2*f(i))/2*(p(i)-p(i+1)) c2=max(0.d0,c2) ! fc_p(i)=c2/(1-c2)*(1+f_p(i)) fc_p(i)=c2*dc elseif (met_cond(:8).eq.'vitesses') then if (v_p(i+1).lt.vsed(i+1)) then fc_p(i)=0.d0 else c1=fc_p(i+1) c2=fc_p(i) c2=c1-(f_p(i)*(1-c2)-f(i)-(f_p(i+1)*(1-c1)-f(i+1)))-phi(i)*( & & f_p(i+1)+(1-f_p(i+1))*c1-2*f(i+1)+ & & f_p(i)+(1-f_p(i))*c2-2*f(i))/2*(p(i)-p(i+1)) c2=max(0.d0,c2) ! fc_p(i)=c2/(1-c2)*(1+f_p(i)) fc_p(i)=c2*dc endif else fc_p(i)=fc_p(i+1)-(f_p(i)-f_p(i+1))*dc endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Appelle ETAT c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call etat(p(i),t_p(i),xchim,.false., & & ro_p,drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2,bid) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Calcul des chaleurs specifiques c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cpd_p(i)=t_p(i)*dst if (c_mol(:3).eq.'H2O') then cpv_p(i)=1.87d7 !en erg/g/K cpv(i)=1.87d7 if (phase) then cpc_p(i)=4.19d7 !h2o liquide else cpc_p(i)=2.106d7+0.0073d7*(t_p(i)-273.15) !h2o solide endif elseif (c_mol(:2).eq.'Fe') then cpv_p(i)=0.449d7 cpv(i)=0.449d7 cpc_p(i)=0.449d7 !CRC 5-68 endif ! cpd_p(i)=1.08d8 ! cpv_p(i)=1.08d8 ! cpv(i)=1.08d8 ! cpc_p(i)=1.08d8 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Gradient humide d'apres Emanuel (1994) et Awal & Lunine (1986) c ! On suppose que le gradient exterieur tient compte de f c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (r_p.le.0.) then c1=0. else phi(i)=-0.2/r_p/grav(i)/ro(i) c1=phi(i)*((t_p(i)-t(i))/(t(i)/p(i)*gradad(i))* & & (cp(i)+f(i)*cpv(i))/(cpd_p(i)+f_p(i)*cpv_p(i))+ & & lat_p(i)/(t(i)/p(i)*gradad(i))/ & & (cpd_p(i)+f_p(i)*cpv_p(i))*(f_p(i)-f(i))) endif if (met_grad(:9).eq.'emanuel94') then gradwet(i)=gradad(i)/((1+f(i))/(1+f(i)*cpv(i)/cp(i)))* & & ((1+f_p(i))/(1+f_p(i)*cpv_p(i)/cpd_p(i)))* & & (1+lat_p(i)*f_p(i)*mumoy/granr/t_p(i)-c1)/ & & (1+fc_p(i)*cpc_p(i)/(cpd_p(i)+f_p(i)*cpv_p(i))+ & & ((lat_p(i)/t_p(i))**2*mux/granr *f_p(i)*(1+f_p(i)/eps)/ & & (cpd_p(i)+f_p(i)*cpv_p(i)))) else if (met_grad(:6).eq.'awal94') then !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Gradient humide d'apres Awal & Lunine (1986) c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! gradwet(i)=gradad(i)*(1+lat_p(i)*f_p(i)*mumoy/granr/tv_p(i)- ! & (t(i)-t_p(i))/(t(i)/p(i)*gradad(i))*phi(i)- ! & lat_p(i)/(t(i)/p(i)*gradad(i))/cp(i)*(f_p(i)-f(i))* ! & phi(i))/(1+(lat_p(i)/t_p(i))**2*mux/granr/cp(i)*f_p(i)) gradwet(i)=gradad(i)*(1+lat_p(i)*f_p(i)*mumoy/granr/tv_p(i)- & & c1)/(1+(lat_p(i)/t_p(i))**2*mux/granr/cp(i)*f_p(i)) else write(*,*)'Quelle est cette methode?' pause goto 10 endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Temperature de la plume c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc t_p(i)=t_p(i+1)+(t_p(i+1)*gradwet(i+1)+t_p(i)*gradwet(i))* & & log(p(i)/p(i+1))/2 if (t_p(i).lt.t(i)/2) t_p(i)=t(i)/2 tv_p(i)=t_p(i)*(1+f_p(i)/eps)/(1+f_p(i)) ! tv_p(i)=tv_p(i)-tv(i)*fc_p(i)/(f_p(i)+1) tv_p(i)=tv_p(i)-tv(i)*fc_p(i) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Vitesse ascentionnelle c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! c1=(f_p(i+1)+1)/(f_p(i+1)+fc_p(i+1)+1) ! c2=(f_p(i)+1)/(f_p(i)+fc_p(i)+1) c1=1. c2=1. v_p(i)=max(0.d0,v_p(i+1))+(c1*(tv(i+1)-tv_p(i+1))/tv(i+1)/ & & ro(i+1)+c2*(tv(i)-tv_p(i))/tv(i)/ro(i))*(p(i)-p(i+1)) if (r_p.gt.0.) then v_p(i)=v_p(i)-(phi(i)*v_p(i)+phi(i+1)*max(0.d0,v_p(i+1)))* & & (p(i)-p(i+1)) ! Test: comparaison avec l'expression obtenue par Stoker ! write(*,*)'v_p1:',v_p(i) ! v_p(i)=((granr/p(i)/mumoy*(tv_p(i)-tv(i))+phi(i)* ! & max(0.d0,v_p(i+1)))* ! & exp(-2*phi(i)*(p(i)-p(i+1)))-(granr/p(i)/mumoy)* ! & (tv_p(i)-tv(i)))/phi(i) ! write(*,*)'v_p2:',v_p(i) endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! CAPE: conditional available potential energy c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cape(i)=cape(i+1)+log(p(i+1)/p(i))*granr/mumoy*(tv_p(i)-tv(i)+ & & tv_p(i+1)-tv(i+1))/2 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Erreur commise lors de la presente iteration c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (abs(v_p(i)).lt.1d-2) then err=abs(v_old-v_p(i)) else err=abs((v_old-v_p(i))/v_p(i)) endif err=err+abs(t_old-t_p(i))/t_p(i) if (abs(f_p(i)).lt.1d-8) then err=err+abs(f_old-f_p(i)) else err=err+abs((f_old-f_p(i))/f_p(i)) endif v_old=v_p(i) t_old=t_p(i) f_old=f_p(i) enddo if (.true.) !(i.eq.10*(i/10)+1) & & write(*,'(i4,f8.4,2f8.3,f5.1,1p4d10.3)'),i,p(i)/1d6,tv(i), & & tv_p(i),tv_p(i)-tv(i),f_p(i),fc_p(i),v_p(i), & & sqrt(max(0.d0,v_p(i))) enddo open(unit=11,file='nuages.asc',status='unknown') do i=1,i0 if (i.eq.10*(i/10)+1) write(11,110) 110 format(/,' N',t5,'P',t16,'grav',t27,'Rho',/, & & t5,'T_ext',t16,'Tv_ext',t27,'Gradad',t38,'f_ext',t49, & & 'L_ext',t60,'Cpd_ext',t71,'Cpv_ext',/, & & t5,'T_plu',t16,'Tv_plu',t27,'Gradwet',t38,'f_plu',t49, & & 'L_plu',t60,'Cpd_plu',t71,'Cpv_plu',/, & & t5,'Delta_Tv',t16,'Fb(Tv)',t27,'fc_plu',t38,'Cpc_plu', & & t49,'V_plu',t60,'V_CAPE',t71,'Phi') write(11,111)i,p(i),grav(i),ro(i), & & t(i),tv(i),gradad(i),f(i),lat(i),cp(i),cpv(i), & & t_p(i),tv_p(i),gradwet(i),f_p(i),lat_p(i),cpd_p(i), & & cpv_p(i), & & tv_p(i)-tv(i),grav(i)*(tv_p(i)-tv(i))/tv(i), & & fc_p(i),cpc_p(i),sqrt(max(0.d0,v_p(i))), & & sqrt(max(0.d0,2*cape(i))),phi(i) 111 format(i3,1p3d11.4,/,3x,1p6d11.4,1pd10.3,/,3x,1p6d11.4,1pd10.3, & & /,3x,1p6d11.4,1pd10.3) enddo close(11) write(*,*)'J''ai cree le fichier "nuages.asc"' open(unit=12,file='nuages.des',status='unknown') write(12,*)'**** JOVIAN H2O CLOUDS ****' write(12,*)'* P(bar),Tv(K),Tv_cloud(K),Delta_Tv(K),f,', & & 'f_cloud,fc_cloud,V_cloud(cm/s)' write(12,*)i0 do i=i0,1,-1 write(12,'(f8.4,2f8.3,1x,f8.3,1p4d10.3)'),p(i)/1d6,tv(i), & & tv_p(i),tv_p(i)-tv(i),f(i),f_p(i),fc_p(i), & & sqrt(max(0.d0,v_p(i))) enddo close(12) write(*,*)'J''ai cree le fichier "nuages.des"' goto 10 return end !********************************************************************* subroutine psplcep(bp_S,q,qt_S,knot_S,n, & & m_tc,mt_tc,mr_tc,n_tc,nr_tc,knot_tc,xchims_t,xchims_tc, & & x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,knot_Sm,z_Sm,s_Sm, & & etat,opa,conv,nuc) ! Version: 12/03/95 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Interpolation de l'enveloppe c ! A partir d'un modele calcule par CEPAM, on spline c ! m,r,t,ro,gradconv-gradad,vconv,v_ss en fonction de p c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Entrees: ! bp_S,q,qt_S,knot_S: coefficients de spline pour interpolation de l'enveloppe ! n: nombre de couches ! m_tc,mt_tc,mr_tc,n_tc,nr_tc,knot_tc,xchims_t,xchims_tc: splines pour composition ! chimique ! etat,opa,conv,nuc: subroutines externes (EOS...etc) ! Sorties: ! x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,knot_Sm,z_Sm,s_Sm: splines recherchees implicit none include './cepam.inc' !---- Variables globales ---- integer knot_S,n,n_tc,nr_tc,knot_tc real*8 bp_S(1),q(1),qt_S(1),m_tc(1),mt_tc(1),mr_tc(1), & & xchims_t(1),xchims_tc(1) integer pn_Sm,ps_Sm parameter (pn_Sm=7,ps_Sm=4) real*8 x_Sm(pn),xt_Sm(pn+2*ps_Sm),xr_Sm(pn),z_Sm(pn_Sm*pn), & & s_Sm(pn_Sm*ps_Sm*pn),dfdx(pn_Sm),fx(pn_Sm) integer n_Sm,nr_Sm,m_Sm,knot_Sm,lx external etat,opa,conv,nuc !---- Variables locales ---- integer plnp_S,pds,pn_o parameter (plnp_S=(pn-1)*pm_S+pr_S,pds=(pn_noy-1)*pm_S+pr_S, & & pn_o=1000+pn_noy) integer i,j,k real*8 lnp_S(plnp_S),lnt_S(plnp_S),r2_S(plnp_S),l23_S(plnp_S), & & mu_S(plnp_S), & & p(pn),t(pn),r(pn),l(pn),m(pn), & & gradad(pn),gradconv(pn), & & ro(pn),vconv(pn),v_ss(pn),grav(pn),del_t(pn) real*8 rplanet,qn,interq,interm real*8 xchim(1),dxchim(1),cp, & & gradient,dgradp,dgradt,dgradl,dgradr,dgradm,dgradx, & & epsilo,depsp,depst,depsx,kap,dkapp,dkapt,dkapx, & & delta,deltap,deltat,deltax,dcpp,dcpt,dcpx, & & hp,gradrad,d_grad,q_o real*8 drop,drot,drox,u,s,dsp,dst,dsx,xh,xh2 !---- Communs ---- include './Communs/ctephys.inc' include './Communs/planetes.inc' include './Communs/modelp.inc' !===================================================================== do i=1,(n-1)*m_S+r_S !modele au temps t lnp_S(i)=bp_S(ne*(i-1)+1) lnt_S(i)=bp_S(ne*(i-1)+2) r2_S(i) =bp_S(ne*(i-1)+3) l23_S(i)=bp_S(ne*(i-1)+4) mu_S(i)=bp_S(ne*(i-1)+5) enddo rplanet=sqrt(interm(r2_S,mu_S,q,qt_S,m_S,r_S,knot_S,1.d0,n, & & .false.,qn)) write(*,*)'rplanet=',rplanet do i=1,n q_o=i p(i)=interq(lnp_S,q,qt_S,m_S,r_S,knot_S,q_o,n) t(i)=interq(lnt_S,q,qt_S,m_S,r_S,knot_S,q_o,n) r(i)=interq(r2_S ,q,qt_S,m_S,r_S,knot_S,q_o,n) l(i)=interq(l23_S,q,qt_S,m_S,r_S,knot_S,q_o,n) m(i)=interq(mu_S ,q,qt_S,m_S,r_S,knot_S,q_o,n) p(i)=exp(p(i)) t(i)=exp(t(i)) r(i)=sqrt(r(i)) l(i)=sqrt(l(i)**3) m(i)=sqrt(1.d0-m(i))**3 ! write(*,2001)i,q_o,p(i),t(i),r(i),l(i),m(i) ! 2001 format(i3,6(1x,1pd10.3)) ! Interpolation de la composition chimique en 1-M/Mplanet !------------------------------------------------------------ call pp1dn(nbelem,m_tc,mt_tc,mr_tc,n_tc,nr_tc,mc_S,1.d0-m(i), & & knot_tc,k,dxchim,xchims_t,xchims_tc,xchim,.true.) ! Calcul de rho et autres quantites thermodynamiques !------------------------------------------------------- call thermos(p(i),t(i),xchim,m(i),l(i),r(i),.true.,dxchim, & & ro(i),drop,drot,drox,u,s,dsp,dst,dsx, & & gradient,dgradp,dgradt,dgradl,dgradr,dgradm,dgradx, & & epsilo,depsp,depst,depsx,kap,dkapp,dkapt,dkapx, & & delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, & & hp,gradad(i),gradrad,gradconv(i),d_grad, & & xh,xh2,etat,opa,conv,nuc) grav(i)=g*m(i)*mjup*mplanet/(r(i)*rjup)**2 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Vitesse de convection d'apres la MLT c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! write(*,*)i,p(i),drot/ro(i)*t(i) vconv(i)=sqrt(grav(i)*(-drot)/ro(i)*t(i)*alpha**2*hp/8* & & (gradconv(i)-gradad(i))) enddo do i=1,n v_ss(i)=0.d0 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Vitesses convectives: formulation de Shaviv & Salpeter (1973) c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! trouve a partir de i=r2 le point j=r1 situe a une Hp j=1 do while ((p(j).lt.p(i)*exp(1.0)).and.(j.lt.n)) j=j+1 enddo del_t(j)=0. ! Ecart de temperature !------------------------- do k=j-1,i,-1 del_t(k)=del_t(k+1)+(gradconv(k+1)-gradad(k+1)+ & & gradconv(k)-gradad(k))/2*(log(p(k+1)-p(k))) enddo ! Vitesse convective !----------------------- v_ss(i)=0.d0 do k=j-1,i,-1 v_ss(i)=v_ss(i)+(grav(k)*del_t(k)+grav(k+1)*del_t(k+1))* & & (r(k)-r(k+1))*rjup enddo v_ss(i)=sqrt(v_ss(i)) ! write(*,'(2i4,1p4d10.3)'),i,j,p(i),t(i),vconv(i),v_ss(i) enddo !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! On spline m,r,t,ro,gradconv-gradad,vconv,v_ss en fonction de p c ! Attention: tout est en log sauf m,r,grad,vconv et v_ss c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,n x_Sm(i)=log(p(i)) z_Sm(pn_Sm*(i-1)+1)=m(i)*mjup*mplanet z_Sm(pn_Sm*(i-1)+2)=r(i)*rjup z_Sm(pn_Sm*(i-1)+3)=log(t(i)) z_Sm(pn_Sm*(i-1)+4)=log(ro(i)) z_Sm(pn_Sm*(i-1)+5)=gradconv(i)-gradad(i) z_Sm(pn_Sm*(i-1)+6)=vconv(i) z_Sm(pn_Sm*(i-1)+7)=v_ss(i) enddo n_Sm=n m_Sm=4 !splines d'ordre 4 (cubiques) call pp1dn(pn_Sm,x_Sm,xt_Sm,xr_Sm,n_Sm,nr_Sm,m_Sm,x_Sm(1), & & knot_Sm,lx,dfdx, & & z_Sm,s_Sm,fx,.false.) return end !********************************************************************* subroutine psplatm(xchim,x_Sa,xt_Sa,xr_Sa,n_Sa,nr_Sa,m_Sa,knot_Sa,& & z_Sa,s_Sa,etat) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Interpolation de l'atmosphere c ! D'apres les observations de Voyager (cf. Lindal 1992) c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Entrees: ! xchim(1): fraction massique d'hydrogene ! etat: routine d'equation d'etat ! Sorties: ! x_Sa,xt_Sa,xr_Sa,n_Sa,nr_Sa,m_Sa,knot_Sa,z_Sa,s_Sa: splines pour l'atmosphere implicit none include './cepam.inc' !---- Variables globales ---- real*8 xchim(1) integer pna,pn_Sa,ps_Sa parameter (pna=100,pn_Sa=3,ps_Sa=8) real*8 x_Sa(pna),xt_Sa(pna+2*ps_Sa),xr_Sa(pna),z_Sa(pn_Sa*pna), & & s_Sa(pn_Sa*ps_Sa*pna),dfda(pn_Sa),fa(pn_Sa) integer n_Sa,nr_Sa,m_Sa,knot_Sa,lx external etat !---- Variables locales ---- integer i,j,nf real*8 p(pna),t(pna),grad(pna),dgrad(pna),gradient(pna),tl(pna), & & gradad(pna),ro(pna) real*8 pmbar,tk,drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2 character*40 fichier data fichier/'./Donnees/atmos.dat'/ !===================================================================== !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture du fichier "fichier" c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc open(unit=11,file=fichier,status='old') do i=1,3 read(11,*) enddo read(11,*)nf do i=1,nf read(11,*)pmbar,tk if (tk.gt.0) then j=j+1 p(j)=pmbar*1d3 t(j)=tk endif enddo n_Sa=j close(11) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! On lisse le profil de temperature c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc grad(1)=log(t(1)/t(2))/log(p(1)/p(2)) do i=2,n_Sa-1 grad(i)=log(t(i-1)/t(i+1))/log(p(i-1)/p(i+1)) enddo grad(n_Sa)=log(t(n_Sa-1)/t(n_Sa))/log(p(n_Sa-1)/p(n_Sa)) do i=1,n_Sa call etat(p(i),t(i),xchim,.false., & & ro(i),drop,drot,drox,drott,drotp,drotx, & & u,s,dsp,dst,dsx,dstt,dstp,dspp,dstx,dspx,xh,xh2, & & gradad(i)) dgrad(i)=grad(i)-gradad(i) enddo !!!!! Arbitrairement, on suppose que pour les 3 derniers pts, grad=gradad !!!! !!!!! (ceci est justifie par le faible ecart entre ces deux valeurs) do i=n_Sa,n_Sa-4,-1 dgrad(i)=0.d0 enddo ! m_S=5 est la valeur qui semble lisser le + convenablement call contour(1,p,n_Sa,5,dgrad,xt_Sa,knot_Sa,.false.,lx,p(1), & & fa,dfda) do i=1,n_Sa call contour(1,p,n_Sa,5,dgrad,xt_Sa,knot_Sa,.true.,lx,p(i), & & fa,dfda) gradient(i)=gradad(i)+fa(1) enddo tl(n_Sa)=t(n_Sa) do i=n_Sa-1,1,-1 tl(i)=tl(i+1)*(p(i)/p(i+1))**((gradient(i)+gradient(i+1))/2) enddo ! open(unit=12,file='psplatm.dat',status='unknown') ! write(12,*)n_Sa ! do i=n_Sa,1,-1 ! write(12,'(1p8d11.3)')log10(p(i))-6,t(i),tl(i), ! & grad(i)-gradad(i),gradient(i)-gradad(i) ! enddo ! close(12) ! write(*,*)'J''ai cree le fichier "psplatm.dat"' !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! On spline t,ro,gradient-gradad en fonction de p c ! Attention: tout est en log sauf m,r,grad,vconv et v_ss c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,n_Sa x_Sa(i)=log(p(i)) z_Sa(pn_Sa*(i-1)+1)=log(tl(i)) z_Sa(pn_Sa*(i-1)+2)=log(ro(i)) z_Sa(pn_Sa*(i-1)+3)=gradient(i)-gradad(i) enddo m_Sa=4 !splines d'ordre 4 (cubiques) call pp1dn(pn_Sa,x_Sa,xt_Sa,xr_Sa,n_Sa,nr_Sa,m_Sa,x_Sa(1), & & knot_Sa,lx,dfda, & & z_Sa,s_Sa,fa,.false.) return end