Module READ_MODEL ! Implicit None ! Private Public :: EXTRACT_ASCII, CALCRHOX, CALC_ERROR ! CONTAINS ! Subroutine EXTRACT_ASCII(fin, Teff, grav, metal, mass, ndepth, + xlr_ref, tau5, tauR, temp, prese, presg, + xit, rad, sph, xkapr, ro, abund) ! ! Adapted from P. DeLaverny ! Character(len = *), intent(in) :: fin Integer, intent(out) :: ndepth Real, intent(out) :: Teff, grav, metal, xlr_ref Real, dimension(:), intent(out) :: tau5, tauR, temp, prese, presg, + xit, rad, xkapr, ro Logical, intent(out) :: sph Real, dimension(:), intent(out) :: abund ! Integer :: i, k, idum, ios Integer, parameter :: ndp = 200, imod = 10 Character(len = 50) :: modcode, blabla Logical :: err Real :: radius, mass, GG, xic Real, dimension(ndp) :: tau, t, z, ptot, prad, pg, pturb, pe, + xmass, grad, gravity, pcheck, rr, + geff, gradptur, dp, taus, coldens Real, dimension(30) :: xlr real, dimension(ndp) :: emu, vconv, fconv Real, dimension(155000) :: xlb, w, fluxme ! Inquire(file = fin, exist = err) ! Write(*, *) "fin = ", Trim(fin(110 : )) ! If (.NOT. err) Then STOP "Extract_ascii : PB-00 : UNKNOWN FILE." End If ! Open(imod, file = fin, status = 'OLD', iostat = ios) IF (ios /= 0) STOP "Extract_ascii : PB-01" Read(imod, '(a)', iostat = ios) modcode IF (ios /= 0) STOP "Extract_ascii : PB-02" ! If (modcode(1 : 1) .eq. 'p' .or. modcode(1 : 3) .eq. 'sun') Then Write(*, *) ' This model is PLANE PARALLEL' Else If (modcode(1 : 1) .eq. 's') Then Write(*, *) ' This model is SPHERICALLY SYMMETRIC' Else Write(*, *) ' This model may not be a NewMARCS model!' End If ! sph = (modcode(1 : 1) .eq. 's') xlr_ref = 5000 Read(imod, *, iostat = ios) Teff If (ios .ne. 0) STOP "Extract_ascii : PB-03" Read(imod, *) Read(imod, *, iostat = ios) grav If (ios .ne. 0) STOP "Extract_ascii : PB-04" grav = log10(grav) Read(imod, *, iostat = ios) xic If (ios .ne. 0) STOP "Extract_ascii : PB-05" Read(imod, *, iostat = ios) mass If (ios .ne. 0) STOP "Extract_ascii : PB-06" Read(imod, *, iostat = ios) metal If (ios .ne. 0) STOP "Extract_ascii : PB-07" Read(imod, *, iostat = ios) radius If (ios .ne. 0) STOP "Extract_ascii : PB-08" ! blabla = "" ! Do While (Index(blabla,'abundances') .eq. 0) Read(imod,'(a)', iostat = ios) blabla If (ios .ne. 0) STOP "Extract_ascii : PB-13" End Do ! Do k = 1, 9 Read(imod,'(10F7.2)', iostat = ios)(abund(10*(k-1)+i), i=1, 10) If (ios .ne. 0) STOP "Extract_ascii : PB-14" End Do ! Read(imod,'(2F7.2)', iostat = ios) abund(91), abund(92) If (ios .ne. 0) STOP "Extract_ascii : PB-15" ! Do While (blabla .ne. 'Model structure') Read(imod,'(a)', iostat = ios) blabla If (ios .ne. 0) STOP "Extract_ascii : PB-09" End Do ! Backspace(imod) Backspace(imod) Read(imod, *, iostat = ios) ndepth If (ios .ne. 0) STOP "Extract_ascii : PB-10" Read(imod, *) Read(imod, *) ! Do k = 1, ndepth read(imod, *, iostat = ios) idum, tauR(k), tau5(k), rad(k), & temp(k), Pe(k), Pg(k) If (ios .ne. 0) STOP "Extract_ascii : PB-11" prese(k) = log10(Pe(k)) presg(k) = log10(Pg(k)) xit(k) = xic rad(k) = radius - rad(k) End Do ! Read(imod,*) ! Do k = 1, ndepth Read(imod,*, iostat = ios) idum, tauR(k), xkapr(k), ro(k), & emu(k), Vconv(k), Fconv(k) If (ios .ne. 0) STOP "Extract_ascii : PB-12" xkapr(k) = log10(xkapr(k)) End Do ! Close(imod) ! End Subroutine EXTRACT_ASCII ! !----------------------------------------------------------------- ! Subroutine CALC_ERROR(xinf, xsup, yinf, ysup, zinf, zsup, + teff_ref, logg_ref, z_ref) ! Real, intent(in) :: xinf, xsup, yinf, ysup, zinf, zsup, + teff_ref,logg_ref,z_ref ! real :: error_T, error_Pe, error_Pg, error_kappa real :: errorTeffT, errorloggT, errorzT, + errorTeffPe, errorloggPe, errorzPe, + errorTeffPg, errorloggPg, errorzPg, + errorTeffkappa, errorloggkappa, errorzkappa ! ! Values read out of the figures of the manual and ! scaled down o the according step ! errorTeffT = 0.055 / 32 errorloggT = 0.008 / 5 errorzT = 0.015 / 4 errorTeffPe = 0.65 / 32 errorloggPe = 0.4 / 5 errorzPe = 0.38 / 4 errorTeffPg = 0.25 / 32 errorloggPg = 0.23 / 5 errorzPg = 0.35 / 4 errorTeffkappa = 0.8 / 32 errorloggkappa = 0.36 / 5 errorzkappa = 0.38 / 4 ! error_T = min((xsup -teff_ref), (teff_ref-xinf))/100*errorTeffT + + min((ysup-logg_ref), (logg_ref-yinf))*errorloggT + + min((zsup-z_ref), (z_ref-zinf))*errorzT Write(*, 1409) 'Estimated max error on T =', error_T*100, '%' ! error_Pe = min((xsup-teff_ref), (teff_ref-xinf))/100*errorTeffPe + + min((ysup-logg_ref), (logg_ref-yinf))*errorloggPe + + min((zsup-z_ref), (z_ref-zinf))*errorzPe Write(*, 1409) 'Estimated max error on Pe =', error_Pe*100, '%' ! error_Pg= min((xsup-teff_ref), (teff_ref-xinf))/100*errorTeffPg + + min((ysup-logg_ref), (logg_ref-yinf))*errorloggPg + + min((zsup-z_ref), (z_ref-zinf))*errorzPg Write(*, 1409) 'Estimated max error on Pg =', error_Pg*100, '%' ! error_kappa = min((xsup-teff_ref), (teff_ref-xinf))/100 + *errorTeffkappa + + min((ysup-logg_ref), (logg_ref-yinf))*errorloggkappa + + min((zsup-z_ref), (z_ref-zinf))*errorzkappa Write(*, 1409) 'Estimated max error on kappa =', error_kappa*100, + '%' ! 1409 format(a30, f5.1, a2) ! End Subroutine CALC_ERROR ! !----------------------------------------------------------------- ! Subroutine CALCRHOX(tau, kappa, ndepth, rhox) ! ! 2 ways to calculate rhox : int(ro*dx) or int(1/kappa*dtau) ! A&A 387, 595-604 (2002) ! Integer, intent(in) :: ndepth Real, dimension(:), intent(in) :: tau, kappa Real, dimension(:), intent(out) :: rhox ! Integer :: I real :: first, tot real, dimension(ndepth) :: f, x real, external :: rinteg ! f = (1/10**kappa) x = 10**(tau) first = x(1)*f(1) tot = rinteg(x, f, rhox, ndepth, first) ! Do i = 2, ndepth rhox(i) = rhox(i-1) + rhox(i) End Do ! End Subroutine CALCRHOX ! End Module READ_MODEL ! !***************************************************************** ! Program interpol_modeles ! !***************************************************************** ! interpolation of model atmosphere ! parameter space of interpolation {Teff,logg,z} ! 8 MARCS binary models as input required: ! ! the order of input models matters ! ! ! only standard composition has been tested, ! ! no guaranty for peculiar composition ! ! turbospectrum/babsma format compatible output ! ! Interpolation scheme of model atmosphere(@) ! in {Teff,logg,z} stellar parameter space ! ^ z ! |_ _ _ _ _ ! | /| ! /| / ! _ |_ _ _ / | ! | | | ! | | ! | | @ | ! |--------------> logg ! |/ | / ! /_ _ _ _ ./ ! / ! / ! Teff ! ! Each structure component of each model ! (T,Pe,Pg,xit,optical depth, kappaross) ! is first resampled on a common optical depth basis ! tau(5000 or ross) (see resample routine). ! ! Then each component of the atmosphere (@) ! (T,Pe,Pg,xit,kappa,geom depth +tau5000 and tauross)) ! is interpolated individually along each dimension between ! each input model value(*) ! (physically it is impossible to ensure a simple relationship between models in more than ! one dimension at the same time, that's why the input models MUST form a "cube" ! in the stellar parameter space {Teff,logg,z}). ! The interpolation is successively done at each optical depth (tau) ! + interpolation point weighted by an empirical value (see TM thesis + manual) ! ! ! ^ T,Pe,Pg,xit or radius ! | ! | * ! | /| ! | @ ! | /| | ! | / / ! | / / | | ! ^ | / * ! | |/ | | | ! | -----------------------> Teff, logg and z ! | / low ref up ! | /* / / / ! | / |\ ! | / \ / / / ! | / | \ ! | / /@ / / ! |/ | |\ ! / _ _ /_ _ /*_ _ _ _ _ ! / low ref up ! / ! tau{5000 or Ross} ! !***************************************************************** ! ! TM 07/2003 ! ! 07/2004 resampling of each model on a common optical depth basis ! 06/2006 works for spherical geometry models ! 09/2007 new calibration of free parameter alpha ! + modified to read Uppsala ascii models ! + kappa interpolation ! + rhox calculation ! + 2 outputs (babsma and ATLAS/MOOG) ! 10/2007 error estimates ! 07/2010 Conversion in fortran95 (Thibault Merle) ! + interpolated model in MOOG/MULTI format !***************************************************************** ! ! pathf95 -C -Wall -o INTERMARCS interpol_modeles.f ! USE READ_MODEL, ONLY : EXTRACT_ASCII, CALCRHOX, CALC_ERROR ! Implicit None ! Integer :: ios Integer :: f_in, k, ndepth_ref, out, nlinemod Integer, parameter :: ndp = 200, nfile = 10 Logical :: verif, check, test, err Logical :: extrapol, binary,optimize Real :: lambda_ref, temp_ref, logg_ref, z_ref Character(len=2) :: comp Real :: radius_interp, mass_interp Real :: x, y, z Real :: xinf, xsup, yinf, ysup, zinf, zsup Real :: D_z, D_a Real :: teffpoint, loggpoint, metpoint Character(len = 256), dimension (nfile) :: FILE_IN Character(len = 256) :: FILE_ALT Real, dimension (:, :), allocatable :: taus, tauR, T, Pe, Pg, + xit, rr, xkapref,rhox, ro Real, dimension (:, :), allocatable :: taus_aux, tauR_aux, T_aux, + Pe_aux, Pg_aux, xit_aux, + rr_aux, xkapref_aux,ro_aux, + abund Integer, dimension (nfile) :: ndepth Real, dimension (nfile) :: xlr, teff, logg, metal, mass Logical, dimension (nfile) :: sph External :: blend_103 Real, external :: inf, sup Real, dimension(9, 3) :: lin_dif, power Character(len = 256) :: fout ! Solar abundances Real, dimension(16), parameter :: solabu = (/12.00, 10.93, 7.51, + 7.53, 6.37, 7.45, 8.39, 6.17, 7.14, 5.08, 6.31, 6.23, 5.64, 7.78, + 8.66, 7.84/) ! ! Option for MULTI format output Character(len = 10), parameter :: MODELTYPE = "TAU5000" ! TAU5000 OR MASS Logical, parameter :: LOGASCALE = .TRUE. ! FOR ABUNDANCES Logical :: type1, type2, type3, type4, type5 ! INTERFACE reec subroutine resample(taus, tauR, T, Pe, Pg, xit, rr, xkapref) real, dimension(:, :) :: taus, tauR, T, Pe, Pg, xit, rr, xkapref end END INTERFACE reec ! Write(*, *) '******************************' Write(*, *) '* Start of Interpolation *' Write(*, *) '******************************' ! ! You can choose here to switch off the "optimization" and ! prefer simple linear interpolation ! optimize = .TRUE. ! ! Read 8 models, put in tables, ! Check number of layers, reference optical depth, ! and geometry compatibility ! out = 9 ! Write(*, *) 'Interpolation between :' ! Do f_in = 1, 9 Read(*, *, iostat = ios) FILE_IN(f_in) IF (ios /= 0) STOP "Main : PB-00" End Do ! Read(*, *) ! Added to ensure compatibility with interpol_models from bacchus. Read(*, *, iostat = ios) temp_ref IF (ios /= 0) STOP "Main : PB-01" read(*, *, iostat = ios) logg_ref IF (ios /= 0) STOP "Main : PB-02" read(*, *, iostat = ios) z_ref IF (ios /= 0) STOP "Main : PB-03" read(*, *, iostat = ios) comp IF (ios /= 0) STOP "Main : PB-031" read(*, *, iostat = ios) test IF (ios /= 0) STOP "Main : PB-04" ! verif = .TRUE. check = .TRUE. nlinemod = ndp ! Allocate(taus_aux(nlinemod, nfile), tauR_aux(nlinemod, nfile), + T_aux(nlinemod, nfile), Pe_aux(nlinemod, nfile), + Pg_aux(nlinemod, nfile), xit_aux(nlinemod, nfile), + rr_aux(nlinemod, nfile), xkapref_aux(nlinemod, nfile), + ro_aux(nlinemod, nfile), abund( 92, nfile), + stat = ios) IF (ios /= 0) STOP "Main : PB-05" ! Read(*, *, iostat = ios) binary IF (ios /= 0) STOP "Main : PB-06" ! If (binary) Then ! Do f_in = 1, 8 CALL extract_bin(FILE_IN(f_in), + teff(f_in), logg(f_in), metal(f_in), + ndepth(f_in), xlr(f_in), taus_aux(:, f_in), + tauR_aux(:, f_in), T_aux(:, f_in), + Pe_aux(:, f_in), Pg_aux(:, f_in), + xit_aux(:, f_in), rr_aux(:, f_in), + sph(f_in), xkapref_aux(:, f_in)) verif = verif .AND. (ndepth(f_in) == ndepth(1)) check = check .AND. (xlr(f_in) == xlr(1)) If ( .NOT. (((sph(1) .AND. sph(f_in))) .OR. + ((.NOT. (sph(1))) .AND. (.NOT. (sph(f_in)))))) Then Write(*, *) 'Geometry compatibility problem with' Write(*, 78) f_in, teff(f_in), logg(f_in), metal(f_in) STOP End If Write(*, 78) f_in, teff(f_in), logg(f_in), metal(f_in) End Do ! Else ! Do f_in = 1, 8 CALL extract_ascii(FILE_IN(f_in), + teff(f_in), logg(f_in), metal(f_in), + mass(f_in), + ndepth(f_in), xlr(f_in), taus_aux(:, f_in), + tauR_aux(:, f_in), T_aux(:, f_in), + Pe_aux(:, f_in), Pg_aux(:, f_in), + xit_aux(:, f_in), rr_aux(:,f_in), + sph(f_in), xkapref_aux(:, f_in), + ro_aux(:, f_in), abund(:, f_in)) verif = verif .AND. (ndepth(f_in) == ndepth(1)) check = check .AND. (xlr(f_in) == xlr(1)) If ( .NOT. (((sph(1) .AND. sph(f_in))) .OR. + ((.NOT. (sph(1))) .AND. (.NOT. (sph(f_in)))))) Then Write(*, *) 'Geometry compatibility problem with' Write(*, 78) f_in, teff(f_in), logg(f_in), metal(f_in) STOP End If Write(*, 78) f_in, teff(f_in), logg(f_in), metal(f_in) End Do ! End If ! 78 Format('Model', i2,' Teff =', f8.0,' logg =', f5.2,' z =', f6.2) ! Write(*, 79) temp_ref, logg_ref, z_ref 79 Format('Interpolation point : Teff =', f8.0,' logg =', f5.2, + ' z =', f6.2) ! ! Check if files are length and depth ref compatible ! If (.NOT.(check)) Then write(*, *) 'All the models do not have the same lambda ref.' write(*, *) 'No interpolation done.' STOP Else If (.NOT.(verif)) Then Write(*, *) "WARNING: All the models don't have the same number" Write(*, *) "of layers, resampling to", ndepth(1), "layers." End If ! ndepth_ref = ndepth(1) lambda_ref = xlr(1) ! ! Check if the models have the same mass ! If( All(mass(1:8)/= Minval(mass(1:8))) ) Then STOP "Main: PB-28" Else mass_interp = Minval(mass(1:8)) End If ! ! Calculation of the interpolation point(x, y, z) in {Teff, logg, z} space ! Allocate(taus(ndepth_ref, nfile), tauR(ndepth_ref, nfile), + T(ndepth_ref, nfile), Pe(ndepth_ref, nfile), + Pg(ndepth_ref, nfile), xit(ndepth_ref, nfile), + rr(ndepth_ref, nfile), xkapref(ndepth_ref, nfile), + ro(ndepth_ref, nfile), + stat = ios) IF (ios /= 0) STOP "Main : PB-07" ! taus = taus_aux(1 : ndepth_ref, :) tauR = tauR_aux(1 : ndepth_ref, :) T = T_aux(1 : ndepth_ref, :) Pe = Pe_aux(1 : ndepth_ref, :) Pg = Pg_aux(1 : ndepth_ref, :) xit = xit_aux(1 : ndepth_ref, :) rr = rr_aux(1 : ndepth_ref, :) xkapref = xkapref_aux(1 : ndepth_ref, :) ro = ro_aux(1 : ndepth_ref, :) ! xinf = minval(teff(1 : 8)) yinf = minval(logg(1 : 8)) zinf = minval(metal(1 : 8)) xsup = maxval(teff(1 : 8)) ysup = maxval(logg(1 : 8)) zsup = maxval(metal(1 : 8)) ! If (xsup == xinf) Then teffpoint = 0 Else teffpoint = (temp_ref-xinf) / (xsup-xinf) End If ! If (ysup == yinf) Then loggpoint = 0 Else loggpoint = (logg_ref-yinf) / (ysup-yinf) End if ! If (zsup == zinf) Then metpoint = 0 Else metpoint=(z_ref-zinf)/(zsup-zinf) End If ! extrapol = ((teffpoint < 0) .OR. (teffpoint > 1) .OR. + (loggpoint < 0) .OR. (loggpoint > 1) .OR. + (metpoint < 0) .OR. (metpoint > 1)) ! If (extrapol) Then Write(*, *) '!!! WARNING : extrapolation !!!' End If ! ! Resample each layer of each input model on a common depth basis ! (tau5000 or tauRoss, see resample routine) ! If you don't want to resample all the model to the same depth scale, ! just comment the following line ! CALL resample(taus, tauR, T, Pe, Pg, xit, rr, xkapref) ! ! Initialisation of empirical constants for optimized interpolation ! (see TM thesis) ! lin_dif(1, 1) = 0 !tau5000 vs Teff lin_dif(1, 2) = 0 ! ... vs logg lin_dif(1, 3) = 0 ! ... vs z lin_dif(2, 1) = 0 !tauross vs Teff lin_dif(2, 2) = 0 ! ... vs logg lin_dif(2, 3) = 0 ! ... vs z lin_dif(3, 1) = 0.15 !T vs Teff lin_dif(3, 2) = 0.3 ! ... vs logg lin_dif(3, 3) = 1-(temp_ref/4000)**2.0 ! ... vs z lin_dif(4, 1) = 0.15 !logPe vs Teff lin_dif(4, 2) = 0.06 ! ... vs logg lin_dif(4, 3) = 1-(temp_ref/3500)**2.5 ! ... vs z lin_dif(5, 1) = -0.4 !logPg vs Teff lin_dif(5, 2) = 0.06 ! ... vs logg lin_dif(5, 3) = 1-(temp_ref/4100)**4 ! ... vs z lin_dif(6, 1) = 0 !xit vs Teff lin_dif(6, 2) = 0 ! ... vs logg lin_dif(6, 3) = 0 ! ... vs z lin_dif(7, 1) = 0 !rr vs Teff lin_dif(7, 2) = 0 ! ... vs logg lin_dif(7, 3) = 0 ! ... vs z lin_dif(8, 1) = -0.15 !logxkapref vs Teff lin_dif(8, 2) = -0.12 ! ... vs logg lin_dif(8, 3) = 1-(temp_ref/3700)**3.5 ! ... vs z ! !ADD TM lin_dif(9, 1) = 0 !ro vs Teff lin_dif(9, 2) = 0 ! ... vs logg lin_dif(9, 3) = 0 ! ... vs z !END ADD TM ! If (optimize) Then Write(*, fmt = '(A)', advance = 'no') ' Optimized interpolation' Write(*, *) ' applied for standard composition models.' Else lin_dif = 0. Write(*,*) 'Linear interpolation applied.' End if ! ! These constants are calibrated on a broad range of stellar parameters; ! scale them now to the present one. ! power(:, 1) = 1 - (lin_dif(:, 1) * (abs(xsup-xinf)/(7000-3800))) power(:, 2) = 1 - (lin_dif(:, 2) * (abs(ysup-yinf)/(5-0.0))) power(:, 3) = 1 - (lin_dif(:, 3) * (abs(zsup-zinf)/(0-(-4)))) ! ! Interpolation of each component of the atmosphere ! (taus, teff, Pe, Pg, microt, rr) and at each layer. ! Do k = 1, ndepth_ref x = teffpoint**power(1, 1) y = loggpoint**power(1, 2) z = metpoint**power(1, 3) CALL blend_103(x, y, z, + taus(k, 1), taus(k, 2), taus(k, 3), taus(k, 4), + taus(k, 5), taus(k, 6), taus(k, 7), taus(k, 8), + taus(k, out)) x = teffpoint**power(2, 1) y = loggpoint**power(2, 2) z = metpoint**power(2, 3) CALL blend_103(x, y, z, + tauR(k, 1), tauR(k, 2), tauR(k, 3), tauR(k, 4), + tauR(k, 5), tauR(k, 6), tauR(k, 7), tauR(k, 8), + tauR(k, out)) x = teffpoint**power(3,1) y = loggpoint**power(3,2) z = metpoint**power(3,3) CALL blend_103(x, y, z, + T(k, 1), T(k, 2), T(k, 3), T(k, 4), + T(k, 5), T(k, 6), T(k, 7), T(k, 8), + T(k, out)) x = teffpoint**power(4,1) y = loggpoint**power(4,2) z = metpoint**power(4,3) CALL blend_103(x, y, z, + Pe(k, 1), Pe(k, 2), Pe(k, 3), Pe(k, 4), + Pe(k, 5), Pe(k, 6), Pe(k, 7), Pe(k, 8), + Pe(k, out)) x = teffpoint**power(5,1) y = loggpoint**power(5,2) z = metpoint**power(5,3) CALL blend_103(x, y, z, + Pg(k, 1), Pg(k, 2), Pg(k, 3), Pg(k, 4), + Pg(k, 5), Pg(k, 6), Pg(k, 7), Pg(k, 8), + Pg(k, out)) x = teffpoint**power(6,1) y = loggpoint**power(6,2) z = metpoint**power(6,3) CALL blend_103(x, y, z, + xit(k, 1), xit(k, 2), xit(k, 3), xit(k, 4), + xit(k, 5), xit(k, 6), xit(k, 7), xit(k, 8), + xit(k, out)) x = teffpoint**power(7,1) y = loggpoint**power(7,2) z = metpoint**power(7,3) CALL blend_103(x, y, z, + rr(k, 1), rr(k, 2), rr(k, 3), rr(k, 4), + rr(k, 5), rr(k, 6), rr(k, 7), rr(k, 8), + rr(k, out)) x = teffpoint**power(8,1) y = loggpoint**power(8,2) z = metpoint**power(8,3) CALL blend_103(x, y, z, + xkapref(k, 1), xkapref(k, 2), xkapref(k, 3), xkapref(k, 4), + xkapref(k, 5), xkapref(k, 6), xkapref(k, 7), xkapref(k, 8), + xkapref(k, out)) x = teffpoint**power(9,1) y = loggpoint**power(9,2) z = metpoint**power(9,3) !ADD TM CALL blend_103(x, y, z, + ro(k, 1), ro(k, 2), ro(k, 3), ro(k, 4), + ro(k, 5), ro(k, 6), ro(k, 7), ro(k, 8), + ro(k, out)) ! END ADD TM End Do ! ndepth(out) = ndepth_ref xlr(out) = lambda_ref ! ! Now calculate rhox ! Write(*, *) 'Now calculate rhox.' ! Allocate(rhox(ndepth_ref, nfile), stat = ios) IF (ios /= 0) STOP "Main : PB-08" ! Do f_in = 1, 9 CALL calcrhox(tauR(:, f_in), xkapref(:, f_in), + ndepth_ref, rhox(:, f_in)) End Do ! ! Calculate estimated error ! If (optimize) Then Write(*, *) 'Now calculate error.' CALL calc_error(xinf, xsup, yinf, ysup, zinf, zsup, + temp_ref, logg_ref, z_ref) End If ! ! Compute the stellar radius of interpolated model (in sph case) ! radius_interp = sqrt(mass_interp*10**(4.44-logg_ref)) ! ! Identification of format output ! type1 = Index(FILE_IN(9), 'TURBO') /= 0 type2 = Index(FILE_IN(9), 'ATLAS') /= 0 type3 = Index(FILE_IN(9), 'MOOG') /= 0 type4 = Index(FILE_IN(9), 'MULTI') /= 0 ! ! Write interpolated model in file nber out ! Write(*, *) 'Now write result.' ! If (type1) Then Call WRITE_TURBO() Else If (type2) Then Call WRITE_ATLAS() Else If (type3) Then Call WRITE_MOOG() Else If (type4) Then Call WRITE_MULTI() Else If (type5) Then Call WRITE_COMP() ! Complete model Else STOP "Main : FORMAT OUTPUT UNKNOWN" End If ! ! Write a file compatible for sm (used for a control plot) ! Open (unit = 24, file = 'model_atmosphere.sm', iostat = ios) IF (ios /= 0) STOP "Main : PB-24" Do f_in = 1, 8 Do k = 1, ndepth_ref Write(24, *, iostat = ios) taus(k,f_in),T(k,f_in),Pe(k,f_in), + Pg(k,f_in),taur(k,f_in) IF (ios /= 0) STOP "Main : PB-25" End Do End Do ! ! Case of a 10th comparison model ! Read(*, *, iostat = ios) FILE_IN(10) IF (ios /= 0) STOP "Main : PB-26" ! Inquire(file = FILE_IN(10), exist = err) If (.NOT. err) Then test = .FALSE. Write(*, *) "The file ", Trim(FILE_IN(10)), " doesn't exist." End If ! If (test) Then f_in = 10 If (binary) Then CALL extract_bin(FILE_IN(f_in), + teff(f_in), logg(f_in), metal(f_in), + ndepth(f_in), xlr(f_in), taus_aux(:, f_in), + tauR_aux(:, f_in), T_aux(:, f_in), + Pe_aux(:, f_in), Pg_aux(:, f_in), + xit_aux(:, f_in), rr_aux(:,f_in), + sph(f_in), xkapref_aux(:,f_in)) Else CALL extract_ascii(FILE_IN(f_in), + teff(f_in), logg(f_in), metal(f_in), + mass(f_in), + ndepth(f_in), xlr(f_in), taus_aux(:, f_in), + tauR_aux(:, f_in), T_aux(:, f_in), + Pe_aux(:, f_in), Pg_aux(:, f_in), + xit_aux(:, f_in),rr_aux(:, f_in), + sph(f_in), xkapref_aux(:, f_in), + ro_aux(:,f_in), abund(:, f_in)) End If ! Do k = 1, ndepth(f_in) Write(24, *, iostat = ios) taus_aux(k, f_in), T_aux(k, f_in), + Pe_aux(k, f_in), Pg_aux(k, f_in), + tauR_aux(k, f_in) IF (ios /= 0) STOP "Main : PB-27" End Do ! End If ! Close(24) ! If (extrapol) Then Write (*, *) 'Extrapolation done !' Else write (*, *) 'Interpolation done !' End If ! Deallocate(taus, tauR, T, Pe, Pg, xit, taus_aux, tauR_aux, T_aux) Deallocate(Pe_aux, Pg_aux, xit_aux, rr_aux, rr, rhox) ! CONTAINS ! !----------------------------------------------------------------- ! Subroutine WRITE_COMP() ! Open(26, file = FILE_IN(9), iostat = ios) IF (ios /= 0) STOP "WRITE_COMP: PB-00." Write(26, *) "k log(tau5) log(tauR) kappaR z m T Ne rho" Write(26, *) " [cm2/g] [cm] [g/cm2] [K] [1/cm3] [g/cm3]" ! Do k = 1, ndepth_ref Write(26, '(I2, 8ES12.4)') k, taus(k, out), tauR(k, out), + xkapref(k, out), rr(k, out), + rhox(k, out), T(k ,out), + 10**(Pe(k, out)) / (1.38065E-16 * T(k, out)), ro(k, out) End Do ! Close(26) End Subroutine WRITE_COMP ! !----------------------------------------------------------------- ! Subroutine WRITE_MULTI() ! ! atmos.marcs ! Open(26, file = "atmos." // FILE_IN(9), iostat = ios) IF (ios /= 0) STOP "WRITE_MULTI: PB-00." Write(26, '(A)') "MARCS" ! If (Trim(MODELTYPE) == "TAU5000") Then Write(26, '(A)') "TAU5000 SCALE" Else If (Trim(MODELTYPE) == "MASS") Then Write(26, '(A)') "MASS SCALE" Else STOP"WRITE_MULTI: MODELTYPE UNKNOWN, ONLY 'TAU5000' OR 'MASS'" End If ! Write(26, '(A)') "*" Write(26, '(A)') "* INTERPOLATED MODEL BY MASSERON CODE + (TM VERSION)" Write(26, '(A, F7.1, A, F5.2, A, F5.2)') + "* Teff_min =", xinf, " loggf_min =", yinf, " [Fe/H]min =", zinf Write(26, '(A, F7.1, A, F5.2, A, F5.2)') + "* Teff_max =", xsup, " loggf_max =", ysup, " [Fe/H]max =", zsup ! Write(26, '(A, F3.1, A, F5.1, A)') "* M = ", mass_interp, + " Msun R = ", radius_interp, " Rsun" Write(26, '(A)') "*" Write(26, '(A)') "* LOG G" Write(26, '(F6.2)') logg_ref Write(26, '(A)') "*" Write(26, '(A)') "* NDEP" Write(26, '(I0)') ndepth_ref Write(26, '(A)') "*" If (MODELTYPE == "TAU5000") Then Write(26,'(A)') "* LOG(Tau5000) T[K] Ne [1/cm3] Vmacro[k +m/s] Vmicro[km/s]" Else If (MODELTYPE == "MASS") Then Write(26,'(A)') "* LOG(M) T[K] Ne [1/cm3] Vmacro[k +m/s] Vmicro[km/s]" End If ! Do k = 1, ndepth_ref If (Trim(MODELTYPE) == "TAU5000") Then Write(26,'(5ES12.4)')taus(k, out),T(k, out), + 10**(Pe(k, out))/(1.38065E-16 * T(k, out)), + 0.0, xit(k,out) Else If (Trim(MODELTYPE) == "MASS") Then Write(26, '(5ES12.4)') alog10(rhox(k, out)), T(k, out), + 10**(Pe(k, out))/(1.38065E-16 * T(k, out)), + 0.0, xit(k,out) End If End Do ! Close(26) ! ! dscale.marcs ! Open(26, file = "dscale." // FILE_IN(9), iostat = ios) IF(ios /= 0) STOP "WRITE_MULTI: PB-01." Write(26, '(A)') "MARCS" ! If (MODELTYPE == "TAU5000") Then Write(26, '(A)') "TAU5000 SCALE" Write(26, '(I4, F10.6)', iostat = ios)ndepth(out), rhox(1,out) IF (ios /= 0) STOP "WRITE_MULTI: PB-02." ! Do k = 1, ndepth_ref Write(26, '(F10.4)', iostat = ios) taus(k, out) IF (ios /= 0) STOP "WRITE_MULTI: PB-03." End Do ! Else If (MODELTYPE == "MASS") Then Write(26, '(A)') "MASS SCALE" Write(26, '(I4, F10.6)', iostat = ios)ndepth(out), taus(1,out) IF (ios /= 0) STOP "WRITE_MULTI: PB-04." ! Do k = 1, ndepth_ref Write(26, '(F10.4)', iostat = ios) alog10(rhox(k, out)) IF (ios /= 0) STOP "WRITE_MULTI: PB-05." End Do ! Else STOP"WRITE_MULTI: MODELTYPE UNKNOWN, ONLY 'TAU5000' OR 'MASS'" End If ! Close(26) ! ! abund.marcs ! Open(26, file = "abund." // FILE_IN(9), iostat = ios) IF(ios /= 0) STOP "WRITE_MULTI: PB-06" ! If (comp == 'st') Then If ( z_ref > -1.0 .AND. z_ref < 0.0) Then D_a = -0.4*z_ref Else If (z_ref <= -1.0) Then D_a = +0.4 Else If (z_ref >= 0.0) Then D_a = 0.0 End If Else If (comp == 'ap') Then D_a = 0. End If ! If (LOGASCALE) Then Write(26, '(A3, F5.2)') 'H ', solabu(1) Write(26, '(A3, F5.2)') 'HE ', solabu(2) Write(26, '(A3, F5.2)') 'SI ', solabu(3) + z_ref + D_a Write(26, '(A3, F5.2)') 'MG ', solabu(4) + z_ref + D_a Write(26, '(A3, F5.2)') 'AL ', solabu(5) + z_ref Write(26, '(A3, F5.2)') 'FE ', solabu(6) + z_ref Write(26, '(A3, F5.2)') 'C ', solabu(7) + z_ref Write(26, '(A3, F5.2)') 'NA ', solabu(8) + z_ref Write(26, '(A3, F5.2)') 'S ', solabu(9) + z_ref + D_a Write(26, '(A3, F5.2)') 'K ', solabu(10) + z_ref Write(26, '(A3, F5.2)') 'CA ', solabu(11) + z_ref + D_a Write(26, '(A3, F5.2)') 'NI ', solabu(12) + z_ref Write(26, '(A3, F5.2)') 'CR ', solabu(13) + z_ref Write(26, '(A3, F5.2)') 'N ', solabu(14) + z_ref Write(26, '(A3, F5.2)') 'O ', solabu(15) + z_ref Write(26, '(A3, F5.2)') 'NE ', solabu(16) + z_ref + D_a Else Write(26, '(A4,ES9.3)') 'H ', 10**(solabu(1) - 12) Write(26, '(A4,ES9.3)') 'HE ', 10**(solabu(2) - 12) Write(26, '(A4,ES9.3)') 'SI ', 10**(solabu(3)+z_ref+D_a-12) Write(26, '(A4,ES9.3)') 'MG ', 10**(solabu(4)+z_ref+D_a-12) Write(26, '(A4,ES9.3)') 'AL ', 10**(solabu(5) + z_ref - 12) Write(26, '(A4,ES9.3)') 'FE ', 10**(solabu(6) + z_ref - 12) Write(26, '(A4,ES9.3)') 'C ', 10**(solabu(7) + z_ref - 12) Write(26, '(A4,ES9.3)') 'NA ', 10**(solabu(8) + z_ref - 12) Write(26, '(A4,ES9.3)') 'S ', 10**(solabu(9)+z_ref+D_a-12) Write(26, '(A4,ES9.3)') 'K ', 10**(solabu(10) + z_ref - 12) Write(26, '(A4,ES9.3)') 'CA ', 10**(solabu(11)+z_ref+D_a-12) Write(26, '(A4,ES9.3)') 'NI ', 10**(solabu(12) + z_ref - 12) Write(26, '(A4,ES9.3)') 'CR ', 10**(solabu(13) + z_ref - 12) Write(26, '(A4,ES9.3)') 'N ', 10**(solabu(14) + z_ref - 12) Write(26, '(A4,ES9.3)') 'O ', 10**(solabu(15) + z_ref - 12) Write(26, '(A4,ES9.3)') 'NE ', 10**(solabu(16)+z_ref+D_a+12) End If ! Close(26) ! End Subroutine WRITE_MULTI ! !----------------------------------------------------------------- ! Subroutine WRITE_MOOG() ! Open(26, file = FILE_IN(9), iostat = ios) IF (ios /= 0) STOP "WRITE_MOOG: PB-00." Write(26, '(A, 2F7.1, 4F6.2)') +"WEB2MARC INTERPOLATED MODEL: ", xinf,xsup, yinf,ysup, zinf,zsup Write(26, '(A, F8.0, 2(A, F6.2))') + "Teff = ", temp_ref, " log g = ", logg_ref, " z = ", z_ref Write(26, '(A, I0)') "TAU ", ndepth_ref Write(26, '(F8.0)') xlr(out) ! Do k = 1, ndepth_ref Write(26, '(I2.2, ES11.3, F8.1, 3ES11.3)') k, taus(k, out), + T(k, out), Pe(k, out), Pg(k, out), rhox(k, out) End Do ! Write(26, *) xit(1, out) Write(26, '(A, F5.2)') "NATOMS 0 ", z_ref Write(26, '(A)') "NMOL 0" Write(26, '(A)') " 0.0" ! Close(26) ! End Subroutine WRITE_MOOG ! !----------------------------------------------------------------- ! Subroutine WRITE_TURBO() ! Open(unit = 23, file = FILE_IN(9), iostat = ios) IF (ios /= 0) STOP "WRITE_TURBO: PB-00." ! If (sph(1)) Then Write(*, *) 'Spherical models.' Write(23, 1967, iostat = ios) ndepth_ref, xlr(out), logg_ref IF (ios /= 0) STOP "WRITE_TURBO: PB-01." 1967 Format('''sphINTERPOL''', 1x, i3, f8.0, 2x, f4.2, 1x, '0 0.00') Do k = 1, ndepth_ref Write(23, 1968, iostat = ios) taus(k, out), T(k, out), + Pe(k, out), Pg(k, out), + xit(k, out), rr(k, out), + taur(k, out) IF (ios /= 0) STOP "WRITE_TURBO: PB-02." End Do 1968 Format(f8.4, 1x, f8.2, 3(1x, f8.4), 1x, e15.6, 1x, f8.4) Write(25, 19671, iostat = ios) ndepth_ref, xlr(out) IF (ios /= 0) STOP "WRITE_TURBO: PB-03." 19671 Format('sphINTERPOL', 1x, i3, f8.0) Else Write(*,*) 'Plane parallel models.' Write(23, 1966, iostat = ios) ndepth_ref, xlr(out), logg_ref IF (ios /= 0) STOP "WRITE_TURBO: PB-04." 1966 Format('''ppINTERPOL''',1x,i3,f8.0,2x,f4.2,1x,'0 0.00') Do k = 1, ndepth_ref Write(23, 1965, iostat = ios) taus(k, out), T(k, out), + Pe(k, out), Pg(k, out), + xit(k, out), rr(k, out), + taur(k, out) IF (ios /= 0) STOP "WRITE_TURBO: PB-05." 1965 Format(f8.4,1x,f8.2,3(1x,f8.4),1x,e15.6,1x,f8.4) End Do End If ! Write(23, *, iostat = ios) (TRIM(FILE_IN(f_in)), f_in = 1, 8) IF (ios /= 0) STOP "WRITE_TURBO: PB-06." ! Close(23) ! End Subroutine WRITE_TURBO ! !----------------------------------------------------------------- ! Subroutine WRITE_ATLAS() ! Open(unit = 25, file = FILE_IN(9), iostat = ios) IF (ios /= 0) STOP "WRITE_ATLAS: PB-00." ! If (sph(1)) Then Write(*, *) 'Spherical models.' Write(25, 19671, iostat = ios) ndepth_ref, xlr(out) IF (ios /= 0) STOP "WRITE_TURBO: PB-01." 19671 Format('sphINTERPOL', 1x, i3, f8.0) Write(25, 19672, iostat = ios) IF (ios /= 0) STOP "WRITE_TURBO: PB-02." 19672 Format(' k log(tau) T log(Pe) log(Pg) rhox') Do k = 1, ndepth_ref Write(25,19681, iostat = ios) k, taus(k, out), T(k, out), + Pe(k, out), Pg(k, out), + rhox(k,out) IF (ios /= 0) STOP "WRITE_TURBO: PB-03." 19681 Format(i5, 1x, f8.4, 1x, f8.2, 2(1x,f8.4), 1x, e15.6) End Do Else Write(*,*) 'Plane parallel models.' Write(25, 19661, iostat = ios) xlr(out) IF (ios /= 0) STOP "WRITE_TURBO: PB-04." 19661 Format('ppINTERPOL',f8.0) Write(25,19672, iostat = ios) IF (ios /= 0) STOP "WRITE_TURBO: PB-05." Do k = 1, ndepth_ref Write(25, 19681, iostat = ios) k, taus(k, out), T(k, out), + Pe(k, out), Pg(k, out), + rhox(k, out) IF (ios /= 0) STOP "WRITE_TURBO: PB-06." End Do End If ! Write(25, *, iostat = ios) (TRIM(FILE_IN(f_in)), f_in = 1, 8) IF (ios /= 0) STOP "WRITE_TURBO: PB-07." ! Close(25) ! End Subroutine WRITE_ATLAS ! !----------------------------------------------------------------- ! End Program interpol_modeles ! ************************************************************************ c----------------------------------------------------------------------- subroutine blend_103 (r,s,t,x000,x001,x010,x011,x100,x101,x110, & x111, x ) ! !*********************************************************************** !from http://www.math.iastate.edu/burkardt/f_src/ ! !! BLEND_103 extends scalar point data into a cube. ! ! ! Diagram: ! ! 011--------------111 ! | | ! | | ! | | ! | | ! | | ! 001--------------101 ! ! ! *---------------* ! | | ! | | ! | rst | ! | | ! | | ! *---------------* ! ! ! 010--------------110 ! | | ! | | ! | | ! | | ! | | ! 000--------------100 ! ! ! Formula: ! ! Written as a polynomial in R, S and T, the interpolation map has the ! form: ! ! X(R,S,T) = ! 1 * ( + x000 ) ! + r * ( - x000 + x100 ) ! + s * ( - x000 + x010 ) ! + t * ( - x000 + x001 ) ! + r * s * ( + x000 - x100 - x010 + x110 ) ! + r * t * ( + x000 - x100 - x001 + x101 ) ! + s * t * ( + x000 - x010 - x001 + x011 ) ! + r * s * t * ( - x000 + x100 + x010 + x001 - x011 - x101 - x110 + x111 ) ! ! Reference: ! ! William Gordon, ! Blending-Function Methods of Bivariate and Multivariate Interpolation ! and Approximation, ! SIAM Journal on Numerical Analysis, ! Volume 8, Number 1, March 1971, pages 158-177. ! ! William Gordon and Charles Hall, ! Transfinite Element Methods: Blending-Function Interpolation over ! Arbitrary Curved Element Domains, ! Numerische Mathematik, ! Volume 21, Number 1, 1973, pages 109-129. ! ! William Gordon and Charles Hall, ! Construction of Curvilinear Coordinate Systems and Application to ! Mesh Generation, ! International Journal of Numerical Methods in Engineering, ! Volume 7, 1973, pages 461-477. ! ! Joe Thompson, Bharat Soni, Nigel Weatherill, ! Handbook of Grid Generation, ! CRC Press, 1999. ! ! Modified: ! ! 11 October 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, T, the coordinates where an interpolated value ! is desired. ! ! Input, real X000, X001, X010, X011, X100, X101, X110, X111, the ! data values at the corners. ! ! Output, real X, the interpolated data value at (R,S,T). ! implicit none ! real r real s real t real x real x000 real x001 real x010 real x011 real x100 real x101 real x110 real x111 ! ! Interpolate the interior point. ! x = & 1.0E+00 * ( + x000 ) & + r * ( - x000 + x100 ) & + s * ( - x000 + x010 ) & + t * ( - x000 + x001 ) & + r * s * ( + x000 - x100 - x010 + x110) & + r * t * ( + x000 - x100 - x001 + x101 ) & + s * t * ( + x000 - x010 - x001 + x011 ) & + r * s * t * ( - x000 + x100 + x010 + x001 - x011 - x101 - x110+ & x111 ) return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine extract_bin(FILE,TEFF,grav,metal,ndepth,xlr_ref,tau5, & tauR,temp,prese,presg,xit,rad,sph,kappa) !extracted from osplot.f 07/2003, to get tau,T,Pe,Pg,microturb from a model implicit none Integer :: ios integer :: ndp,ndepth,k,nlp,nlb parameter(ndp=200) CHARACTER*117 ADUM CHARACTER*256 FILE,file2 CHARACTER*30 COMMENT real metal,grav,xlr_ref,TEFF,GG,radius logical :: sph real :: tau(ndp),t(ndp),z(ndp),ptot(ndp),prad(ndp), & pg(ndp),pturb(ndp),pe(ndp),ro(ndp),xmass(ndp),xkapr(ndp) real :: gradp(ndp),gravity(ndp),pcheck(ndp),rr(ndp), & xit(ndp),geff(ndp),gradptur(ndp),dp(ndp),taus(ndp),xlr(30), & coldens(ndp) real :: tau5(ndp),tauR(ndp),temp(ndp),prese(ndp), & presg(ndp),rad(ndp),kappa(ndp) real :: xlb(155000),w(155000),fluxme(155000) real :: presmo(30,ndp),ptio(ndp) real :: bPPR(NDP),bPPT(NDP),bPP(NDP),bGG(NDP), & bZZ(NDP),bDD(NDP), & bVV(NDP),bFFC(NDP),bPPE(NDP),bTT(NDP), & bTAULN(NDP),erad(ndp) integer :: NbTAU,IbTER common /struct/ tau,t,z,ptot,prad,pg,pturb,pe,ro,rr,taus,xlr, & nlp,xkapr common /spectr/ nlb,xlb,w,fluxme common /pressure/ presmo,ptio common /binstruc/bPPR,bPPT,bPP,bGG, & bZZ,bDD, & bVV,bFFC,bPPE,bTT, & bTAULN,NbTAU,IbTER,erad common /radius/ radius OPEN(UNIT=10,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED', & iostat = ios) IF (ios /= 0) STOP "Main : PB-23" c & convert='big_endian') ccc & RECL=152600) * CALL READMO(10,NDEPTH,TEFF,GG,metal,sph) c open(21,file=file2,status='unknown', iostat = ios) IF (ios /= 0) STOP "Main : PB-24" do k=1,ndepth tau5(k)=log10(taus(k)) tauR(k)=log10(tau(k)) temp(k)=t(k) prese(k)=log10(pe(k)) presg(k)= log10(pg(k)) kappa(k)=log10(xkapr(k)) end do xit=2.0 xlr_ref=xlr(nlp) grav=log10(GG) * rad=rr if(.not.sph) radius=0.0 rad=radius-z c write(21,1966, iostat = ios) ndepth,xlr(nlp),log10(GG) IF (ios /= 0) STOP "Main : PB-25" c1966 format('''INTERPOL''',1x,i3,f8.0,2x,f4.2,1x,'0 0.00') c do k=1,ndepth c write(21,1965, isotat=ios) log10(taus(k)),t(k),log10(pe(k)) c & log10(pg(k)), xit c1965 format(f8.4,1x,f8.2,3(x,f8.4)) IF (ios /= 0) STOP "Main : PB-26" c enddo c close(21) close(10) END C SUBROUTINE READMO(IARCH,JTAU,TEFF,G,metal,spherical) C THIS ROUTINE READS ONE MODEL, TO GET INFO ON PRAD C ( All features taken from listmo ) PARAMETER (NDP=200) C DIMENSION ABUND(16),TKORRM(NDP),FCORR(NDP),TAU(NDP),TAUS(NDP), *T(NDP),PE(NDP),PG(NDP),PRAD(NDP),PTURB(NDP),XKAPR(NDP),RO(NDP), *CP(NDP),CV(NDP),AGRAD(NDP),Q(NDP),U(NDP),V(NDP),ANCONV(NDP), *PRESMO(30,NDP),FCONV(NDP),RR(NDP),Z(NDP),EMU(NDP),HNIC(NDP) *,NJ(16),XLR(30),IEL(16),ANJON(16,5),PART(16,5),PROV(50,20+1), *ABSKA(50),SPRIDA(50),XLB(155000),FLUXME(155000),FLUMAG(155000), & PEP(16), * ABNAME(50),SOURCE(50),PTOT(NDP) DIMENSION W(155000),UW(12),BW(21),VW(25) CHARACTER*10 DAG,NAME,NAMEP,KLOCK CHARACTER*8 ABNAME,SOURCE DIMENSION WAVFLX(10) dimension PTIO(NDP) real*8 dluminosity real abSc,abTi,abV,abMn,abCo,metal logical spherical common /binstruc/ dummy(11*ndp+2),erad(ndp) common /struct/ tau,t,z,ptot,prad,pg,pturb,pe,ro,rr,taus,xlr, & nlp,xkapr common /spectr/ nlb,xlb,w,fluxme common /pressure/ presmo,ptio common /radius/ radius DATA UW/0.145,0.436,0.910,1.385,1.843,2.126,2.305,2.241,1.270, *0.360,0.128,0.028/,BW/0.003,0.026,0.179,0.612,1.903,2.615,2.912, *3.005,2.990,2.876,2.681,2.388,2.058,1.725,1.416,1.135,0.840,0.568, *0.318,0.126,0.019/,VW/0.006,0.077,0.434,1.455,2.207,2.703,2.872, *2.738,2.505,2.219,1.890,1.567,1.233,0.918,0.680,0.474,0.312,0.200, *0.132,0.096,0.069,0.053,0.037,0.022,0.012/ DATA NAME/'LOCAL'/,NAMEP/'PARSONS'/ DATA A,B/.34785485,.65214515/ IREAD=5 C REWIND IARCH READ(IARCH) erad=-1.e30 READ(IARCH) INORD,DAG,KLOCK READ(IARCH) TEFF,FLUX,G,PALFA,PNY,PY,PBETA,ILINE,ISTRAL, & MIHAL,IDRAB1,IDRAB2,IDRAB3,IDRAB4,IDRAB5,IDRAB6, & ITMAX,NEL,(ABUND(I),I=1,NEL),abSc,abTi,abV,abMn,abCo GLOG=ALOG10(G) FNORD=0.1*INORD C CONVERT TO 'PHYSICAL FLUX' FLUX=3.14159*FLUX DO 2 I=1,NEL 2 ABUND(I)=ALOG10(ABUND(I))+12. metal=abund(15)-7.50 READ(IARCH)JTAU,NCORE,DIFLOG,TAUM,RADIUS,(RR(K),K=1,JTAU) if (jtau.gt.ndp) then print*, 'ERROR !!! Jtau (number of depths of model) = ',jtau print*, ' is larger than ndp!! Increase NDP.' stop endif if (radius.le.2.) then spherical = .false. else spherical = .true. rr=radius-rr endif READ(IARCH)JTAU,(TKORRM(I),I=1,JTAU),(FCORR(K),K=1,JTAU) NTPO=0 DO 3 K=1,JTAU READ(IARCH) KR,TAU(K),TAUS(K),Z(K),T(K),PE(K),PG(K),PRAD(K), & PTURB(K),XKAPR(K),RO(K),EMU(K),CP(K),CV(K), & AGRAD(K),Q(K),U(K),V(K),ANCONV(K),HNIC(K),NMOL, & (PRESMO(J,K),J=1,NMOL),ptio(k) TAUK=ALOG10(TAU(K))+10.01 KTAU=TAUK IF(ABS(TAUK-KTAU).GT.0.02) GO TO 31 IF(KTAU.EQ.10) K0=K NTPO=NTPO+1 31 CONTINUE 3 CONTINUE c Z0=Z(K0) c DO 5 I=1,JTAU c Z(I)=Z(I)-Z0 c i1=min(i+1,jtau) c PTOT(I)=PG(I)+PRAD(I)+0.5*(pturb(i)+pturb(i1)) 5 CONTINUE *** READ(IARCH)(NJ(I),I=1,NEL),NLP,(XLR(I),I=1,NLP) & ,NPROV,NPROVA,NPROVS,(ABNAME(KP),SOURCE(KP),KP=1,NPROV) c DO 22 KTAU=1,NTPO c DO 20 IE=1,NEL c NJP=NJ(IE) c READ(IARCH) KR,TAUI,TI,PEI,IEL(IE),ABUND(IE), c & (ANJON(IE,JJ),JJ=1,NJP),(PART(IE,JJ),JJ=1,NJP) c 20 CONTINUE c DO 21 KLAM=1,NLP c READ(IARCH) KR,TAUIL,(PROV(J,KLAM),J=1,NPROV), c & ABSKA(KLAM),SPRIDA(KLAM) 21 CONTINUE 22 continue c READ(IARCH) NLB,(XLB(J),FLUXME(J),J=1,NLB),(W(J),J=1,NLB) C CONVERT TO 'PHYSICAL' FLUXES c DO 24 J=1,NLB c 24 FLUXME(J)=3.14159*FLUXME(J) c dluminosity=0. c do 25 j=1,nlb c dluminosity=dluminosity+fluxme(j)*w(j) 25 continue c dluminosity=dluminosity*4.*3.14159*radius**2/3.82d33 c ddddd=real(dluminosity) c print*,'luminosity: ',dluminosity*3.82d33,' erg/s = ',ddddd, c & ' solar luminosities' *** RETURN END ************************************************************************ subroutine resample(taus,tauR,T,Pe,Pg,xit,rr,xkapref) implicit none integer :: nlinemod,file,nfile,k,i real,dimension(:,:) :: taus,tauR,T,Pe,Pg,xit,rr,xkapref real,dimension(:,:),allocatable :: taubas real,dimension(:),allocatable :: tauresample,taustemp,tauRtemp, & Ttemp,Petemp,Pgtemp,xittemp,rrtemp,xkapreftemp INTERFACE function SevalSingle(u,x,y) REAL,INTENT(IN) :: u ! abscissa at which the spline is to be evaluated REAL,INTENT(IN),DIMENSION(:) :: x ! abscissas of knots REAL,INTENT(IN),DIMENSION(:):: y ! ordinates of knots real SevalSingle end END INTERFACE nlinemod = size(taus, 1) nfile=size(taus,2)-3 allocate(tauresample(nlinemod),taustemp(nlinemod), & tauRtemp(nlinemod),Ttemp(nlinemod) & ,Petemp(nlinemod),Pgtemp(nlinemod),xittemp(nlinemod), & rrtemp(nlinemod),xkapreftemp(nlinemod), & taubas(nlinemod,nfile)) !!!!choose here the depth basis for interpolation (tau5000,tauRoss)!!!!! taubas = tauR(:, : nfile) Write(*, *) 'Resample models on common depth basis: tauRoss.' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call common_depth_basis(tauresample,taubas,nlinemod,nfile) c now do the resampling with the common tau do file=1,nfile do k=1,nlinemod taustemp(k)=SevalSingle(tauresample(k),taubas(:,file) & ,taus(:,file)) tauRtemp(k)=SevalSingle(tauresample(k),taubas(:,file) & ,tauR(:,file)) Ttemp(k)=SevalSingle(tauresample(k),taubas(:,file),T(:,file)) Petemp(k)=SevalSingle(tauresample(k),taubas(:,file),Pe(:,file)) Pgtemp(k)=SevalSingle(tauresample(k),taubas(:,file),Pg(:,file)) xittemp(k)=SevalSingle(tauresample(k),taubas(:,file),xit(:,file)) rrtemp(k)=SevalSingle(tauresample(k),taubas(:,file),rr(:,file)) xkapreftemp(k)= & SevalSingle(tauresample(k),taubas(:,file),xkapref(:,file)) end do taus(:,file)=taustemp tauR(:,file)=tauRtemp T(:,file)=Ttemp Pe(:,file)=Petemp Pg(:,file)=Pgtemp xit(:,file)=xittemp rr(:,file)=rrtemp xkapref(:,file)=xkapreftemp end do deallocate(tauresample,taustemp,tauRtemp,Ttemp,Petemp &,Pgtemp,xittemp,rrtemp,xkapreftemp,taubas) end !*********************************************************************** subroutine common_depth_basis(tauresample,tau,nlinemod,nfile) implicit none integer :: file,nlinemod,nfile real,dimension(nlinemod,nfile) :: tau real,dimension(nlinemod) :: tauresample tauresample=0 c initialize the common tau(5000) with min depth = max of the min depth of the models c and max depth = min of the max depth of the models c essential for the resampling with cubic spline tauresample(1)=tau(1,1) tauresample(nlinemod)=tau(nlinemod,1) do file=2,nfile if (tauresample(1).lt.tau(1,file)) then tauresample(1)=tau(1,file) end if if (tauresample(nlinemod).gt.tau(nlinemod,file)) then tauresample(nlinemod)=tau(nlinemod,file) end if end do call blend_i_0d1 ( tauresample, nlinemod ) end !*********************************************************************** subroutine blend_i_0d1 ( x, m ) ! ! !! BLEND_I_0D1 extends indexed scalar data at endpoints along a line. !http://orion.math.iastate.edu/burkardt/f_src/f_src.html ! ! Diagram: ! ! ( X1, ..., ..., ..., ..., ..., XM ) ! ! Reference: ! ! William Gordon, ! Blending-Function Methods of Bivariate and Multivariate Interpolation ! and Approximation, ! SIAM Journal on Numerical Analysis, ! Volume 8, Number 1, March 1971, pages 158-177. ! ! William Gordon and Charles Hall, ! Transfinite Element Methods: Blending-Function Interpolation over ! Arbitrary Curved Element Domains, ! Numerische Mathematik, ! Volume 21, Number 1, 1973, pages 109-129. ! ! William Gordon and Charles Hall, ! Construction of Curvilinear Coordinate Systems and Application to ! Mesh Generation, ! International Journal of Numerical Methods in Engineering, ! Volume 7, 1973, pages 461-477. ! ! Joe Thompson, Bharat Soni, Nigel Weatherill, ! Handbook of Grid Generation, ! CRC Press, 1999. ! ! Modified: ! ! 15 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X(M). ! ! On input, X(1) and X(M) contain scalar values which are to be ! interpolated through the entries X(2) through X(M). It is assumed ! that the dependence of the data is linear in the vector index I. ! ! On output, X(2) through X(M-1) have been assigned interpolated ! values. ! ! Input, integer M, the number of entries in X. implicit none ! integer m ! integer i real r real x(m) ! do i = 2, m - 1 r = real ( i - 1 ) / real ( m - 1 ) call blend_101 ( r, x(1), x(m), x(i) ) end do return end subroutine blend_101 ( r, x0, x1, x ) ! !*********************************************************************** ! !! BLEND_101 extends scalar endpoint data to a line. !http://orion.math.iastate.edu/burkardt/f_src/f_src.html ! ! Diagram: ! ! 0-----r-----1 ! ! Reference: ! ! William Gordon, ! Blending-Function Methods of Bivariate and Multivariate Interpolation ! and Approximation, ! SIAM Journal on Numerical Analysis, ! Volume 8, Number 1, March 1971, pages 158-177. ! ! William Gordon and Charles Hall, ! Transfinite Element Methods: Blending-Function Interpolation over ! Arbitrary Curved Element Domains, ! Numerische Mathematik, ! Volume 21, Number 1, 1973, pages 109-129. ! ! William Gordon and Charles Hall, ! Construction of Curvilinear Coordinate Systems and Application to ! Mesh Generation, ! International Journal of Numerical Methods in Engineering, ! Volume 7, 1973, pages 461-477. ! ! Joe Thompson, Bharat Soni, Nigel Weatherill, ! Handbook of Grid Generation, ! CRC Press, 1999. ! ! Modified: ! ! 14 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, the coordinate where an interpolated value is desired. ! ! Input, real X0, X1, the data values at the ends of the line. ! ! Output, real X, the interpolated data value at (R). ! implicit none ! real r real x real x0 real x1 ! x = ( 1.0E+00 - r ) * x0 + r * x1 return end !---------------------------------------------------------------------------- REAL FUNCTION SevalSingle(u,x,y) ! --------------------------------------------------------------------------- !http://www.pdas.com/fmm.htm ! PURPOSE - Evaluate the cubic spline function ! Seval=y(i)+b(i)!(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3 ! where x(i) <= u < x(i+1) ! NOTES- if ux(n), i=n is used REAL,INTENT(IN) :: u ! abscissa at which the spline is to be evaluated REAL,INTENT(IN),DIMENSION(:) :: x ! abscissas of knots REAL,INTENT(IN),DIMENSION(:):: y ! ordinates of knots REAL, DIMENSION(:),allocatable :: b,c,d ! linear,quadratic,cubic coeff INTEGER, SAVE :: i=1 INTEGER :: j, k, n REAL:: dx INTERFACE subroutine FMMsplineSingle(x, y, b, c, d) REAL,DIMENSION(:), INTENT(IN) :: x ! abscissas of knots REAL,DIMENSION(:), INTENT(IN) :: y ! ordinates of knots REAL,DIMENSION(:), INTENT(OUT) :: b ! linear coeff REAL,DIMENSION(:), INTENT(OUT) :: c ! quadratic coeff. REAL,DIMENSION(:), INTENT(OUT) :: d ! cubic coeff. end END INTERFACE !------------------------------------------------------------------------ n=SIZE(x) allocate(b(n),c(n),d(n)) call FMMsplineSingle(x, y, b, c, d) !.....First check if u is in the same interval found on the ! last call to Seval............................................. IF ( (i<1) .OR. (i >= n) ) i=1 IF ( (u < x(i)) .OR. (u >= x(i+1)) ) THEN i=1 ! binary search j=n+1 DO k=(i+j)/2 IF (u < x(k)) THEN j=k ELSE i=k END IF IF (j <= i+1) EXIT END DO END IF dx=u-x(i) ! evaluate the spline SevalSingle=y(i)+dx*(b(i)+dx*(c(i)+dx*d(i))) RETURN deallocate(b,c,d) END Function SevalSingle ! -------------------------------------- SUBROUTINE FMMsplineSingle(x, y, b, c, d) !----------------------------------------------------------------------- !http://www.pdas.com/fmm.htm ! PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline ! so that the interpolated value is given by ! s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3 ! when x(k) <= x <= x(k+1) ! The end conditions match the third derivatives of the interpolated curve to ! the third derivatives of the unique polynomials thru the first four and ! last four points. ! Use Seval or Seval3 to evaluate the spline. REAL,DIMENSION(:), INTENT(IN) :: x ! abscissas of knots REAL,DIMENSION(:), INTENT(IN) :: y ! ordinates of knots REAL,DIMENSION(:), INTENT(OUT) :: b ! linear coeff REAL,DIMENSION(:), INTENT(OUT) :: c ! quadratic coeff. REAL,DIMENSION(:), INTENT(OUT) :: d ! cubic coeff. INTEGER:: k,n REAL:: t,aux REAL,PARAMETER:: ZERO=0.0, TWO=2.0, THREE=3.0 !----------------------------------------------------------------------- n=SIZE(x) IF (n < 3) THEN ! Straight line - special case for n < 3 b(1)=ZERO IF (n == 2) b(1)=(y(2)-y(1))/(x(2)-x(1)) c(1)=ZERO d(1)=ZERO IF (n < 2) RETURN b(2)=b(1) c(2)=ZERO d(2)=ZERO RETURN END IF !.....Set up tridiagonal system......................................... !. b=diagonal, d=offdiagonal, c=right-hand side d(1)=x(2)-x(1) c(2)=(y(2)-y(1))/d(1) DO k=2,n-1 d(k)=x(k+1)-x(k) b(k)=TWO*(d(k-1)+d(k)) c(k+1)=(y(k+1)-y(k))/d(k) c(k)=c(k+1)-c(k) END DO !.....End conditions. third derivatives at x(1) and x(n) obtained !. from divided differences....................................... b(1)=-d(1) b(n)=-d(n-1) c(1)=ZERO c(n)=ZERO IF (n > 3) THEN c(1)=c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1)) c(n)=c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3)) c(1)=c(1)*d(1)*d(1)/(x(4)-x(1)) c(n)=-c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3)) END IF DO k=2,n ! forward elimination t=d(k-1)/b(k-1) b(k)=b(k)-t*d(k-1) c(k)=c(k)-t*c(k-1) END DO c(n)=c(n)/b(n) ! back substitution ( makes c the sigma of text) DO k=n-1,1,-1 c(k)=(c(k)-d(k)*c(k+1))/b(k) END DO !.....Compute polynomial coefficients................................... b(n)=(y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n)) DO k=1,n-1 b(k)=(y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k)) d(k)=(c(k+1)-c(k))/d(k) c(k)=THREE*c(k) END DO c(n)=THREE*c(n) d(n)=d(n-1) RETURN END Subroutine FMMsplineSingle ! --------------------------------------------------- SUBROUTINE NaturalSplineSingle(x,y,b,c,d) ! --------------------------------------------------------------------------- ! PURPOSE - Construct the natural spline thru a set of points ! NOTES - A natural spline has zero second derivative at both endpoints. REAL,INTENT(IN),DIMENSION(:):: x,y ! coordinates of knots REAL,INTENT(OUT),DIMENSION(:):: b,c,d ! cubic coeff. INTEGER:: k,n REAL,PARAMETER:: ZERO=0.0, TWO=2.0, THREE=3.0 !----------------------------------------------------------------------- n=SIZE(x) IF (n < 3) THEN ! Straight line - special case for n < 3 b(1)=ZERO IF (n == 2) b(1)=(y(2)-y(1))/(x(2)-x(1)) c(1)=ZERO d(1)=ZERO b(2)=b(1) c(2)=ZERO d(2)=ZERO RETURN END IF d(1:n-1) = x(2:n)-x(1:n-1) ! Put the h-array of the text into array d !.....Set up the upper triangular system in locations 2 thru n-1 of ! arrays b and c. B holds the diagonal and c the right hand side. b(2)=TWO*(d(1)+d(2)) c(2)=(y(3)-y(2))/d(2)-(y(2)-y(1))/d(1) DO k=3,n-1 b(k)=TWO*(d(k-1)+d(k))-d(k-1)*d(k-1)/b(k-1) c(k)=(y(k+1)-y(k))/d(k)-(y(k)-y(k-1))/d(k-1)-d(k-1)*c(k-1)/b(k-1) END DO c(n-1)=c(n-1)/b(n-1) ! Back substitute to get c-array DO k=n-2,2,-1 c(k)=(c(k)-d(k)*c(k+1))/b(k) END DO c(1)=ZERO c(n)=ZERO ! c now holds the sigma array of the text !.....Compute polynomial coefficients .................................. b(n)=(y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n)) DO k=1,n-1 b(k)=(y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k)) d(k)=(c(k+1)-c(k))/d(k) c(k)=THREE*c(k) END DO c(n)=THREE*c(n) d(n)=d(n-1) RETURN END Subroutine NaturalSplineSingle real function rinteg(x,f,fint,n,start) ************************************************************************ * This routine is from ATLAS6 ************************************************************************ implicit none integer :: n,i,n1 real x(5000), f(5000), fint(5000) real a(5000), b(5000), c(5000) real :: start call parcoe (f,x,a,b,c,n) fint(1) = start rinteg = start n1 = n - 1 do 10 i=1,n1 fint(i+1)= (a(i)+b(i)/2.*(x(i+1)+x(i))+ . c(i)/3.*((x(i+1)+x(i))*x(i+1)+x(i)*x(i)))*(x(i+1)-x(i)) 10 rinteg = rinteg + fint(i+1) return end ************************************************************************ subroutine parcoe(f,x,a,b,c,n) implicit none integer :: n,n1,j,j1 real f(5000), x(5000), a(5000), b(5000), c(5000) real :: d,wt c(1)=0. b(1)=(f(2)-f(1))/(x(2)-x(1)) a(1)=f(1)-x(1)*b(1) n1=n-1 c(n)=0. b(n)=(f(n)-f(n1))/(x(n)-x(n1)) a(n)=f(n)-x(n)*b(n) if(n.eq.2)return do 1 j=2,n1 j1=j-1 d=(f(j)-f(j1))/(x(j)-x(j1)) c(j)=f(j+1)/((x(j+1)-x(j))*(x(j+1)-x(j1)))-f(j)/((x(j)-x(j1))* 1(x(j+1)-x(j)))+f(j1)/((x(j)-x(j1))*(x(j+1)-x(j1))) b(j)=d-(x(j)+x(j1))*c(j) 1 a(j)=f(j1)-x(j1)*d+x(j)*x(j1)*c(j) c(2)=0. b(2)=(f(3)-f(2))/(x(3)-x(2)) a(2)=f(2)-x(2)*b(2) c(3)=0. b(3)=(f(4)-f(3))/(x(4)-x(3)) a(3)=f(3)-x(3)*b(3) if(c(j).eq.0.)go to 2 j1=j+1 wt=abs(c(j1))/(abs(c(j1))+abs(c(j))) a(j)=a(j1)+wt*(a(j)-a(j1)) b(j)=b(j1)+wt*(b(j)-b(j1)) c(j)=c(j1)+wt*(c(j)-c(j1)) 2 continue a(n1)=a(n) b(n1)=b(n) c(n1)=c(n) return end