program cepamdon ! Version: 29/11/99 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture/reecriture du fichier cepam.don c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Auteur: T. Guillot, OCA implicit none include './cepam.inc' integer i,j,irep,ii,nfic,irad integer ind(5) real*8 x0,dx0,mnoyau0,p_ice0,req_ct,j2_ct,j4_ct,j6_ct,dj4,yhe real*8 dx(5),x(5),mc(5) character*40 dir,fic,ficdes,ficdon(5) integer*4 leng real*8 mplanet,mnoyau,p_ice,omega,precix, & & agemax,age_list(pnlist),mdt,dtmax, & & nucleo(pnelem),alpha,vsal, & & ctep,ctet,cter,ctel,ctem, & & pc0,tc0,ltot0,rtot0 integer silent,ncouche,m_S,mc_S,iter_max logical rapid,en_masse,ppt character*31 fich(8) character*40 fichier_ctes character*31 opa1,opa2,opa3,opa4 include './Communs/planetes.inc' namelist/nl_sorties/silent namelist/nl_espace/mplanet,mnoyau,p_ice,omega,ncouche,rapid,m_S, & & mc_S,en_masse,ppt,precix,iter_max namelist/nl_ctes/fichier_ctes namelist/nl_temps/agemax,age_list,mdt,dtmax namelist/nl_chim/nucleo namelist/nl_conv/alpha,vsal namelist/nl_etat/fich namelist/nl_opa/opa1,opa2,opa3,opa4 namelist/nl_repartit/ctep,ctet,cter,ctel,ctem namelist/nl_initial/pc0,tc0,ltot0,rtot0 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Formules pour approximation dJ4=0 c ! Attention: ici, DX=Xmol-Xmet donc generallement DX>0 c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Jupiter real*8 dxjup,xjup,mcjup,dxjupmax dxjup(dj4)=-dj4/30.9 xjup(dx0)=0.6016*dx0 mcjup(dx0)=-0.244*dx0 dxjupmax(mnoyau0)=mnoyau0/0.244 ! dxjup(dj4)=-dj4/35.4 ! xjup(dx0)=0.5018*dx0 ! mcjup(dx0)=-0.2605*dx0 ! dxjupmax(mnoyau0)=mnoyau0/0.2605 ! Saturne real*8 dxsat,xsat,mcsat,dxsatmax dxsat(dj4)=-dj4/2.68 xsat(dx0)=0.1694*dx0 mcsat(dx0)=-0.339*dx0 dxsatmax(mnoyau0)=mnoyau0/0.339 ! dxsat(dj4)=-0.30332+sqrt(0.09200-0.33910*dj4) ! xsat(dx0)=dx0*(0.04373+dx0*0.10486) ! mcsat(dx0)=dx0*(-0.2035-dx0*0.3757) ! dxsatmax(mnoyau0)=(-0.2035+sqrt(0.2035**2+4*0.3757*mnoyau0))/2/ ! & 0.3757 !---------------------------------------------------------------------- write(*,100) 100 format(/,'LECTURE/ECRITURE D''UN FICHIER DE TYPE CEPAM.DON',/, & & '1 -> Modification de masse pour calculs d''evolution',/, & & '2 -> Modification de DX pour optimisation de Jupiter',/, & & '3 -> Modification de DX pour optimisation de Saturne',/) read(*,*)irep if (irep.eq.1) then fic='cepam.don' elseif (abs(irep).eq.2) then write(*,101)'xjq0.don' 101 format('Indiquez le repertoire ou rechercher le fichier ',a) read(*,'(a)')dir fic=dir(:leng(dir))//'/xjq0.don' ficdes=dir(:leng(dir))//'/xjq0.des' elseif (abs(irep).eq.3) then write(*,101)'xsq0.don' read(*,'(a)')dir fic=dir(:leng(dir))//'/xsq0.don' ficdes=dir(:leng(dir))//'/xsq0.des' else write(*,*)'Entrer un chiffre de 1 a 3' stop endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture du fichier "cepam.don" c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc open(unit=3,form='formatted',status='old',file=fic) read(3,nl_sorties) read(3,nl_espace) read(3,nl_ctes) read(3,nl_temps) read(3,nl_chim) read(3,nl_conv) read(3,nl_etat) read(3,nl_opa) read(3,nl_repartit) read(3,nl_initial) close(unit=3) write(*,12)mplanet*mjup/1.89861120e30, & & mplanet*mjup/5.976e27,mnoyau*mjup/5.976e27,p_ice,omega, & & nucleo(1),nucleo(2),alpha,PPT 12 format(30('*'),' CEPAM.DON initial ',30('*'),/, & & ' mplanet=',1p,d12.5,' Mjup',/, & & ' mplanet=',0p,f10.4,' Mt; mnoyau=',f9.4,' Mt; p_ice=',f7.4, & & '; omega=',f9.4,/, & & ' xmol=',f11.8,'; dx=',f11.8,'; alpha=',f10.6,'; PPT=',l1,/, & & 79('*')) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture du fichier des constantes c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call constantes(fichier_ctes) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! IREP=2 ou 3: c ! lecture du fichier .des et calcul du dx c ! elimination des fichiers "xjq.don" redondants c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (abs(irep).ge.2) then write(*,*)'Modele adiabatique (->0) ou radiatif (->1)?' read(*,*)irad if (abs(irep).eq.3) then write(*,*)'Modele de Saturne: entrer Yhe (ex: 0.06 ou 0.11)' read(*,*)yhe endif open(unit=26,form='formatted',status='old',file=ficdes) read(26,148)x0,dx0,mnoyau0,p_ice0,req_ct,j2_ct,j4_ct,j6_ct 148 format(////, & & 'C X=',f11.8,' DX=',f11.8,' Mcoeur=',f11.8, & & ' frac_glaces=',f11.8,/, & & /////,'C Cal/Th',t10,4(1x,1pd13.6)) close(unit=26) dj4=(j4_ct-1.d0)*j4jup/sig_j4 if (dx0.ne.0.) then write(*,*)'Probleme: dx<>0 :',dx ! pause stop 1 endif if (abs(irep).eq.2) then dx0=dxjupmax(mnoyau0) dx(1)=min(dx0,dxjup(dj4)) dx(2)=min(dx0,dxjup(dj4+1)) dx(3)=min(dx0,dxjup(dj4-1)) dx(4)=min(dx0,dxjup(dj4+1.8)) !rot. solide pour Jupiter dx(5)=min(dx0,dxjup(dj4-0.2)) !rot. solide pour Jupiter do i=1,5 x(i)=x0+xjup(dx(i)) mc(i)=max(0.d0,mnoyau0+mcjup(dx(i))-1d-6) enddo elseif (abs(irep).eq.3) then dx0=dxsatmax(mnoyau0) dx(1)=min(dx0,dxsat(dj4)) dx(2)=min(dx0,dxsat(dj4+1)) dx(3)=min(dx0,dxsat(dj4-1)) dx(4)=min(dx0,dxsat(dj4+1.5)) !rot. solide pour Saturne dx(5)=min(dx0,dxsat(dj4-0.5)) !rot. solide pour Saturne do i=1,5 x(i)=x0+xsat(dx(i)) mc(i)=max(0.d0,mnoyau0+mcsat(dx(i))-1d-6) enddo endif write(*,110)dj4,x0,mnoyau0,'Optimal',dx(1),x(1),mc(1), & & 'Rot. diff -',dx(2),x(2),mc(2), & & 'Rot. diff +',dx(3),x(3),mc(3), & & 'Rot. solide -',dx(4),x(4),mc(4), & & 'Rot. solide +',dx(5),x(5),mc(5) 110 format('Parametres lus:',t25,'dJ4',t40,'X',t55,'Mnoyau',/, & & t23,f8.4,t38,f12.8,t53,f12.8,/, & & 'Parametres deduits:',t25,'dX',t40,'X',t55,'Mnoyau',/, & & 5(1x,a,t23,f12.8,t38,f12.8,t53,f12.8,/)) nfic=0 do i=1,5 dx(i)=-dx(i) ! Jupiter ------------------------------ if (abs(irep).eq.2) then if (mc(i).gt.0) then if (-dx(i).ge.0.) then j=int(-dx(i)*100) ii=nint(-dx(i)*1000)-10*int(-dx(i)*100) if (ii.ge.10) then j=j+1 ii=ii-10 endif write(ficdon(nfic+1),120)dir(:leng(dir)),j,ii 120 format(a,'/xjq',i1,',',i1) else j=int(dx(i)*100) ii=nint(dx(i)*1000)-10*int(dx(i)*100) if (ii.ge.10) then j=j+1 ii=ii-10 endif write(ficdon(nfic+1),121)dir(:leng(dir)),j,ii 121 format(a,'/xjq-',i1,',',i1) endif else ficdon(nfic+1)=dir(:leng(dir))//'/xjq-nc' endif ! Saturne ------------------------------ elseif (abs(irep).eq.3) then if (mc(i).gt.0) then ii=int(log10(abs(-dx(i)*100)))+1 if (-dx(i).lt.0.) ii=ii+1 ! write(ficdon(nfic+1),122)dir(:leng(dir)),nint(-dx(i)*100) ! 122 format(a,'/xsq',i) if (ii.eq.1) write(ficdon(nfic+1),1221) & & dir(:leng(dir)),nint(-dx(i)*100) if (ii.eq.1) write(ficdon(nfic+1),1222) & & dir(:leng(dir)),nint(-dx(i)*100) if (ii.eq.1) write(ficdon(nfic+1),1223) & & dir(:leng(dir)),nint(-dx(i)*100) if (ii.eq.1) write(ficdon(nfic+1),1224) & & dir(:leng(dir)),nint(-dx(i)*100) 1221 format(a,'/xsq',i1) 1222 format(a,'/xsq',i2) 1223 format(a,'/xsq',i3) 1224 format(a,'/xsq',i4) else ficdon(nfic+1)=dir(:leng(dir))//'/xsq-nc' endif endif ! Elimination des fichiers redondants ---- ind(nfic+1)=i j=1 do while ((j.ne.nfic).and.(ficdon(j).ne.ficdon(nfic+1))) j=j+1 enddo if ((nfic.eq.0).or.(j.eq.nfic).and. & & (ficdon(j).ne.ficdon(nfic+1))) nfic=nfic+1 enddo ! do i=1,nfic ! write(*,*)i,ficdon(i) ! enddo endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Reecriture d'un nouveau "cepam.don" c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (irep.eq.1) then write(*,*)'Entrer le nouveau mplanet (ex: 1.0)' read(*,*)mplanet write(*,13)mplanet*mjup/1.89861120e30 13 format('-----> mplanet=',1p,d12.5,' Mjup') ! open(unit=3,form='formatted',status='old',file='cepam.don', ! & delim='quote') open(unit=3,form='formatted',status='old',file='cepam.don') write(3,nl_sorties) write(3,nl_espace) write(3,nl_ctes) write(3,nl_temps) write(3,nl_chim) write(3,nl_conv) write(3,nl_etat) write(3,nl_opa) write(3,nl_repartit) write(3,nl_initial) close(unit=3) else open(unit=4,form='formatted',status='unknown',file='ex_optc.exe') write(4,444)dir(:leng(dir)) 444 format('#* Optimisations',/, & & '#* Executable cree par le code cepamdon.f',/, & & 'unalias cp',/, & & 'cp ',a,'/etat.don .') write(*,*)'Ecritures:' do i=1,nfic nucleo(1)=x(ind(i)) nucleo(2)=dx(ind(i)) mnoyau=mc(ind(i)) write(*,143)ficdon(i)//'.don',nucleo(1),nucleo(2),mnoyau 143 format(a30,'-> X=',f11.8,'; DX=',f11.8,'; Mc=',f11.8) if (irep.ge.2) then ! open(unit=3,form='formatted',status='unknown', ! & file=ficdon(i)(:leng(ficdon(i)))//'.don',delim='quote') open(unit=3,form='formatted',status='unknown', & & file=ficdon(i)(:leng(ficdon(i)))//'.don') write(3,nl_sorties) write(3,nl_espace) write(3,nl_ctes) write(3,nl_temps) write(3,nl_chim) write(3,nl_conv) write(3,nl_etat) write(3,nl_opa) write(3,nl_repartit) write(3,nl_initial) close(unit=3) if (mnoyau.gt.0) then write(4,144)irad,ficdon(i)(:leng(ficdon(i))),-nucleo(2)*100 144 format('opt ',i1,1x,a,1x,f12.8) else write(4,146)irad,ficdon(i)(:leng(ficdon(i))),-nucleo(2)*100 146 format('optnc ',i1,1x,a,1x,f12.8) endif endif write(4,243)ficdon(i)(:leng(ficdon(i)))//'.don' 243 format('cp ',a,' cepam.don') write(4,244)irad,ficdon(i)(:leng(ficdon(i))),yhe 244 format('zel ',i1,1x,a,1x,f10.8) enddo ! write(4,147)(ficdon(i)(:leng(ficdon(i)))//'.optim',i=1,nfic) !147 format(/,'echo " "',/,'echo ',70('0'),/, ! & 'mtail',(1x,a),' ') 147 format(/,'echo " "',/,'echo ',70('0'),/, & & 'mtail',10(1x,a),' ') write(*,145) 145 format(/,'Pour lancer l''execution, taper:',/, & & 'chmod u+x ex_optc.exe',/,'ex_optc.exe') endif end