module varie integer,save::nlam real,save,allocatable:: wave(:) integer,save::nz=1500 ! stuff connected to use od SSP libraries real, save, allocatable::fluxes(:,:,:),ztab(:), ages(:,:) integer, save, allocatable:: nmax(:) integer, parameter:: maxlib=50 integer, save:: nages=150,ntab,ntab1,imf_flag(maxlib) real,save,allocatable:: varflu(:),prova(:) end module varie ! ******************************************* program leggiSSP_e_interpola !!!passare eta' e metallicita' (ed eventualmente ssp_mass) use varie implicit none real:: ssp_mass,age_SSP,met_SSP integer::jj,ll !N.B: SSPs have Mlow_orig=0.15Mo, Mup=120Mo !for Salpeter, Scalo, Miller e Scalo, Kennicut !ssp_mass_orig respectivlely: 5.0153, 0.09714, 3.4375, 2.9501 !(i.e. it is the mass integral fro Mlow to Mup of the IMF, eg int_0.15^120 M^-1.35 dm) !This quantity is written in the second row of the SSP files. !For small variations of Mlow set the corresponding ssp_mass to rescale the SSP SED: !ssp_mass=IMF mass integral from the required Mlow to Mup write(*,*)'ssp_mass=',5.02,0.097,3.43,2.95,5.81 write(*,*)'respectively for Salpeter, Scalo, Miller&Scalo, Kennicutt (-1.5 slope), and Kennicutt type (-1 slope)' ssp_mass=5.02 !write(*,*) 'now ssp_mass is set to=',ssp_mass write(*,*) 'Give me ssp_mass' read(*,*) ssp_mass call READFLUSSI_asci(ssp_mass) age_SSP=0.01 !!!in Gyr met_SSP=0.02 write(*,*) 'Give me SSP age in Gyr' read(*,*) age_SSP write(*,*) 'Give me SSP metallicity (e.g. 0.02 for solar metallicity)' read(*,*) met_SSP call fluxmetal(age_SSP,met_SSP) !!! a questo punto con use varie o passandoli direttamente !!! qui si hanno wave(nlam)[Ang], varflu(nlam)[1e30 erg/s/Ang/M_sun] !!! dove varflu e' lo spettro della SSP interpolata alla data eta' e metallicita' write(*,*) varflu(1), varflu(nlam) write(*,*) wave(1), wave(nlam) open(15,file='SED_SSP.out') write(15,*) '# wave(Ang) Lum(1e30 erg/s/A/M_sun)' do ll=1,nlam write(15,*) wave(ll),varflu(ll) enddo close(15) end program leggiSSP_e_interpola !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine READFLUSSI_asci(ssp_mass) !,do_all_flag) ! if do_all_flag=.true. then all the subroutine is executed. ! otherwise only the first part in which the SSP section of .par ! is read. This is to know if a double IMF is required when ! reading star formation files ! GL:On 24/01/03 added the possibility to read SSP files in which the wavelenghts ! are directly written as free form flaoting points, to allow high resolution !use modname use varie,only:nlam,wave,nages,fluxes,ztab,ages,ntab,ntab1,nmax,maxlib,imf_flag implicit none real,parameter:: Sunl=3.826e33, Sunl30=3.826e3, Bolsun=4.77 real,parameter:: sigma=5.669e-5,pigreco=3.141593 real,allocatable:: HNU(:),wksp(:) integer,allocatable:: h_int(:),w_int(:),ind(:),iwksp(:) real, allocatable:: norma(:) Character (len=80) TITLE,titleb,fileasc(maxlib), & fileasca,nomfil Character (len=200) filein Character (len=80), allocatable:: swksp(:) real:: teff,agelog,fflum,ssp_mass,ssp_zeta,freq,ssp_mass_orig integer:: jz,jjz,jok,il,ia,err_alloc,nlamold,err real, parameter::safe=1.0 logical::do_all_flag, free_form_wav !ssp_mass_orih e' l' integrale dell' IMF (in massa) vera delle !ssp di cui si leggono i flussi. Questi flussi sono normalizzati !a 1Mo, allora per poter effettuare piccole variazioni della Mlow !dell' IMF si devono denormalizzare i flussi con ssp_mass_orig !e rinormalizzarli con l' integrale relativo alla Mlow desiderata. !N.B: le ssp che abbiamo a disposizione hanno di solito Mlow_orig=0.15Mo, Mup=120Mo !e sono: Salpeter, Scalo, Miller e Scalo, Kennicut, con ssp_mass_orig !risp: 5.0153, 0.09714, 3.4375, 2.9501 open(unit=15,file='per_SSP.dat') read(15,*) read(15,*) jz=0 do jjz=1,maxlib read(15,end=901,err=901,fmt=*)jok,fileasca if(jok.eq.1.or.jok.eq.2) then jz=jz+1 fileasc(jz)=fileasca !imf_flag(jz)=jok endif enddo 901 close(15) ntab=jz ! if (do_all_flag) then allocate(norma(ntab),ztab(ntab),ind(ntab),wksp(ntab),iwksp(ntab), & swksp(ntab),STAT=ERR_ALLOC) IF (ERR_ALLOC /= 0) THEN !call gserror("ALLOCATION ERROR in readflussi") write(*,*) "ALLOCATION ERROR in readflussi" ENDIF norma=0;ztab=0;ind=0;wksp=0;iwksp=0;swksp='' ! first reading of SSP header to ensure everething is OK ! and in particular to arrange fileasc in the correct order ! of metallicies do jz=1,ntab filein=fileasc(jz) open(15,file=filein,status='old',iostat=err) if (err/=0) then print *,'ERROR in opening required SSP file ',filein print *,'May be you are using a wrong-old format of SSP section in .par' stop endif read(15,'(a)') titleb ! information on IMF read(15,*) ssp_mass_orig ! ssp_mass written in ssp if (ssp_mass_orig /= ssp_mass) then norma(jz)=ssp_mass_orig/ssp_mass else norma(jz)=1.0 endif !if (imf_flag(jz)==1) then ! renormalization may apply only to "fundamental" set of SSP! ! norma(jz)=ssp_mass_orig/ssp_mass !else ! norma(jz)=1.0 !endif read(15,*) ztab(jz) ! Zeta of ssp read(15,*) nlam ! number of used wavelenghts if (jz==1) nlamold=nlam if (nlamold/=nlam) then !call gserror('required SSP files have different wavelenghts') write(*,*) 'required SSP files have different wavelenghts' endif nlamold=nlam close(15) enddo if (any(norma(:)/=norma(1)) ) then print *, 'WARNING: you are using SSPs with different IMFs' endif ! arrange relevant arrays in increasing order of metallicity+imf_flag ! in such a way we end up with (possible) different sets of SSP ! each ordered by metallicity !call indexx(ntab,ztab+imf_flag,ind) !wksp(1:ntab)=ztab(1:ntab) ; ztab(1:ntab)=wksp(ind(1:ntab)) !wksp(1:ntab)=norma(1:ntab) ; norma(1:ntab)=wksp(ind(1:ntab)) !iwksp(1:ntab)=imf_flag(1:ntab) ; imf_flag(1:ntab)=iwksp(ind(1:ntab)) !swksp(1:ntab)=fileasc(1:ntab) ; fileasc(1:ntab)=swksp(ind(1:ntab)) ! set ntab1, the index of last ssp with imf_flag=1 !do, jz=1,ntab ! if (imf_flag(jz)==1) ntab1=jz !enddo ntab1=ntab IF(nlam < 0) then nlam=abs(nlam) free_form_wav = .true. else free_form_wav = .false. endif allocate(wave(nlam),hnu(nlam),fluxes(nlam,nages,ntab),& w_int(nlam),h_int(nlam),& ages(nages,ntab),nmax(ntab), STAT=ERR_ALLOC) IF (ERR_ALLOC /= 0) THEN !call gserror( "ALLOCATION ERROR in readflussi") print *, 'ALLOCATION ERROR in readflussi' ENDIF wave=0;hnu=0;fluxes=0;w_int=0;h_int=0;ages=0;nmax=0 ! now full reading of SSP files ! indice ia corre per eta' ! i=1 sta per eta' piu' vecchia mentre ! i=nmax sta per eta' piu' giovane ! tutte le tabelle hanno le stesse dimensioni nmax do jz=1,ntab filein=fileasc(jz) open(15,file=filein,status='old',iostat=err) !now the header can be jumped read(15,'(a)') titleb ! information on IMF read(15,*) read(15,*) read(15,*) ! reading the first library file, need to allocate arrays if(free_form_wav) then !S03 read(15,*) wave !free form else read(15,'(20I6)') W_INT wave(:)=EXP(w_int(:)/10000.0) endif ! read(15,'(20i6)')w_int ! wave(:)=exp(w_int(:)/10000.0) ! transform from nanometer to angstrom wave(:)=wave(:)*10. do ia=1,1000 if(ia .gt. nages) then write(*,*) ' recompile with greater NAGES' stop endif read(15,'(a)',END=90)TITLE read(TITLE,*)TEFF,agelog,fflum ! fflum e' l' integrale indefinito dell' IMF in numero relativo all' ultima ! stella letta nell' isocrona di eta' e metallicita' considerata. ! Invertendo fflum si ottiene la massa ! corrente (non iniziale) della stella (la piu' massiccia della ssp) ages(ia,jz)=agelog ! read the flux in 10**30 ergs/s/hz/ READ(15,'(20i5)') H_int !hnu(:)=exp(h_int(:)/1000.0-100) hnu(:)=h_int(:)/1000.0-100 where (hnu>(-range(hnu(1))+safe)*log(10.0)) hnu=exp(hnu) elsewhere hnu=0 end where ! transform in ergs/sec/Angstrom/ster ! multiply by 4 to restore compatibility with ! planckian: B(T)= integr(Flux) ! remember that this is the way it is saved (nb wave is in angstrom) ! renormalization for different ssp_mass is performed as well fluxes(:,ia,jz)=4*norma(jz)*hnu(:)*2.99792458e18/wave(:)/wave(:) enddo 90 continue nmax(jz)=ia-1 close(15) enddo ages(:,:)=(10**ages(:,:))/1.e9 deallocate(w_int,h_int,hnu,norma,ind,wksp,iwksp,swksp) ! endif ! do_all _flag end Subroutine READFLUSSI_asci !!!!!!!!!!!!!! ! ******************************************************************************************** Subroutine FLUXMETAL(age_SSP,met_SSP) ! Interpola in metallicita' (zzeta) ed eta' (time) i flussi delle SSP ! result is put in variable varflu (defined in prointegrando). ! time e' l' eta', Tg-tbirth, della ssp nata all' istante tbirth. Si ! cercano i valori ztab(i), ztab(i+1) entro cui cade la metallicita' ! z di formazione della ssp; Per entrambe le ztab, si cercano le eta' ! delle ssp tra cui interpolare il flusso (zflu), Poi si interpolano ! questa quantita' tra le due metallicita': varflu(nlam) ! il passaggio da hunt90 e liin con hunt e liin_s per ! il problema con il PGF90 qui rallenta assai su PC. ! Per ora la lascio commentata e vediamo se se ne puo' fare ! a meno su linux! ! use nr_gian, only:liin_s,liin,hunt90 ! use prointegrando ! use remind,only:nlam,wave,nz,icall,tbirth ! use glotauspe, only:fluxes,ztab,ages,ntab,nmax,zzeta,sfr,ntab1, & ! n_bur,t_bur,t_end_bur ! use parametri, only:tgal use varie,only:nlam,wave,nz,fluxes,ztab,ages,ntab,nmax,ntab1,varflu implicit none integer:: err_alloc,il,io,ii,dim,iz1,iz2,it real::t,wk,age,zzeta real, allocatable, save:: zflu(:,:) integer::zero_tab,last_tab ! required to deal with possible use ! of different SSP during bursts real:: age_SSP, met_SSP age=age_SSP !tgal-TBIRTH(ICALL) zzeta=met_SSP ! on first call allocate arrays if (.not.allocated(varflu)) then allocate(zflu(nlam,2),varflu(nlam),STAT=ERR_ALLOC) IF (ERR_ALLOC /= 0) THEN PRINT *, "ALLOCATION ERROR in fluxmetal" STOP ENDIF varflu=0; zflu=0 endif !! here the range of z-indexes of SSP to use is decided ! if ((ntab==ntab1).or.(n_bur==0)) then ! use of different SSP for busts ! ! is NOT required, either because only one SSP set is given or because ! ! there are not bursts. In this second case it is important that ! ! last_tab is set to ntab1, not ntab2 (in the former case is the same) ! zero_tab=0;last_tab=ntab1 ! else ! use of different SSP for busts IS required ! if (any(tbirth(icall)>t_bur(:) .AND. tbirth(icall) ages(1,iz1)) then zflu(:,1)=0 ! set to 0 SSP for ages greater than max available else t=max(age,ages(nmax(iz1),iz1)) ! use youngest SSP available for younger ages call hunt(ages(1:nmax(iz1),iz1),nmax(iz1),t,it) it=max(1,it) ; it=min(it,nmax(iz1)-1) wk=( t-ages(it,iz1) ) / ( ages(it+1,iz1)-ages(it,iz1) ) zflu(:,1)=fluxes(:,it,iz1)*(1-wk)+fluxes(:,it+1,iz1)*wk ! this interpolation in a single instruction over all lambda ! speeds up a lot!! endif if (iz1==iz2) then varflu(:)=zflu(:,1) else if (age > ages(1,iz2)) then zflu(:,1)=0 ! set to 0 SSP for ages greater than max available else t=max(age,ages(nmax(iz2),iz2)) ! use youngest SSP available for younger ages call hunt(ages(1:nmax(iz2),iz2),nmax(iz2),t,it) it=max(1,it) ; it=min(it,nmax(iz2)-1) wk=( t-ages(it,iz2) ) / ( ages(it+1,iz2)-ages(it,iz2) ) zflu(:,2)=fluxes(:,it,iz2)*(1-wk)+fluxes(:,it+1,iz2)*wk endif varflu=zflu(:,1)+(zzeta-ztab(iz1))/(ztab(iz2)-ztab(iz1))*(zflu(:,2)-zflu(:,1)) endif end Subroutine FLUXMETAL ! ******************************************************************************************* SUBROUTINE HUNT(XX,N,X,JLO) DIMENSION XX(N) LOGICAL ASCND ASCND=XX(N).GT.XX(1) IF(JLO.LE.0.OR.JLO.GT.N)THEN JLO=0 JHI=N+1 GO TO 3 ENDIF INC=1 IF(X.GE.XX(JLO).EQV.ASCND)THEN 1 JHI=JLO+INC IF(JHI.GT.N)THEN JHI=N+1 ELSE IF(X.GE.XX(JHI).EQV.ASCND)THEN JLO=JHI INC=INC+INC GO TO 1 ENDIF ELSE JHI=JLO 2 JLO=JHI-INC IF(JLO.LT.1)THEN JLO=0 ELSE IF(X.LT.XX(JLO).EQV.ASCND)THEN JHI=JLO INC=INC+INC GO TO 2 ENDIF ENDIF 3 IF(JHI-JLO.EQ.1)RETURN JM=(JHI+JLO)/2 IF(X.GT.XX(JM).EQV.ASCND)THEN JLO=JM ELSE JHI=JM ENDIF GO TO 3 END SUBROUTINE hunt