program optimise ! Version: 03-Mar-2011 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! optimisation for "CEPAM" c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Author: T.Guillot, Departement Cassiopee, Observatoire de la Cote d'Azur ! Loops on several parameters (p_ice,Yrocks,pd(EOS)) in order to obtain a ! series of optimized models (so that the function "f" that depends on Req, J2, J4...etc.) ! is minimized. ! The results are first put in a .tmp file, then a .opt file. ! The parameters of the minimization are put in the file "cepam_optimise.don" !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Example of a file "cepam_optimise.don": ! &NL_VERSION ! iversion=1/ ! &NL_REFRADJS ! refradjs=6.330d9/ ! &NL_IOPTFCT ! ioptfct=3/ ! &NL_PARAMS ! yn=0.265 0.275, ! yatm=0.233 0.238 0.243, ! zrocks=0.,0.02,0.04, ! pice=0. 1., ! peos=0. 12.5 13.0 13.5 14.0, ! dro=-0.02 0.02/ ! &NL_NPARAMS ! nyn=2, ! nyatm=3, ! nzrocks=3, ! npice=2, ! npeos=2, ! ndro=2/ ! ! ! This provides parameters for the CEPAM/OPTIMISE code ! ! NL_REFRADJS ! ! refradjs is the reference radius for the Js (0 => the boundary CEPAM radius is used instead) ! ! NL_IOPTFCT ! ! ioptfct describes the type of function that is minimized (1-> Req, 2-> Req and J2, 3-> Req, J2 and J4, 31-> ...) ! ! NL_PARAMS ! ! yn: values of the global (protosolar) helium to hydrogen mass mixing ratio ! ! yatm: values of (measured) the atmospheric helium mass mixing ratio ! ! zrocks: mass fraction of rocks in the envelope (molecular+metallic) ! ! pice: fractions of ice in the core (by mass) that are considered ! ! peos: logP (in cgs units) values for modifying the density ! ! dro: fractional density variations that are considered ! ! NL_NPARAMS ! ! nyn...etc.: number of cases that are considered (from the first in the list) [eg. npeos=2 implies that only peos=0. and 12.5 are considered] !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ implicit none include './cepam.inc' !---------------------------------------- ! COMMONS: real*8 req,j2,j4,j6,j8,j10,nu0 common/gravmmts/req,j2,j4,j6,j8,j10,nu0 include './Communs/modelp.inc' include './Communs/planetes.inc' real*8 xpd,xdp,xdro common /nl_densite/xpd,xdp,xdro ! real*8 xrefradjs,xdj2diff(2),xdj4diff(2),xdj6diff(2),xdj8diff(2) ! common /nl_gravopt/xrefradjs,xdj2diff,xdj4diff,xdj6diff,xdj8diff real*8 xrefradjs common /nl_gravopt/xrefradjs integer kioptfct common /nl_igravopt/kioptfct !---------------------------------------- !.... NAMELISTS ..................................... ! iversion is the version of the cepam_optimise.don file integer iversion namelist/NL_VERSION/iversion data iversion/0/ ! refradjs is the reference radius for the Js (0 => the boundary CEPAM radius is used instead) real*8 refradjs namelist/NL_REFRADJS/refradjs data refradjs/0/ ! ioptfct describes the type of function that is minimized (1-> Req, 2-> Req and J2, 3-> Req, J2 and J4, 31-> ...) integer ioptfct namelist/NL_IOPTFCT/ioptfct data ioptfct/2/ ! yn: values of the global (protosolar) helium to hydrogen mass mixing ratio ! yatm: values of (measured) the atmospheric helium mass mixing ratio ! zrocks: mass fraction of rocks in the envelope (molecular+metallic) ! pice: fractions of ice in the core (by mass) that are considered ! eos: logP (in cgs units) values for modifying the density ! dro: fractional density variations that are considered integer pi_pice,pi_yn,pi_yatm,pi_zrocks,pi_peos,pi_dro parameter (pi_pice=2,pi_yn=2,pi_yatm=3,pi_zrocks=3,pi_peos=5, & & pi_dro=2) real*8 yn(pi_yn),yatm(pi_yatm),zrocks(pi_zrocks),pice(pi_pice), & & peos(pi_peos),dro(pi_dro) namelist/NL_PARAMS/yn,yatm,zrocks,pice,peos,dro data yn/0.265,0.275/,yatm/0.233,0.238,0.243/,zrocks/0.,0.02,0.04/,& & pice/0.,1./,peos/0.,12.5,13.0,13.5,14.0/,dro/-0.02,0.02/ ! NL_NPARAMS ! nyn...etc.: number of cases that are considered (from the first in the list) [eg. ieos=2 implies that only eos=0. and 12.5 are considered] integer nyn,nyatm,nzrocks,npice,npeos,ndro namelist/NL_NPARAMS/nyn,nyatm,nzrocks,npice,npeos,ndro data nyn/2/,nyatm/3/,nzrocks/3/,npice/2/,npeos/2/,ndro/2/ !...................................................................... ! Variables locales integer i,iappels,i_pice,i_yn,i_yatm,i_rocks,i_eos,i_dro, & & nx,nn,ireq,ij2,ij4,pp,ierr,nnx real*8 j4min,j4max,xreq,xj2,xj4,f,refj2,refj4,refj6,refj8,refj10 !$$$ real*8 v_pice(pi_pice),v_yn(pi_yn),v_yatm(pi_yatm), & !$$$ & v_rocks(pi_rocks),v_eos(pi_eos),v_dro(pi_dro) !$$$! Jupiter: v_yatm/0.233,0.238,0.243/,v_rocks/0.,0.02,0.04/, !$$$! Saturn: v_yatm/0.18,0.21,0.25/,v_rocks/0.,0.02,0.04/, !$$$ data v_pice/0.,1./,v_yn/0.265,0.275/, & !$$$ & v_yatm/0.18,0.21,0.25/,v_rocks/0.,0.03,0.06/, & !$$$ & v_eos/0.,12.5,13.0,13.5,14.0/,v_dro/-0.02,0.02/ character*1 rep ! Fonction externe integer*4 leng !------------------------------------------------------------------ write(*,100) 100 format(/,'***************** OPTIMIZATION OF A CEPAM MODEL ', & & '*****************',//, & & ' Warning: using files "cepam.don"',/, & & ' and "calib.don"',/, & & ' calculation directly on .bin',/, & & ' results are in file .', & & 'optim',//,'Minimization method:',/, & & 3x,'*** Simplex ***',/) write(*,*)'Name of the to be optimized' read(*,'(a)')modele write(*,101) 101 format('What should we do? (o) Continue from existing loop ', & & '(.tmp model)',/, & & 19x,'(n) start a new loop',/, & & 19x,'(x) no loop: single optimisation.') read(*,'(a)')rep !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture des namelists c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call lit_nl(0) call constantes(fichier_ctes) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Reading file "cepam_optimise.don" ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! open(unit=55,file='cepam_optimise.don',status='old', & & form='formatted',err=50) read(55,NL_VERSION) if (iversion.ne.1) then write(*,*)'optimise: invalid version, iversion=',iversion,'<>1' write(*,*)'check file cepam_optimise.don' stop endif read(55,NL_REFRADJS) read(55,NL_IOPTFCT) read(55,NL_PARAMS) read(55,NL_NPARAMS) close(55) goto 55 50 continue !if there was an error with 'cepam_optimise.don' write(*,*)'optimise: cepam_optimise.don has not been found ', & & 'USING DEFAULT VALUES FOR THE OPTIMISATION' 55 continue xrefradjs=refradjs kioptfct=ioptfct !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Valeurs convenables de J4 c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !$$$ if (fichier_ctes(:11).eq.'data/ctes_j') then !$$$ j4min=-5.89d-4 !$$$ j4max=-5.75d-4 !$$$ v_yatm(1)=0.233 !$$$ v_yatm(2)=0.238 !$$$ v_yatm(3)=0.243 !$$$ elseif (fichier_ctes(:10).eq.'ctes_jup06') then !$$$ j4min=-5.89d-4 !$$$ j4max=-5.75d-4 !$$$ v_yatm(1)=0.233 !$$$ v_yatm(2)=0.238 !$$$ v_yatm(3)=0.243 !$$$ elseif (fichier_ctes(:11).eq.'data/ctes_s') then !$$$ j4min=-9.564d-4 !$$$ j4max=-8.590d-4 !$$$ v_yatm(1)=0.18 !$$$ v_yatm(2)=0.21 !$$$ v_yatm(3)=0.25 !$$$ elseif (fichier_ctes(:10).eq.'ctes_sat06') then !$$$ j4min=-9.386d-4 !$$$ j4max=-9.330d-4 !$$$ v_yatm(1)=0.18 !$$$ v_yatm(2)=0.21 !$$$ v_yatm(3)=0.25 !$$$ else !$$$ j4min=0. !$$$ j4max=0. !$$$ endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! UNE SEULE OPTIMISATION (rep.eq.'x') c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if ((rep.eq.'x').or.(rep.eq.'X')) then write(*,715) 715 format(/,'Optimisation simple a partir des parametres ', & & 'des fichiers .don',/,60('-'),/) !>>>>>>>>>>>>>>>>>>>>>>>>>> call optmod(f) !<<<<<<<<<<<<<<<<<<<<<<<<<< xreq=min(2.d0,abs((req-req_jup)/sig_req)) xj2=min(2.d0,abs((j2-j2jup)/sig_j2)) xj4=min(2.d0,abs((j4-j4jup)/sig_j4)) write(*,716)int(xreq+0.5),int(xj2+0.5),int(xj4+0.5) 716 format(/,60('-'),/,'RESULTATS: IFIT=',i1,i1,i1) stop endif !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! OPTIMISATIONS MULTIPLES c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc nn=0 if ((rep.eq.'o').or.(rep.eq.'O')) then open(unit=12,file='nbpts.tmp',status='old',err=380) read(12,*)nn close(12) write(*,*)'Nombre de points deja calcules: ',nn ! call optmod endif 380 continue !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! On boucle sur les differents parametres c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc nx=0 nnx=nn do i_pice=1,npice p_ice=pice(i_pice) do i_yn=1,nyn yz(1)=yn(i_yn) do i_yatm=1,nyatm yz(2)=yatm(i_yatm) do i_rocks=1,nzrocks yz(5)=zrocks(i_rocks) ! do i_eos=1,pi_eos ! (here pi_eos => 1 => no loop on eos intrinsic variations) do i_eos=1,npeos xpd=peos(i_eos) pp=ndro if (i_eos.eq.1) pp=1 !pour eviter de calculer 2x avec Pd=0 do i_dro=pp,1,-1 !pour modif de densite positif d'abord ! do i_dro=1,pp ! do i_dro=1,1 xdro=dro(i_dro) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! APPEL POUR OPTIMISATION c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc nx=nx+1 if (nx.gt.nn) then if (nx.eq.nn+1) then ! On reinitialise tous les parametres (modifies par lecture des .don) write(*,555)nx 555 format(//,70('-'),/, & & 'ON REINITIALISE LES PARAMETRES; nx=',i4,/,70('-')) call cepam(1,ierr) p_ice=pice(i_pice) yz(1)=yn(i_yn) yz(2)=yatm(i_yatm) yz(5)=zrocks(i_rocks) xpd=peos(i_eos) xdro=dro(i_dro) if (ierr.gt.0) then write(*,*)'pb: on reappelle CEPAM /ierr=',ierr call cepam(1,ierr) if (ierr.ne.0) stop endif endif iappels=0 666 iappels=iappels+1 !>>>>>>>>>>>>>>>>>>>>>>>>>> call optmod(f) !<<<<<<<<<<<<<<<<<<<<<<<<<< ! call cepam(3,ierr) write(*,*)'OPTIMISE: Req, Nu0=',req,nu0 if (ierr.ne.0) stop ! Normalizing the Js to the reference radius (if there is one) if (xrefradjs.gt.0) then refj2=(req/xrefradjs)**2*j2 refj4=(req/xrefradjs)**4*j4 refj6=(req/xrefradjs)**6*j6 refj8=(req/xrefradjs)**8*j8 refj10=(req/xrefradjs)**10*j10 else refj2=j2 refj4=j4 refj6=j6 refj8=j8 refj10=j10 endif xreq=min(9.d0,abs((req-req_jup)/sig_req)) xj2=min(9.d0,abs((refj2-j2jup)/sig_j2)) xj4=min(9.d0,abs((refj4-j4jup)/sig_j4)) ! This part added to be able to write *all steps* (in sub_sep) if (nnx.gt.0) then open(unit=12,file='nbpts.tmp',status='unknown') read(12,*)nnx close(12) endif nnx=nnx+1 open(unit=11,file=modele(:leng(modele))//'.tmp', & & status='unknown') do i=1,nnx-1 read(11,*) enddo write(11,111)xpd,xdp,xdro,p_ice,yz(1),yz(2),yz(5), & & yz(3),mnoyau,req,refj2,refj4,refj6,refj8, & & f,int(xreq),int(xj2),int(xj4),nu0*1d6 write(*,112)xpd,xdp,xdro,p_ice,yz(1),yz(2),yz(5), & & yz(3),mnoyau,req,refj2,refj4,refj6,refj8, & & f,int(xreq),int(xj2),int(xj4),nu0*1d6 111 format(3f7.3,4f7.4,2f11.8,1p,d13.6,4d14.6,d9.2,1x,i1,i1,i1,& & 0p,f10.4) 112 format(/,3(72('#'),/), & & 3f7.3,4f7.4,2f11.8,1p,d13.6,4d14.6,d9.2,1x,i1,i1,i1, & & 0p,f10.4,/,3(72('#'),/),/) 113 format(t2,'logPd',t9,'dlogP',t16,'dRho',t23,'Pice',t30, & & t29,'YN/(X+Y)',t37,'Ymol',t44,'Zrocks',t53, & & 'Zices',t64,'Mcore',t75,'Req',t88,'J2',t102,'J4', & & t116,'J6',t130,'J8',t143,'F(min)',t151,'Fit?', & & t156,'Nu0') close(11) open(unit=12,file='nbpts.tmp',status='unknown') write(12,*)nnx close(12) else write(*,*)'Point: ',nx,' ... passe' endif enddo enddo enddo enddo enddo enddo close(11) open(unit=11,file=modele(:leng(modele))//'.tmp',status='old') open(unit=13,file=modele(:leng(modele))//'.opt',status='unknown') write(13,113) write(13,*)nx do i=1,nx read(11,111)xpd,xdp,xdro,p_ice,yz(1),yz(2),yz(5), & & yz(3),mnoyau,req,j2,j4,j6,j8,f,ireq,ij2,ij4,nu0 write(13,111)xpd,xdp,xdro,p_ice,yz(1),yz(2),yz(5), & & yz(3),mnoyau,req,j2,j4,j6,j8,f,ireq,ij2,ij4,nu0 enddo close(13) write(*,*)'J''ai cree le fichier ',modele(:leng(modele))//'.opt' end