From a0892a289d5e1488a1d8fa1cebaef06220aee4bc Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 11 Dec 2024 11:53:59 +1100 Subject: [PATCH] Removing unused variables --- src/ML_prodimo.f90 | 5 ++--- src/SPH2mcfost.f90 | 4 ++-- src/benchmarks.f90 | 2 +- src/diffusion.f90 | 5 +---- src/disk_physics.f90 | 2 +- src/dust_prop.f90 | 4 ++-- src/gas/abo.f90 | 19 +++++++++--------- src/gas/atom_type.f90 | 24 +++++++++++------------ src/gas/broad.f90 | 16 +++++++-------- src/gas/elements_type.f90 | 4 +--- src/gas/gas_contopac.f90 | 39 +++++++++++++++++-------------------- src/gas/lte.f90 | 27 +++++++++++-------------- src/gas/opacity_atom.f90 | 35 +++++++++++++++++---------------- src/gas/wavelengths_gas.f90 | 19 +++++++++--------- src/input.f90 | 2 +- src/mcfost_env.f90 | 2 +- src/molecular_emission.f90 | 2 +- src/readVTK.f90 | 5 +---- src/read_fargo3d.f90 | 3 ++- src/read_idefix.f90 | 13 +++---------- src/read_param.f90 | 2 +- src/read_phantom.f90 | 2 +- src/read_pluto.f90 | 1 + src/read_spherical_grid.f90 | 8 ++++---- src/scattering.f90 | 38 +++++++++++++++++++----------------- src/stars.f90 | 8 +++----- src/voigts.f90 | 23 +++++++++++----------- 27 files changed, 144 insertions(+), 170 deletions(-) diff --git a/src/ML_prodimo.f90 b/src/ML_prodimo.f90 index be2624221..de50db641 100644 --- a/src/ML_prodimo.f90 +++ b/src/ML_prodimo.f90 @@ -95,9 +95,8 @@ subroutine save_J_ML(lambda, lISM) ! avant et apres le calcul du champ ISM use dust_ray_tracing, only : n_phot_envoyes - use radiation_field, only : xJ_abs, xN_abs + use radiation_field, only : xJ_abs use prodimo, only : chi_ISM, R_ISM, T_ISM_stars, Wdil, Tcmb - use parametres, only : Rmax use temperature, only : E_disk use cylindrical_grid, only : volume use wavelengths, only : tab_lambda @@ -151,7 +150,7 @@ subroutine xgb_compute_features() use optical_depth, only : compute_column use density, only : densite_pouss use cylindrical_grid, only : r_grid, z_grid - use molecular_emission, only : densite_gaz, Tcin + use molecular_emission, only : densite_gaz use temperature, only : Tdust use Voronoi_grid, only : Voronoi diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 54df20aa3..b69feef6c 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -222,9 +222,9 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g real, allocatable, dimension(:) :: a_SPH, log_a_SPH, rho_dust real(dp), allocatable, dimension(:) :: gsize, grainsize_f, dN_ds, N_monomers, rho_monomers - real(dp), dimension(4) :: lambsol, lambguess + real(dp), dimension(4) :: lambsol - real(dp) :: mass, somme, Mtot, Mtot_dust, facteur, a, mass_factor + real(dp) :: mass, Mtot, Mtot_dust, facteur, a, mass_factor real :: f, limit_threshold, density_factor integer :: icell, l, k, iSPH, n_force_empty, i, id_n, ierr, N_pb diff --git a/src/benchmarks.f90 b/src/benchmarks.f90 index 9efef7162..89affda4c 100644 --- a/src/benchmarks.f90 +++ b/src/benchmarks.f90 @@ -14,7 +14,7 @@ module benchmarks subroutine init_Pascucci_benchmark() - character(len=32) :: filename = "Pascucci_optSi.dat" + !character(len=32) :: filename = "Pascucci_optSi.dat" write(*,*) "------------------------------------------------" write(*,*) "| Setting up the Pascucci et al 2004 benchmark |" diff --git a/src/diffusion.f90 b/src/diffusion.f90 index 3021dde9d..ac48996cb 100644 --- a/src/diffusion.f90 +++ b/src/diffusion.f90 @@ -640,15 +640,12 @@ subroutine compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity) use wavelengths, only : n_lambda, tab_lambda, tab_delta_lambda use Temperature, only : Tdust use dust_prop, only : kappa - use Voronoi_grid, only : Voronoi - use cylindrical_grid, only : volume - use density, only : masse_gaz, densite_gaz integer, intent(in) :: icell real(dp), intent(out) :: Planck_opacity,rec_Planck_opacity ! cm2/g (ie per gram of gas) integer :: lambda - real(dp) :: somme, somme2, cst, cst_wl, B, dB_dT, coeff_exp, wl, delta_wl, norm, T + real(dp) :: somme, somme2, cst, cst_wl, B, coeff_exp, wl, delta_wl, norm, T !dB_dT integer, pointer :: p_icell integer, target :: icell0 diff --git a/src/disk_physics.f90 b/src/disk_physics.f90 index 2ec9db25c..c8bf862e0 100644 --- a/src/disk_physics.f90 +++ b/src/disk_physics.f90 @@ -185,7 +185,7 @@ subroutine equilibre_hydrostatique() implicit none real, dimension(nz) :: rho, ln_rho - real :: dz, dz_m1, dTdz, fac, fac1, fac2, M_etoiles, M_mol, somme, cst + real :: dz, dz_m1, dTdz, fac1, fac2, M_etoiles, M_mol, somme, cst integer :: i,j, k, icell, icell_m1 real, parameter :: gas_dust = 100 diff --git a/src/dust_prop.f90 b/src/dust_prop.f90 index 3cb740570..9ca16e590 100644 --- a/src/dust_prop.f90 +++ b/src/dust_prop.f90 @@ -553,7 +553,7 @@ subroutine prop_grains(lambda) implicit none integer, intent(in) :: lambda - real :: a, wavel, x, qext, qsca, gsca, amu1, amu2, amu1_coat, amu2_coat, Cext, Csca + real :: a, wavel, x, qext, qsca, gsca, amu1, amu2, amu1_coat, amu2_coat integer :: k, pop qext=0.0 @@ -795,7 +795,7 @@ subroutine opacite(lambda, p_lambda, no_scatt) integer, intent(in) :: lambda, p_lambda logical, intent(in), optional :: no_scatt - integer :: icell, k, thetaj + integer :: icell, k real(kind=dp) :: density, fact, k_abs_RE, k_abs_LTE, k_abs_tot, k_sca_tot, rho0 logical :: lcompute_obs, ldens0, compute_scatt diff --git a/src/gas/abo.f90 b/src/gas/abo.f90 index 759fab0cd..f2d6713da 100644 --- a/src/gas/abo.f90 +++ b/src/gas/abo.f90 @@ -73,7 +73,7 @@ end function gammln subroutine dealloc_barklem(barklem) type(barklem_cross_data), intent(inout) :: barklem deallocate(barklem%neff1, barklem%neff2, barklem%cross,barklem%alpha) - return + return end subroutine dealloc_barklem @@ -168,9 +168,8 @@ subroutine get_Barklem_cross_data(atom, kr, res) real(kind=dp) :: neff1, neff2, nefftmp real(kind=dp) :: findex1, findex2, mu, vmean, sqa0 real(kind=dp) :: cross, E_Rydberg2, deltaEi, deltaEj, expo - integer :: Li, Lj, i, j, Z, ic, k, index, index1 + integer :: Li, Lj, i, j, Z, ic, index, index1 character(len=ATOM_LABEL_WIDTH) :: label_i, label_j - real(kind=dp) :: testcc i = Atom%lines(kr)%i j = Atom%lines(kr)%j @@ -211,7 +210,7 @@ subroutine get_Barklem_cross_data(atom, kr, res) write(*,*) "WARNING: Barklem, cannot parse label, using Unsold !" res = .false. return - endif + endif !write(*,'("Si="(1F2.2)", Li="(1I2)", Ji="(1F2.2))') Si, Li, Ji !write(*,'("Sj="(1F2.2)", Lj="(1I2)", Jj="(1F2.2))') Sj, Lj, Jj @@ -227,7 +226,7 @@ subroutine get_Barklem_cross_data(atom, kr, res) end if if (.not. res ) then - write(*,*) "WARNING: Barklem, cannot get data, using Unsold !" + write(*,*) "WARNING: Barklem, cannot get data, using Unsold !" return endif lread_from_table = .true. @@ -262,7 +261,7 @@ subroutine get_Barklem_cross_data(atom, kr, res) !write(*,*) neff2, bs%neff2(1), bs%neff2(bs%N2) write(*,*) "neff outside domain, use Unsold" res = .false. - return + return end if index = locate(bs%neff1,neff1)!minloc(abs(bs%neff1 - neff1), 1) @@ -293,15 +292,15 @@ subroutine get_Barklem_cross_data(atom, kr, res) atom%lines(kr)%cvdWaals(2) = & bilinear(bs%N1,bs%neff1,index1,bs%N2,bs%neff2,index,bs%alpha,neff1,neff2) - + endif !read from table !comptue the actual cross-section with T factorised out. cross = 2.0*sqa0 * (4.0/pi)**(0.5*atom%lines(kr)%cvdWaals(2)) * & - vmean * (vmean/vref)**(-atom%lines(kr)%cvdWaals(2)) + vmean * (vmean/vref)**(-atom%lines(kr)%cvdWaals(2)) - atom%lines(kr)%cvdWaals(1) = atom%lines(kr)%cvdWaals(1) * cross * & + atom%lines(kr)%cvdWaals(1) = atom%lines(kr)%cvdWaals(1) * cross * & exp(gammln(2.0-0.5*atom%lines(kr)%cvdWaals(2))) res = .true. @@ -309,4 +308,4 @@ subroutine get_Barklem_cross_data(atom, kr, res) return end subroutine get_barklem_cross_data -end module abo \ No newline at end of file +end module abo diff --git a/src/gas/atom_type.f90 b/src/gas/atom_type.f90 index 3a0f24db5..a618aff61 100644 --- a/src/gas/atom_type.f90 +++ b/src/gas/atom_type.f90 @@ -42,7 +42,7 @@ module atom_type real(kind=dp) :: twohnu3_c2, gij real(kind=dp) :: vmax, damp_max, damp_min !m/s real :: qwing - real(kind=dp), allocatable, dimension(:) :: Rij, Rji, Cij, Cji + real(kind=dp), allocatable, dimension(:) :: Rij, Rji, Cij, Cji real(kind=dp), dimension(4) :: cvdWaals real(kind=dp), dimension(:), allocatable :: a real(kind=dp), dimension(:,:,:,:,:), allocatable :: map !2d flux in the line @@ -61,7 +61,7 @@ module atom_type logical :: lany_gauss_prof integer :: initial, nTrans_raytracing !initial: - !0->LTE; 1->OLD_POPULATIONS; 2->ZERO_RADIATION; 3->CSWITCH; 4->SOBOLEV/CEP + !0->LTE; 1->OLD_POPULATIONS; 2->ZERO_RADIATION; 3->CSWITCH; 4->SOBOLEV/CEP !only for lines ! integer :: j_Trans_rayTracing(Nmax_Trans_raytracing), i_Trans_rayTracing(Nmax_Trans_raytracing) real :: vmax_rt !km/s !! @@ -108,12 +108,12 @@ module atom_type real(kind=dp), parameter :: cswitch_val = 1d6!1d12 real(kind=dp), parameter :: cswitch_down_scaling_factor = 10.0!ceil(exp(log(cswitch_val)/cswitch_niter)) logical :: lcswitch_enabled = .false. !if only one atom has cswith, it is True. - + contains function trans_number(atom, ilevel,jlevel) !from lower level i and upper level j - !return the transition number from the first line 1 + !return the transition number from the first line 1 !to the last continuum = atom%Ntrans ! ! @@ -121,7 +121,7 @@ function trans_number(atom, ilevel,jlevel) ! define an independent function not relying on atom%ij_to_trans type (AtomType), intent(in) :: atom integer, intent(in) :: ilevel, jlevel - integer :: i, j, k , trans_number + integer :: trans_number trans_number = atom%ij_to_trans(ilevel,jlevel) @@ -405,13 +405,13 @@ subroutine parse_label(label, g, S, L,J, determined) logical, intent(out) :: determined integer, intent(out) :: L real, intent(out) :: J, S - - + + logical :: unknown_orbital character(len=ATOM_LABEL_WIDTH+1) :: multiplet character(len=1) :: orbit integer :: i, Nl, parity, dk - + multiplet = trim(label) Nl = len(multiplet) @@ -519,16 +519,16 @@ function vbroad(temp, w, xi) real(kind=dp), intent(in) :: temp, xi real, intent(in) :: w real(kind=dp) :: vbroad - + vbroad = sqrt( Vtherm*temp/w + xi*xi ) - + return end function vbroad function ntotal_ions(icell) real(kind=dp) :: ntotal_ions integer, intent(in) :: icell - integer :: n, j, il + integer :: n, il ntotal_ions = 0.0 do n=1,n_atoms @@ -542,4 +542,4 @@ function ntotal_ions(icell) return end function ntotal_ions -end module atom_type \ No newline at end of file +end module atom_type diff --git a/src/gas/broad.f90 b/src/gas/broad.f90 index 77fcaea9d..6602b9eb8 100644 --- a/src/gas/broad.f90 +++ b/src/gas/broad.f90 @@ -18,11 +18,10 @@ function line_damping(icell, line) !for each transition individually integer, intent(in) :: icell type (AtomicLine), intent(in) :: line - integer :: k real(kind=dp) :: radfreq_to_vel real(kind=dp) :: adamp, vdw, vth real(kind=dp) :: line_damping - real(kind=dp) :: Qelast, Gj + line_damping = line%Grad vth = vbroad(T(icell), line%atom%weight, vturb(icell)) @@ -80,9 +79,9 @@ subroutine VanderWaals_line(icell, line, GvdW) real(kind=dp), intent(out) :: GvdW integer, intent(in) :: icell type (AtomicLine), intent(in) :: line - integer :: i, j, ic, Z, k + integer :: i, j, ic, Z real(kind=dp) :: vrel35_H, vrel35_He, fourPIeps0, deltaR - real(kind=dp) :: cross, gammaCorrH, gammaCorrHe, C625 + real(kind=dp) :: cross, C625 j = line%j i = line%i @@ -142,7 +141,7 @@ subroutine Stark_line(icell,line, GStark) real(kind=dp), intent(out) :: GStark integer, intent(in) :: icell type (AtomicLine), intent(in) :: line - integer :: k, ic, Z + integer :: ic, Z real(kind=dp) :: C4, C, Cm, melONamu, cStark23, cStark real(kind=dp) :: neff_u, neff_l, vrel, ERYDBERG @@ -195,9 +194,8 @@ subroutine StarkLinear_line(icell, line, GStark) real(kind=dp), intent(out) :: GStark integer, intent(in) :: icell type (AtomicLine), intent(in) :: line - real(kind=dp) :: n_upper, n_lower, Gstark2 - real :: Z, K, dn, a1, gs - real(kind=dp) :: C, nelectric + real(kind=dp) :: n_upper, n_lower, nelectric + real :: Z, dn, a1, gs Z = 1.0 + line%atom%stage(line%i) @@ -300,4 +298,4 @@ end subroutine StarkLinear_line ! return ! end subroutine resonance_broadening -end module broad \ No newline at end of file +end module broad diff --git a/src/gas/elements_type.f90 b/src/gas/elements_type.f90 index 62194701d..152915687 100644 --- a/src/gas/elements_type.f90 +++ b/src/gas/elements_type.f90 @@ -78,11 +78,9 @@ module elements_type 'Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es'/ - contains subroutine dealloc_elements() - integer :: n deallocate(Tpf) deallocate(elems) @@ -139,7 +137,7 @@ subroutine read_abundance() ! ! Also read partition function - integer :: EOF, n, blocksize, unit, i, j, syst_status + integer :: EOF, n, blocksize, unit, i, j integer :: NAXIST(1), naxis_found, hdutype, Nread ! integer :: cursor_init character(len=256) :: some_comments diff --git a/src/gas/gas_contopac.f90 b/src/gas/gas_contopac.f90 index 1eb5cca45..53622244a 100644 --- a/src/gas/gas_contopac.f90 +++ b/src/gas/gas_contopac.f90 @@ -20,7 +20,7 @@ module gas_contopac use messages, only : error implicit none - + real, parameter :: lambda_limit_HI_rayleigh = 91.1763062831680!102.6 !nm real, parameter :: lambda_limit_HeI_rayleigh = 140.0 ! real(kind=dp), dimension(:), allocatable :: Hray_lambda, HeIray_lambda @@ -30,7 +30,7 @@ module gas_contopac real(kind=dp), dimension(N_geltman) :: lambdai_geltman, alphai_geltman real(kind=dp), dimension(:,:), allocatable :: john_hminus_term2, john_hminus_term1 real(kind=dp) :: john_an(6), john_an2(4), john_bn(6), john_bn2(4),john_cn(6), john_cn2(4), & - john_dn(6), john_dn2(4), john_en(6), john_en2(4), john_fn(6), john_fn2(4) + john_dn(6), john_dn2(4), john_en(6), john_en2(4), john_fn(6), john_fn2(4) real(kind=dp), dimension(:), allocatable :: j0_theta_bell_berr integer :: i0_lam_bell_berr real(kind=dp), allocatable :: alphai_bell_berr_part(:,:) @@ -217,7 +217,7 @@ function Thomson(icell) end function Thomson function HI_rayleigh(icell) - integer, intent(in) :: icell + integer, intent(in) :: icell real(kind=dp) :: HI_rayleigh(size(Hray_lambda)) HI_rayleigh = sigma_e * Hray_lambda(:) * Hydrogen%n(1,icell) @@ -240,7 +240,7 @@ subroutine alloc_gas_contopac(N, lambda) !could be para, at least the cell part. integer, intent(in) :: N real(kind=dp), intent(in) :: Lambda(N) - real(kind=dp), allocatable, dimension(:) :: tab_lambda_ang + real(kind=dp), dimension(N) :: tab_lambda_ang real(kind=dp) :: theta, lam integer :: icell, la, i @@ -256,7 +256,7 @@ subroutine alloc_gas_contopac(N, lambda) (148.d0/lambda)**4d0)*(96.6d0/lambda)**4d0 elsewhere Hray_lambda = 0.0 - end where + end where if (associated(helium)) then allocate(HeIray_lambda(N)) @@ -265,7 +265,7 @@ subroutine alloc_gas_contopac(N, lambda) (64.1d0/lambda)**4d0)*(37.9d0/lambda)**4d0 elsewhere HeIray_lambda = 0.0 - end where + end where endif !cheap ! @@ -283,7 +283,6 @@ subroutine alloc_gas_contopac(N, lambda) do i=1,11 alphai_bell_berr_part(:,i) = linear_1D_sorted(23,lambdai_bell_berr,alphai_bell_berr(:,i),n,tab_lambda_ang) enddo - deallocate(tab_lambda_ang) i0_lam_bell_berr = max(locate(lambdai_bell_berr, 10*lambda(1)),2)!all wavelengths are contiguous. We only need the index of the first one !we only need term1 of John H- free-free since, bell_berr routine is used. @@ -314,7 +313,7 @@ subroutine alloc_gas_contopac(N, lambda) ! write(*,'("allocate "(1F6.4)" GB for exp(-hc/kT)")') real(sizeof(exphckT))/1024./1024./1024. ! write(*,'(" -> max(ehnukt)="(1ES14.5E3)"; min(ehnukt)="(ES14.5E3))') maxval(exphckT(:)), minval(exphckT(:),mask=icompute_AtomRT>0) - return + return end subroutine alloc_gas_contopac subroutine dealloc_gas_contopac @@ -335,7 +334,7 @@ subroutine dealloc_gas_contopac return end subroutine dealloc_gas_contopac - + elemental function Gaunt_bf(u, n_eff) ! M. J. Seaton (1960), Rep. Prog. Phys. 23, 313 ! See also Menzel & Pekeris 1935 @@ -381,7 +380,7 @@ end function Gaunt_ff elemental function H_bf_Xsection(cont, lambda) result(alpha) Type (AtomicContinuum), intent(in) :: cont real(kind=dp), intent(in) :: lambda - real(kind=dp) :: n_eff, g_bf, u, Z, u0, g_bf0, alpha, u1 + real(kind=dp) :: n_eff, g_bf, u, Z, alpha Z = real(cont%atom%stage(cont%i) + 1,kind=dp) @@ -411,11 +410,11 @@ elemental function H_bf_Xsection(cont, lambda) result(alpha) return end function H_bf_Xsection - + elemental function H_ff_Xsection(Z, T, lambda) result(alpha) real(kind=dp), intent(in) :: lambda, T - real(kind=dp) :: g_ff, nu3, alpha, K0 + real(kind=dp) :: g_ff, nu3, alpha integer, intent(in) :: Z ! K0 = (Q_ELECTRON**2)/(4.0*PI*EPSILON_0) / sqrt(M_ELECTRON) @@ -440,8 +439,7 @@ subroutine Hydrogen_ff(icell, N, lambda, chi) integer, intent(in) :: icell, N real(kind=dp), intent(in) :: lambda(N) real(kind=dp), intent(out) :: chi(N) - integer :: la - real(kind=dp) :: stim, np, arg_exp, exp_val, C0, alpha + real(kind=dp) :: np np = Hydrogen%n(Hydrogen%Nlevel,icell) !nH+=Nion H !should be nstar instead of %n ? @@ -470,7 +468,7 @@ subroutine Hydrogen_ff(icell, N, lambda, chi) return end subroutine Hydrogen_ff - + subroutine atom_ff return @@ -489,7 +487,7 @@ subroutine Hminus_bf_john(icell, N, lambda, chi, eta) real(kind=dp) :: lam, lambda0, alpha, sigma, flambda, Cn(6) real(kind=dp) :: diff, stm, funit, cte, pe integer :: la - real(kind=dp) :: arg_exp, nH + real(kind=dp) :: nH real(kind=dp), dimension(N), intent(out) :: chi, eta chi(:) = 0.0_dp @@ -644,7 +642,7 @@ elemental function Hminus_ff_john_lam(icell, lambda) result(chi) !tab 3b, lambda > 0.1823 micron and < 0.3645 micron, size of 4 isntead 5, because the last two rows are 0 real(kind=dp), dimension(4) :: An2, Bn2, Cn2, Dn2, En2, Fn2 real(kind=dp) :: funit, K, kappa - integer :: la, n + integer :: n real(kind=dp) :: chi real(kind=dp) :: lam, theta, nH @@ -700,9 +698,8 @@ subroutine Hminus_ff_john(icell, N, lambda, chi) !----------------------------------------------------------------- integer, intent(in) :: icell, N real(kind=dp), dimension(N), intent(in) :: lambda - integer :: la real(kind=dp), dimension(N), intent(out) :: chi - real(kind=dp) :: lam, theta!, nH + real(kind=dp) :: theta chi(:) = 0.0_dp @@ -761,7 +758,7 @@ subroutine Hminus_ff_bell_berr(icell, N, lambda, chi) integer :: la real(kind=dp), dimension(N), intent(out) :: chi integer :: i0, j0 - real(kind=dp) :: lam, stm, sigma(N), theta, pe, nH, tt + real(kind=dp) :: lam, sigma(N), theta, pe, nH, tt chi(:) = 0.0_dp theta = 5040. / T(icell) @@ -807,4 +804,4 @@ subroutine Hminus_ff_bell_berr(icell, N, lambda, chi) end subroutine Hminus_ff_bell_berr -end module gas_contopac \ No newline at end of file +end module gas_contopac diff --git a/src/gas/lte.f90 b/src/gas/lte.f90 index 667a8daf4..28f772e40 100644 --- a/src/gas/lte.f90 +++ b/src/gas/lte.f90 @@ -67,7 +67,6 @@ function BoltzmannEq4dot20b(temp, Ei, gi, gi1) real(kind=dp) :: BoltzmannEq4dot20b real(kind=dp) :: kT real(kind=dp), intent(in) :: gi, gi1, Ei, temp - integer :: k kT = KB * temp BoltzmannEq4dot20b = (gi1/gi) * exp(-Ei/kT) !ni1/ni @@ -80,7 +79,7 @@ end function BoltzmannEq4dot20b function nH_minus(icell) !Saha Ionisation equation applied to H- integer :: icell - real(kind=dp) :: nH_minus, nH, UH, UHm, n00 + real(kind=dp) :: nH_minus, nH, UH, UHm UH = 2.0 !UH = get_pf(elements(1), 1, T(icell)) @@ -106,10 +105,8 @@ subroutine LTEpops_H_loc (k) ! TO DO: introduce nH- ! -------------------------------------------------------------- ! integer, intent(in) :: k - logical :: locupa_prob - real(kind=dp) :: dEion, dE, sum, c2, phik, phiHmin, ntotal - real(kind=dp) :: n_eff, wocc, chi0, wocc0, E00, E, Egs - integer :: Z, dZ, i, m + real(kind=dp) :: dE, sum, phik, ntotal, wocc, chi0, E00, E, Egs + integer :: dZ, i, m E00 = 1.0 * 3e-11 * electron_charge ! Joules Egs = hydrogen%E(1) @@ -194,7 +191,7 @@ subroutine LTEpops_H_loc (k) else hydrogen%ni_on_nj_star(hydrogen%continua(i)%i,k) = 0.0 endif - enddo + enddo return end subroutine LTEpops_H_loc @@ -209,10 +206,9 @@ subroutine LTEpops_atom_loc(k,atom,debye) integer, intent(in) :: k type (Atomtype), intent(inout) :: atom logical, intent(in) :: debye - logical :: locupa_prob, print_diff - real(kind=dp) :: dEion, dE, sum, c2, phik, phiHmin - real(kind=dp) :: n_eff, wocc, chi0, wocc0, ntotal - integer :: Z, dZ, i, m + logical :: locupa_prob + real(kind=dp) :: dEion, dE, sum, c2, phik, n_eff, wocc, ntotal + integer :: dZ, i, m ! debye shielding activated: ! It lowers the ionisation potential taking into account the charge of all levels @@ -343,7 +339,7 @@ subroutine LTEpops_atom_loc(k,atom,debye) else atom%ni_on_nj_star(atom%continua(i)%i,k) = 0.0 endif - enddo + enddo return end subroutine LTEpops_atom_loc @@ -429,7 +425,7 @@ subroutine LTEpops_H() logical :: locupa_prob real(kind=dp) :: dEion, dE, sum, c2, phik, phiHmin real(kind=dp) :: n_eff, wocc, chi0, wocc0, E00, E, Egs - integer :: Z, dZ, k, i, m + integer :: dZ, k, i, m logical :: print_diff real(kind=dp) :: max_nstar(hydrogen%Nlevel), min_nstar(hydrogen%Nlevel) @@ -576,8 +572,8 @@ subroutine LTEpops_atom(atom, debye) logical, intent(in) :: debye logical :: locupa_prob, print_diff real(kind=dp) :: dEion, dE, sum, c2, phik, phiHmin - real(kind=dp) :: n_eff, wocc, chi0, wocc0 - integer :: Z, dZ, k, i, m + real(kind=dp) :: n_eff, wocc + integer :: dZ, k, i, m real(kind=dp), dimension(atom%Nlevel) :: max_nstar, min_nstar !old_values ! debye shielding activated: @@ -906,7 +902,6 @@ subroutine write_ltepops_file(unit, atom, delta_k) integer, intent(in) :: unit integer, intent(in), optional :: delta_k integer :: dk, k, i, j, Nstage, Nj, ip - real(kind=dp) :: wocc if (present(delta_k)) then dk = delta_k diff --git a/src/gas/opacity_atom.f90 b/src/gas/opacity_atom.f90 index 1185d9f9b..93ac75ea7 100644 --- a/src/gas/opacity_atom.f90 +++ b/src/gas/opacity_atom.f90 @@ -90,7 +90,7 @@ subroutine deactivate_lines() atoms(nact)%p%lines(kr)%lcontrib = .false. enddo enddo - return + return end subroutine deactivate_lines subroutine activate_lines() integer :: nact, kr @@ -99,7 +99,7 @@ subroutine activate_lines() atoms(nact)%p%lines(kr)%lcontrib = .true. enddo enddo - return + return end subroutine activate_lines subroutine deactivate_continua() integer :: nact, kr @@ -108,7 +108,7 @@ subroutine deactivate_continua() atoms(nact)%p%continua(kr)%lcontrib = .false. enddo enddo - return + return end subroutine deactivate_continua subroutine activate_continua() integer :: nact, kr @@ -117,7 +117,7 @@ subroutine activate_continua() atoms(nact)%p%continua(kr)%lcontrib = .true. enddo enddo - return + return end subroutine activate_continua !could be parralel @@ -137,6 +137,7 @@ subroutine alloc_atom_opac(N,x,limage) if (limit_mem > 0) then !computed and stored on tab_lambda_cont and interpolated locally on tab_lambda_nm !-> on a small grid and interpolated later np = n_lambda_cont + allocate(xp(np)) xp = tab_lambda_cont if (limit_mem==1) then contopac_atom_loc => contopac_atom_loc_mem1 @@ -149,6 +150,7 @@ subroutine alloc_atom_opac(N,x,limage) !a linear interpolation is faster than the evaluation of the continua (and they are mostly flat...). else !computed and stored on tab_lambda_nm np = n + allocate(xp(np)) xp = x contopac_atom_loc => contopac_atom_loc_mem0 endif @@ -336,8 +338,8 @@ subroutine contopac_atom_loc_mem0(icell,N,lambda,chi,snu) chi = chi_cont(:,icell) !n_lambda snu = eta_cont(:,icell) !n_lambda - return - + return + end subroutine contopac_atom_loc_mem0 subroutine contopac_atom_loc_mem1(icell,N,lambda,chi,snu) @@ -378,7 +380,7 @@ subroutine contopac_atom_loc_mem2(icell,N,lambda,chi,snu) call background_continua_lambda(icell, n_lambda_cont, tab_lambda_cont, chic, Snuc) !Snu = Snu + scat(lambda, icell) * Jnu(:,icell) !accumulate b-f - call opacity_atom_bf_loc(icell, n_lambda_cont, tab_lambda_cont, chic, Snuc) + call opacity_atom_bf_loc(icell, n_lambda_cont, tab_lambda_cont, chic, Snuc) !linear interpolation from tab_lambda_cont to lambda i0 = 2 @@ -489,16 +491,16 @@ subroutine opacity_atom_bb_loc(id, icell, iray, x, y, z, x1, y1, z1, u, v, w, l_ vth = vbroad(T(icell),Atom%weight, vturb(icell)) - !any line of that atom with a Gaussian profile ? + !any line of that atom with a Gaussian profile ? !use the same profile for all !TO DO: before each mode (nlte, image/flux) if (atom%lany_gauss_prof) then - if (lnon_lte_loop) then + if (lnon_lte_loop) then !completely linear for Gaussian with the sum of Ngauss_log and Ngauss_lin (see line_lambda_grid in gas/wavelengths_gas.f90) nvel_rt = N_gauss peak_g = peak_gauss_limit / sqrt(pi) / vth vg_max = vth * sqrt(-log(peak_g * sqrt(pi) * vth)) - else + else !image mode or flux, linear grid in hv nvel_rt = atom%n_speed_rt vg_max = 1d3 * atom%vmax_rt @@ -698,7 +700,7 @@ subroutine cross_coupling_cont(id,icell,at) ! at%continua(kr)%alpha(:) * at%continua(kr)%twohnu3_c2(:) * ni_on_nj_star * exp(hnu_k(Nb:Nr)/T(icell)) ! eta_atoms(Nb:Nr,nact,id) = eta_atoms(Nb:Nr,nact,id) + & ! at%continua(kr)%alpha(:) * at%continua(kr)%twohnu3_c2(:) * ni_on_nj_star * exp(hnu_k(Nb:Nr)/T(icell)) * at%n(j,icell) - + enddo cont_loop return @@ -710,7 +712,7 @@ subroutine cross_coupling_cont_i(id,icell,at) type(AtomType), intent(in), pointer :: at integer :: nact, j, i, kr, Nb, Nr, la, Nl integer :: i0, la0, N1, N2 - real(kind=dp) :: gij, chicc, wl, ni_on_nj_star, wt + real(kind=dp) :: gij, wl, ni_on_nj_star, wt real(kind=dp) :: term1(Nlambda_max_cont), term2(Nlambda_max_cont), term3(Nlambda_max_cont) nact = at%activeindex @@ -899,7 +901,7 @@ function profile_art_i(line,id,icell,iray,lsubstract_avg,N,lambda, x,y,z,x1,y1,z l_void_before,l_contrib !physical length of the cell integer :: Nvspace real(kind=dp), dimension(NvspaceMax) :: Omegav - real(kind=dp) :: norm, vth + real(kind=dp) :: vth real(kind=dp) :: v0, v1, delta_vol_phi, xphi, yphi, zphi, & dv, omegav_mean type (AtomicLine), intent(in) :: line @@ -1024,7 +1026,7 @@ function gauss_profile(id,icell,iray,lsubstract_avg,N,vel,vth,x,y,z,x1,y1,z1,u,v if (lnon_lte_loop) omegav(1:Nvspace) = omegav(1:Nvspace) - vlabs(iray,id) endif - + u0sq(:) = u0(:)*u0(:) !Note: u1 = (u0 - omegav(nv)/vth)**2 u1(:) = u0sq(:) + (omegav(1)/vth)*(omegav(1)/vth) - 2*u0(:) * omegav(1)/vth @@ -1045,12 +1047,11 @@ subroutine write_opacity_emissivity_bin(Nlambda,lambda) integer, intent(in) :: Nlambda real(kind=dp), intent(in) :: lambda(Nlambda) integer :: unit, unit2, status = 0 - integer :: alloc_status, id, icell, m, Nrec - real(kind=dp), allocatable, dimension(:,:,:) :: chi_tmp, eta_tmp, rho_tmp + integer :: alloc_status, id, icell, Nrec + real(kind=dp), allocatable, dimension(:,:,:) :: chi_tmp, eta_tmp real(kind=dp), allocatable, dimension(:,:,:) :: chic_tmp, etac_tmp character(len=11) :: filename_chi="chi.b" character(len=50) :: filename_eta="eta.b" - character(len=18) :: filename_rho="magnetoopt.b" real(kind=dp) :: zero_dp, u, v, w zero_dp = 0.0_dp diff --git a/src/gas/wavelengths_gas.f90 b/src/gas/wavelengths_gas.f90 index ad6cfa0ab..4b2d195f3 100644 --- a/src/gas/wavelengths_gas.f90 +++ b/src/gas/wavelengths_gas.f90 @@ -61,8 +61,7 @@ function make_sub_wavelength_grid_cont_log_nu(cont) ! real(kind=dp), dimension(Nlambda_cont_log) :: make_sub_wavelength_grid_cont_log_nu type (AtomicContinuum), intent(in) :: cont real(kind=dp), dimension(cont%Nlambda) :: make_sub_wavelength_grid_cont_log_nu - real(kind=dp) :: resol - integer :: la, N1, N2 + integer :: N1, N2 real(kind=dp) :: l0, l1 real(kind=dp) :: nu1, nu0, nu2, dnu @@ -98,7 +97,7 @@ subroutine compute_line_bound(line,limage) ! ------------------------------------------------------------ ! ! Compute the line bounds : lamndamin and lambdamax ! the total extent of line in the atom's frame (no velocity). - ! if limage is .false. : + ! if limage is .false. : ! vmax is computed from the maximum thermal width of the line. ! and from the "natural" width of the line. ! if limage is .true : @@ -170,7 +169,7 @@ function line_lambda_grid(line,Nlambda) real(kind=dp), dimension(Nlambda) :: line_lambda_grid type (atomicline), intent(inout) :: line real(kind=dp) :: vcore, vwing, vth_max - integer :: la, Nmid + integer :: Nmid integer :: N_wing, N_core @@ -244,7 +243,7 @@ function line_lambda_grid_dv(line,Nlambda) real(kind=dp), dimension(Nlambda) :: line_lambda_grid_dv type (atomicline), intent(inout) :: line real(kind=dp) :: vwing - integer :: la, Nmid + integer :: la !not the exact width at a given cell, but the maximum extent! vwing = line%vmax @@ -261,7 +260,7 @@ function line_lambda_grid_dv(line,Nlambda) end function line_lambda_grid_dv subroutine deallocate_wavelengths_gasrt()!(lambda) - use wavelengths, only : n_lambda, tab_lambda, tab_lambda_inf, tab_lambda_sup, tab_delta_lambda, n_lambda2, tab_lambda2 + use wavelengths, only : tab_lambda, tab_lambda_inf, tab_lambda_sup, tab_delta_lambda ! real(kind=dp), intent(inout), allocatable, dimension(:) :: lambda if (allocated(tab_lambda)) deallocate(tab_lambda) !lambda if (allocated(tab_lambda_inf)) then @@ -284,12 +283,12 @@ subroutine make_wavelengths_nlte(lambda,Vmax_overlap) integer :: Ntrans, Ncont, Nlines, Nlambda_cont, Nlambda integer :: n, kr, Nlam, Nmore_cont_freq, Nremoved, Nwaves, check_new_freq type (AtomType), pointer :: atom - integer :: alloc_status, lac, la, nb, krr, max_Nlambda_indiv, Noverlap, ib, ir + integer :: alloc_status, lac, la, max_Nlambda_indiv, Noverlap, ib, ir real(kind=dp), dimension(:), allocatable :: tmp_grid, tmp_grid2, all_l0, all_l1 integer, parameter :: Ngroup_max = 1000 real(kind=dp), dimension(Ngroup_max) :: group_blue_tmp, group_red_tmp, Nline_per_group_tmp integer, dimension(:), allocatable :: sorted_indexes - real(kind=dp) :: max_cont, l0, l1, dvmin, min_lines, min_cont + real(kind=dp) :: max_cont, l0, l1, min_lines, min_cont real :: hv_loc !like hv real, parameter :: prec_vel = 0.0 !remove points is abs(hv_loc-hv)>prec_vel logical :: check_for_overlap @@ -613,7 +612,7 @@ subroutine make_wavelengths_nlte(lambda,Vmax_overlap) Nmore_cont_freq = Nmore_cont_freq + 1 endif - !and below first bound-free if the bluest line wavelength is below + !and below first bound-free if the bluest line wavelength is below !the lower continuum wavelength ! if (l0 < min_cont) then Nmore_cont_freq = Nmore_cont_freq + 1 @@ -1247,7 +1246,7 @@ subroutine make_wavelengths_flux(lambda,lfrom_file) integer :: Ntrans, Ncont, Nlines, Nlambda_cont, Nlambda integer :: n, kr, Nlam, Nmore_cont_freq, Nremoved, Nwaves, check_new_freq type (AtomType), pointer :: atom - integer :: alloc_status, lac, la, nb, krr + integer :: alloc_status, lac, la real(kind=dp), dimension(:), allocatable :: tmp_grid, tmp_grid2, all_l0, all_l1 integer, parameter :: Ngroup_max = 1000 real(kind=dp), dimension(Ngroup_max) :: group_blue_tmp, group_red_tmp, Nline_per_group_tmp diff --git a/src/input.f90 b/src/input.f90 index d6e0c229b..0381802a7 100644 --- a/src/input.f90 +++ b/src/input.f90 @@ -727,7 +727,7 @@ subroutine read_phase_function(phase_function_file) character(len=512), intent(in) :: phase_function_file - integer :: status, readwrite, unit, blocksize,nfound,group,firstpix,nbuffer,npixels, hdutype, i + integer :: status, readwrite, unit, blocksize,nfound,group,firstpix,nbuffer,npixels, i real :: nullval integer, dimension(2) :: naxes logical :: anynull diff --git a/src/mcfost_env.f90 b/src/mcfost_env.f90 index de65a7d82..922e8047d 100644 --- a/src/mcfost_env.f90 +++ b/src/mcfost_env.f90 @@ -5,7 +5,7 @@ module mcfost_env implicit none real, parameter :: mcfost_version = 4.1 - character(8), parameter :: mcfost_release = "4.1.07" + character(8), parameter :: mcfost_release = "4.1.08" real, parameter :: required_utils_version = 4.0 !character(len=128) :: web_server = "https://ipag.osug.fr/public/pintec/" diff --git a/src/molecular_emission.f90 b/src/molecular_emission.f90 index df982e896..1c9063e9b 100644 --- a/src/molecular_emission.f90 +++ b/src/molecular_emission.f90 @@ -687,7 +687,7 @@ function v_proj(icell,x,y,z,u,v,w) ! real(kind=dp), intent(in) :: x,y,z,u,v,w real(kind=dp) :: vitesse, vx, vy, vz, v_r, v_phi, v_theta, v_rcyl, norme, r, phi, rcyl, rcyl2, r2 - real(kind=dp) :: cos_phi, sin_phi, cos_theta, sin_theta, norme2 + real(kind=dp) :: cos_phi, sin_phi, cos_theta, sin_theta if (lVoronoi) then vx = Voronoi(icell)%vxyz(1) diff --git a/src/readVTK.f90 b/src/readVTK.f90 index 114d537bb..b70376816 100644 --- a/src/readVTK.f90 +++ b/src/readVTK.f90 @@ -60,7 +60,6 @@ subroutine readField(file, oldposition, geometry, periodicity, time, newposition integer, intent(in) :: file, oldposition integer, intent(out) :: geometry integer, dimension(3), intent(out) :: periodicity - integer geo1, geo2 real, intent(out) :: time integer, intent(out) :: newposition @@ -113,7 +112,6 @@ subroutine readPoints(file, oldposition, geometry, dimensions, x1, x2, x3, newpo real, allocatable :: pts(:) real, dimension(:,:,:), allocatable :: x, y, z - real q1, q2; integer position, lineSize, nPoints integer i,j,k,n, nx1mid, nx2mid, nx3mid @@ -250,9 +248,8 @@ subroutine readVTK_header(filename, unit, position, geometry, dimensions, time, real, dimension(:), allocatable, intent(out) :: x1, x2, x3 character (:), allocatable :: line - character (3) :: varName integer :: position, newposition - integer :: lineSize, nPoints + integer :: lineSize integer, dimension(3) :: periodicity open(newunit=unit, file=filename,form="unformatted",access="stream",CONVERT='BIG_ENDIAN',status="old") diff --git a/src/read_fargo3d.f90 b/src/read_fargo3d.f90 index 3018d670f..15a7f1a75 100644 --- a/src/read_fargo3d.f90 +++ b/src/read_fargo3d.f90 @@ -202,6 +202,7 @@ subroutine read_fargo3d_files() file_types(2) = "gasvx" file_types(3) = "gasvy" file_types(4) = "gasvz" + iunit=1 do l=1, 4 filename = trim(fargo3d%dir)//"/"//trim(file_types(l))//trim(trim(fargo3d%id))//".dat" write(*,*) "Reading "//trim(filename) @@ -304,7 +305,7 @@ subroutine read_fargo3d_planets(dir, n_planets,x,y,z,vx,vy,vz,Mp,time,Omega_p) integer, intent(out) :: n_planets real(dp), dimension(n_planets_max), intent(out) :: x, y, z, vx, vy, vz, Mp, Omega_p, time - integer :: n_etoiles_old, iunit, ios, n_etoile_old, i, i_planet, id + integer :: iunit, ios, i, i_planet, id character(len=1) :: s character(len=128) :: filename diff --git a/src/read_idefix.f90 b/src/read_idefix.f90 index 153b51712..09fc0903e 100644 --- a/src/read_idefix.f90 +++ b/src/read_idefix.f90 @@ -21,11 +21,10 @@ subroutine read_idefix_parameters(filename) character(len=128) :: name character(len=1) :: sep - integer :: geometry, time, pos1,pos2, id + integer :: pos1,pos2 integer, dimension(3) :: dimensions character(:), allocatable :: origin real, dimension(:), allocatable :: x1, x2, x3 - real(dp) :: dx write(*,*) "Reading header: "//trim(filename) idefix%filename = filename @@ -105,14 +104,8 @@ subroutine read_idefix_model() ! idefix data is ordered in x1 = r, x2 = theta, x3 = phi real, dimension(:,:,:), allocatable :: rho, vx1, vx2, vx3 - integer :: ios, iunit, alloc_status, l, recl, i,j, jj, phik, icell, id, n_planets, n_etoiles_old, i2, i3 - - character(len=128) :: filename - character(len=16), dimension(4) :: file_types - + integer :: alloc_status, i,j, jj, phik, icell, n_planets, i2, i3 real(dp) :: Ggrav_idefix, umass, usolarmass, ulength, utime, udens, uvelocity, ulength_au, mass, facteur, omega - type(star_type), dimension(:), allocatable :: etoile_old - real(dp), dimension(n_planets_max) :: x, y, z, vx, vy, vz, Mp, Omega_p, time ! Planet properties hard coded for now @@ -252,7 +245,7 @@ subroutine read_idefix_planets(dir, id, n_planets,x,y,z,vx,vy,vz,Mp,time) integer, intent(out) :: n_planets real(dp), dimension(n_planets_max), intent(out) :: x, y, z, vx, vy, vz, Mp, time - integer :: n_etoiles_old, iunit, ios, n_etoile_old, i, i_planet + integer :: iunit, ios, i, i_planet real :: timestep character(len=1) :: s diff --git a/src/read_param.f90 b/src/read_param.f90 index dfa0e7719..458cedec7 100644 --- a/src/read_param.f90 +++ b/src/read_param.f90 @@ -4551,7 +4551,7 @@ subroutine read_para211(para) character(len=*), intent(in) :: para - integer :: i, j, alloc_status, ios, tmpint, ind_pop, status, imol + integer :: i, j, alloc_status, ios, ind_pop, status, imol real(kind=dp) :: size_neb_tmp, somme real :: gas_dust diff --git a/src/read_phantom.f90 b/src/read_phantom.f90 index 5e9112a08..353768413 100644 --- a/src/read_phantom.f90 +++ b/src/read_phantom.f90 @@ -46,7 +46,7 @@ subroutine read_phantom_bin_files(iunit,n_files, filenames, x,y,z,h,vx,vy,vz,T_g real(dp), allocatable, dimension(:,:) :: xyzh,xyzmh_ptmass,vxyz_ptmass,dustfrac,vxyzu,nucleation type(dump_h) :: hdr logical :: got_h,got_dustfrac,got_itype,tagged,matched - logical :: got_temperature,got_u,lpotential,do_nucleation + logical :: got_temperature,got_u,lpotential integer :: ifile, np0, ntypes0, np_tot, ntypes_tot, ntypes_max, ndustsmall, ndustlarge ! We first read the number of particules in each phantom file diff --git a/src/read_pluto.f90 b/src/read_pluto.f90 index 925b5bd05..682c103a1 100644 --- a/src/read_pluto.f90 +++ b/src/read_pluto.f90 @@ -186,6 +186,7 @@ subroutine read_pluto_files() file_types(2) = "vx1" file_types(3) = "vx2" file_types(4) = "vx3" + iunit=1 do l=1, 4 filename = trim(pluto%dir)//"/"//trim(file_types(l))//"."//trim(trim(pluto%id))//".dbl" write(*,*) "Reading "//trim(filename) diff --git a/src/read_spherical_grid.f90 b/src/read_spherical_grid.f90 index da98494d2..a22fefd83 100644 --- a/src/read_spherical_grid.f90 +++ b/src/read_spherical_grid.f90 @@ -32,7 +32,7 @@ subroutine read_spherical_grid_parameters(filename) ! in pluto%x1, pluto%x2, pluto%x3. ! ----------------------------------------------------- ! character(len=*), intent(in) :: filename - integer :: ios, Nsize, acc, i, ipos + integer :: ios, Nsize, acc, i real :: dphi lvelocity_file = .true. @@ -42,7 +42,7 @@ subroutine read_spherical_grid_parameters(filename) n_zones = 1 open(unit=1, file=trim(filename), status="old",access="stream",form='unformatted') - !read size along each direction + cell limits (size + 1) + !read size along each direction + cell limits (size + 1) read(1,iostat=ios) pluto%nx1 allocate(pluto%x1(pluto%nx1+1)) read(1,iostat=ios) pluto%x1(:) @@ -138,7 +138,7 @@ subroutine read_spherical_model(filename) ! use grains, only : M_Grain character(len=*), intent(in) :: filename integer :: ios, i, Nsize - + integer :: j, k, jj, icell integer, allocatable :: dz(:,:,:) real, allocatable :: vtmp(:,:,:,:) @@ -271,4 +271,4 @@ subroutine read_spherical_model(filename) return endsubroutine read_spherical_model -end module read_spherical_grid \ No newline at end of file +end module read_spherical_grid diff --git a/src/scattering.f90 b/src/scattering.f90 index 28426cd00..e86b3b9ae 100644 --- a/src/scattering.f90 +++ b/src/scattering.f90 @@ -387,7 +387,7 @@ subroutine Mueller_input(lambda, k_abs,k_sca,gsca) character(len=128) :: string, filename - integer :: EoF, alloc_status,iformat,nlam,nmu,ilam,iline,iang,nang, l + integer :: EoF, alloc_status,iformat,nlam,iang,nang, l logical :: lscat real(dp), dimension(:), allocatable :: angles,wavel,kabs,ksca,g real(dp), dimension(:,:), allocatable :: f11,f12,f22,f33,f34,f44 @@ -426,7 +426,7 @@ subroutine Mueller_input(lambda, k_abs,k_sca,gsca) allocate(wavel(nlam),kabs(nlam),ksca(nlam),g(nlam), stat = alloc_status) ! Read wavelengths, opacities and g - do ilam=1,nlam + do l=1,nlam read(1,*) wavel(l), kabs(l), ksca(l), g(l) enddo @@ -435,14 +435,14 @@ subroutine Mueller_input(lambda, k_abs,k_sca,gsca) if (wl > maxval(wavel)) call error("Mueller matrices: wavelength is larger than tabulated") if (wavel(2) > wavel(1)) then ! increasing order - do ilam=1,nlam + do l=1,nlam if (wavel(l) > wl) then frac = (wavel(l) - wl) / (wavel(l) - wavel(l-1)) exit endif enddo else ! decreasing order - do ilam=1,nlam + do l=1,nlam if (wavel(l) < wl) then frac = (wl - wavel(l)) / (wavel(l-1) - wavel(l)) exit @@ -450,6 +450,7 @@ subroutine Mueller_input(lambda, k_abs,k_sca,gsca) enddo endif if ((frac > 1) .or. (frac < 0)) call error("Mueller matrices: interpolation error") + frac_m1 = 1.0_dp - frac if (lscat) then allocate(angles(nang), f11(nlam,nang),f12(nlam,nang),f22(nlam,nang),& @@ -461,27 +462,28 @@ subroutine Mueller_input(lambda, k_abs,k_sca,gsca) enddo ! Read the scattering matrix. - do ilam=1,nlam + do l=1,nlam do iang=1,nang - read(1,*) f11(ilam,iang),f12(ilam,iang), f22(ilam,iang), f33(ilam,iang), f34(ilam,iang), f44(ilam,iang) + read(1,*) f11(l,iang),f12(l,iang), f22(l,iang), f33(l,iang), f34(l,iang), f44(l,iang) enddo enddo endif close(1) ! Interpolate at correct wavelength - ! todo : find ilam + ! todo : find l + call error("Mueller input: interpolation was never implemented") - k_abs = kabs(ilam-1) * frac + kabs(ilam) * frac_m1 - k_sca = ksca(ilam-1) * frac + ksca(ilam) * frac_m1 - gsca = ksca(ilam-1) * frac + ksca(ilam) * frac_m1 + k_abs = kabs(l-1) * frac + kabs(l) * frac_m1 + k_sca = ksca(l-1) * frac + ksca(l) * frac_m1 + gsca = ksca(l-1) * frac + ksca(l) * frac_m1 - s11(:) = f11(ilam-1,:) * frac + f11(ilam,:) * frac_m1 - s12(:) = f11(ilam-1,:) * frac + f12(ilam,:) * frac_m1 - s22(:) = f11(ilam-1,:) * frac + f22(ilam,:) * frac_m1 - s33(:) = f11(ilam-1,:) * frac + f33(ilam,:) * frac_m1 - s34(:) = f11(ilam-1,:) * frac + f34(ilam,:) * frac_m1 - s44(:) = f11(ilam-1,:) * frac + f44(ilam,:) * frac_m1 + s11(:) = f11(l-1,:) * frac + f11(l,:) * frac_m1 + s12(:) = f11(l-1,:) * frac + f12(l,:) * frac_m1 + s22(:) = f11(l-1,:) * frac + f22(l,:) * frac_m1 + s33(:) = f11(l-1,:) * frac + f33(l,:) * frac_m1 + s34(:) = f11(l-1,:) * frac + f34(l,:) * frac_m1 + s44(:) = f11(l-1,:) * frac + f44(l,:) * frac_m1 ! Deallocate arrays deallocate(wavel,kabs,ksca,g) @@ -507,7 +509,7 @@ subroutine normalise_Mueller_matrix(lambda,taille_grain, s11,s12,s22,s33,s34,s44 ! as we assume it is mostly diffraction that has been missed by the sampling integer :: j - real(kind=dp) :: norm, somme_prob, theta, dtheta + real(kind=dp) :: norm, theta, dtheta ! Integration S11 pour tirer angle if (scattering_method==1) then @@ -1213,7 +1215,7 @@ subroutine update_Stokes(S,u0,v0,w0,u1,v1,w1, M) real(kind=dp), dimension(4,4), intent(in) :: M real(kind=dp), dimension(4), intent(inout) :: S - real :: sinw, cosw, omega, theta, costhet, xnyp, stok_I0, norme, frac_m1 + real :: sinw, cosw, omega, theta, costhet, xnyp real(kind=dp) :: v1pi, v1pj, v1pk, S1_0 real(kind=dp), dimension(4,4) :: ROP, RPO real(kind=dp), dimension(4) :: C, D diff --git a/src/stars.f90 b/src/stars.f90 index 61f3564ea..50142365c 100644 --- a/src/stars.f90 +++ b/src/stars.f90 @@ -265,7 +265,7 @@ subroutine repartition_energie_etoiles() real(kind=dp) :: wl_spectre_max, wl_spectre_min, wl_spectre_avg, wl_deviation - integer :: status, readwrite, blocksize,nfound,group,firstpix,nbuffer,npixels, naxes1_ref + integer :: status, readwrite, blocksize,nfound,group,firstpix,nbuffer,npixels real :: nullval integer, dimension(2) :: naxes logical :: anynull @@ -955,8 +955,6 @@ function is_inshock(id, iray, i_star, icell0, x, y, z, Thp, Tshock, Facc) real(kind=dp), intent(out) :: Thp, Tshock, Facc real(kind=dp) :: x, y, z, rho real(kind=dp) :: Tloc, vaccr, vmod2, rr, sign_z - real :: alpha_1 = 0.75 - real :: alpha_2 = 0.25 is_inshock = .false. if (.not.laccretion_shock) return @@ -995,7 +993,7 @@ function is_inshock(id, iray, i_star, icell0, x, y, z, Thp, Tshock, Facc) if (vaccr < 0.0_dp) then !Facc = 1/2 rho vs^3 Facc = 0.5 * (1d-3 * masseH * rho) * abs(vaccr)**3 - Tloc = ( alpha_1 * Facc / sigma )**0.25 + Tloc = ( 0.75 * Facc / sigma )**0.25 ! is_inshock = (Tloc > 0.5 * etoile(i_star)%T) is_inshock = (T_hp > 1.0_dp * etoile(i_star)%T) Thp = T_hp @@ -1148,7 +1146,7 @@ subroutine compute_stellar_parameters() character(len=1) :: s_age character(len=2) :: SpT - real :: L, R, T, M, minM, maxM, logg, minM_Allard, maxM_Allard, minM_Siess, maxM_Siess + real :: L, R, T, M, maxM, logg, minM_Allard, maxM_Allard, minM_Siess, maxM_Siess real(kind=dp) :: Gcm_to_Rsun integer :: age, k diff --git a/src/voigts.f90 b/src/voigts.f90 index 610ce91a3..503ffbfed 100644 --- a/src/voigts.f90 +++ b/src/voigts.f90 @@ -5,7 +5,7 @@ module voigts use constantes, only : pi, sqrtpi use mcfost_env, only : dp - + implicit none real(kind=dp), dimension(7) :: Aj, Bj, Cj @@ -60,21 +60,20 @@ function VoigtHumlicek(N, a, v, F) w4 = exp(u) - w4/v4 !cdexp, can be optimized by splitting the exp on real and imag parts endif VoigtHumlicek(i) = real(w4,kind=dp)!w4%re - if (present(F)) then + if (present(F)) then F(i) = aimag(w4)!w4%im endif enddo - + return end function VoigtHumlicek - + function VoigtThomson(N, a, v, vbroad) integer, intent(in) :: N real(kind=dp), intent(in) :: a, vbroad, v(N) real(kind=dp) :: VoigtThomson(N) - integer :: i real(kind=dp) :: f, rr, eta, ap ap = a * vbroad @@ -91,7 +90,7 @@ function VoigtThomson(N, a, v, vbroad) !because Thomson is in units of 1/vbroad/sqrtpi already VoigtThomson = eta * ( vbroad * f / pi / sqrtpi / ((vbroad*v)**2 + f**2) ) + & (1.0_dp - eta) * exp(-(v*vbroad/f)**2) * vbroad / f - + return end function VoigtThomson @@ -104,16 +103,16 @@ function max_voigt_profile(vth,a) ap = vth * a f = (vth**5 + 2.69269*vth**4 * ap + 2.42843*vth**3 * ap**2 + & 4.47163*vth**2*ap**3 + 0.07842*vth*ap**4 + ap**5)**(1./5.) - + eta = 1.36603*(ap/f) - 0.47719*(ap/f)**2+0.11116*(ap/f)**3 max_voigt_profile = eta / f / pi + (1.0-eta) / sqrt(pi) / f - return + return end function max_voigt_profile function dmax_voigt(vth,a,eps) !Evaluate the distance in units of vth, at which the voigt profile - ! reaches "eps * d" value : the value of the voigt profile in units of + ! reaches "eps * d" value : the value of the voigt profile in units of ! its peak. The Thomson approximation is used. real(kind=dp) :: dmax_voigt real, intent(in) :: eps @@ -126,9 +125,9 @@ function dmax_voigt(vth,a,eps) f = (vth**5 + 2.69269*vth**4 * ap + 2.42843*vth**3 * ap**2 + & 4.47163*vth**2*ap**3 + 0.07842*vth*ap**4 + ap**5)**(1./5.) - + eta = 1.36603*(ap/f) - 0.47719*(ap/f)**2+0.11116*(ap/f)**3 - + dmax_voigt = eta * sqrt(f/pi/peak_frac - f**2) + (1.0 - eta) * f * sqrt(-log(sqrt(pi)*f*peak_frac)) ! write(*,*) "eta=", eta ! write(*,*) "xp(L) = ", sqrt(f/pi/peak_frac - f**2) @@ -153,4 +152,4 @@ function dmax_lorentz(vth, a, eps) return end function dmax_lorentz -end module voigts \ No newline at end of file +end module voigts