!*************************************************************************** program jacob_cep ! Version: 00/00/03 ! calcul pour calibration du jacobien suivant ! ( dR/dX dR/dMc dR/dGl ) ! JAC= ( dJ2/dX dJ2/dMc dJ2/dGl ) ! ( dJ4/dX dJ4/dMc dJ4/dGl ) ! ou R= Rayon equatorial, J2,J4= moments gravitationnels, ! X= abondance de H, Mc=masse du noyau, Gl= fraction de glaces ! Auteur: T.Guillot, Departement J.D. Cassini, O.C.A., Observatoire de Nice implicit none include './cepam.inc' integer ncouche,m_S,mc_S,iter_max,i,j,irep,silent integer*4 leng real*8 mplanet,mnoyau,p_ice,omega,alpha,mdt,dtmax,dtlist,x0,dx0,z,& & agemax,age_list(pnlist),dt0,precix,precit,d_grav,vsal, & & mnoyau0,p_ice0,nucleo(pnelem),ctep,ctet,cter,ctel,ctem, & & pc0,tc0,ltot0,rtot0,req_ct,j2_ct,j4_ct,j6_ct, & & jac(3*3),jac_er(3*3),a(3,7),b(3,7),j12,j23,h logical en_masse,ppt,rapid character*31 fich(8) character*40 fichier_ctes character*31 opa1,opa2,opa3,opa4 character*80 nom,nom1,no 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 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(*,*)'entrer 1 2 3 4 5 6 7 (8)' write(*,*)' X+h X-h Mc+h Mc-h Gl+h Gl-h jacobien' read(*,*)irep write(*,*)irep if (irep.eq.8) goto 1000 write(*,*)'nom du fichier .des ?' read(*,'(a)')nom write(*,*)nom(:leng(nom)) if (irep.eq.1) then nom1=nom(:leng(nom))//'.des' else nom1=nom(:leng(nom))//'_tempo.des' endif open(unit=25,form='formatted',status='old',file=nom1) read(25,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)) write(*,148)x0,dx0,mnoyau0,p_ice0,req_ct,j2_ct,j4_ct,j6_ct close(unit=25) open(unit=26,form='formatted',status='unknown', & & file='jacob_cep.tempo') 10 continue read(26,*,end=15) goto 10 15 continue if (dx0.eq.0.) then write(26,226)x0,mnoyau0,p_ice0,req_ct,j2_ct,j4_ct,j6_ct else write(26,226)x0,mnoyau0,dx0,req_ct,j2_ct,j4_ct,j6_ct endif 226 format(7f11.8) close(26) nom1=nom(:leng(nom))//'.des' open(unit=25,form='formatted',status='old',file=nom1) read(25,149)x0,dx0,mnoyau0,p_ice0 149 format(////, & & 'C X=',f11.8,' DX=',f11.8,' Mcoeur=',f11.8, & & ' frac_glaces=',f11.8) close(unit=25) 1000 continue if (irep.lt.7) then no=nom(:leng(nom))//'.don' open(unit=25,form='formatted',status='old',file=no) read(25,nl_sorties) read(25,nl_espace) read(25,nl_ctes) read(25,nl_temps) read(25,nl_chim) read(25,nl_conv) read(25,nl_etat) read(25,nl_opa) read(25,nl_repartit) read(25,nl_initial) close(unit=25) ! write(*,*)'Nombre de couches?' ! read(*,*)ncouche ! write(*,*)'Rapid (T/F)?' ! read(*,*)rapid if (irep.lt.3) then h=1.d0+1.d-3 else h=1.d0+1.d-2 endif if (irep.eq.1) then nucleo(1)=x0*h elseif (irep.eq.2) then nucleo(1)=x0/h elseif (irep.eq.3) then mnoyau=mnoyau0*h elseif (irep.eq.4) then mnoyau=mnoyau0/h elseif (irep.eq.5) then if (dx0.eq.0.) then p_ice=p_ice0*h else nucleo(2)=dx0*h endif elseif (irep.eq.6) then if (dx0.eq.0.) then p_ice=p_ice0/h else nucleo(2)=dx0/h endif endif open(unit=25,form='formatted',status='unknown',file='cepam.don') write(25,nl_sorties) write(25,nl_espace) write(25,nl_ctes) write(25,nl_temps) write(25,nl_chim) write(25,nl_conv) write(25,nl_etat) write(25,nl_opa) write(25,nl_repartit) write(25,nl_initial) close(unit=25) else write(*,*)' ' open(unit=26,form='formatted',status='unknown', & & file='jacob_cep.tempo') do i=1,7 read(26,226)(a(j,i),j=1,3),(b(j,i),j=1,3),j6_ct write(*,226)(a(j,i),j=1,3),(b(j,i),j=1,3),j6_ct enddo close(26) do i=1,3 do j=1,3 jac(j+(i-1)*3)=(b(j,2*i)-b(j,2*i+1))/(a(i,2*i)-a(i,2*i+1)) j12=(b(j,2*i)-b(j,1))/(a(i,2*i)-a(i,1)) j23=(b(j,1)-b(j,2*i+1))/(a(i,1)-a(i,2*i+1)) jac_er(j+(i-1)*3)=max(abs(j12/(jac(j+(i-1)*3)+1.d-99)-1.d0), & & abs(j23/(jac(j+(i-1)*3)+1.d-99)-1.d0)) write(*,*)i,j,jac(j+(i-1)*3),j12,j23 enddo enddo open(unit=27,form='formatted',status='unknown', & & file='calib.don') write(27,227)(jac(i),jac_er(i),i=1,9) write(*,227)(jac(i),jac_er(i),i=1,9) 227 format('FICHIER CONTENANT LE JACOBIEN POUR CALIBRATION',/, & & t15,'Derivee',t30,'incertitude',/, & & 'dR/dX',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ2/dX',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ4/dX',t12,1pd12.5,t29,1pd10.3,/, & & 'dR/dMc',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ2/dMc',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ4/dMc',t12,1pd12.5,t29,1pd10.3,/, & & 'dR/d**',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ2/d**',t12,1pd12.5,t29,1pd10.3,/, & & 'dJ4/d**',t12,1pd12.5,t29,1pd10.3) if (dx0.eq.0) then write(27,*)'**=pourcentage de glaces' write(*,*)'**=pourcentage de glaces' else write(27,*)'**=discontinuite de X a la PPT' write(*,*)'**=discontinuite de X a la PPT' endif close(27) endif end