the params require to be read in phys template file phys_init
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(rc_fluid), | intent(inout) | :: | fl |
subroutine radiative_cooling_init(fl) use mod_global_parameters type(rc_fluid), intent(inout) :: fl double precision, dimension(:), allocatable :: t_table double precision, dimension(:), allocatable :: L_table double precision :: ratt, fact1, fact2, fact3, dL1, dL2 double precision :: tstep, Lstep integer :: ntable, i, j logical :: jump call rc_params_read(fl) ! Init for interpolatable tables allocate(fl%tcool(1:fl%ncool), fl%Lcool(1:fl%ncool)) allocate(fl%Yc(1:fl%ncool)) fl%tcool(1:fl%ncool) = zero fl%Lcool(1:fl%ncool) = zero ! Read in the selected cooling curve select case(fl%coolcurve) case('JCcorona') if(mype ==0) & print *,'Use Colgan & Feldman (2008) cooling curve' if(mype ==0) & print *,'This version only till 10000 K, beware for floor T treatment' ntable = n_JCcorona allocate(t_table(1:ntable)) allocate(L_table(1:ntable)) t_table(1:ntable) = t_JCcorona(1:n_JCcorona) L_table(1:ntable) = l_JCcorona(1:n_JCcorona) case('Colgan_DM') if(mype==0)& print *, 'Combination of Colgan (2008) for high temperatures and' if(mype==0)& print *, 'Dalgarno & McCray (1972), DM2, for low temperatures' ntable = n_Colgan + n_DM_2 allocate(t_table(1:ntable)) allocate(L_table(1:ntable)) t_table(1:n_DM_2) = t_DM_2(1:n_DM_2) L_table(1:n_DM_2) = L_DM_2(1:n_DM_2) t_table(n_DM_2+1:ntable) = t_Colgan(1:n_Colgan) L_table(n_DM_2+1:ntable) = l_Colgan(1:n_Colgan) case('TRML_narrow') if(mype==0)& print *, 'TRML narrow cooling' ntable = n_TRML_narrow allocate(t_table(1:ntable)) allocate(L_table(1:ntable)) t_table(1:ntable) = t_TRML_narrow(1:n_TRML_narrow) L_table(1:ntable) = l_TRML_narrow(1:n_TRML_narrow) case('TRML_wide') if(mype==0)& print *, 'TRML wide cooling' ntable = n_TRML_wide allocate(t_table(1:ntable)) allocate(L_table(1:ntable)) t_table(1:ntable) = t_TRML_wide(1:n_TRML_wide) L_table(1:ntable) = l_TRML_wide(1:n_TRML_wide) case default call mpistop("This coolingcurve is unknown") end select ! create cooling table(s) for use in amrvac fl%tcoolmax = t_table(ntable) fl%tcoolmin = t_table(1) ratt = (fl%tcoolmax-fl%tcoolmin)/(dble(fl%ncool-1) + smalldouble) fl%tcool(1) = fl%tcoolmin fl%Lcool(1) = L_table(1) fl%tcool(fl%ncool) = fl%tcoolmax fl%Lcool(fl%ncool) = L_table(ntable) do i=2,fl%ncool ! loop to create one table fl%tcool(i) = fl%tcool(i-1)+ratt do j=1,ntable-1 ! loop to create one spot on a table ! Second order polynomial interpolation, except at the outer edge, ! or in case of a large jump. if(fl%tcool(i) < t_table(j+1)) then if(j.eq. ntable-1 )then fact1 = (fl%tcool(i)-t_table(j+1)) & /(t_table(j)-t_table(j+1)) fact2 = (fl%tcool(i)-t_table(j)) & /(t_table(j+1)-t_table(j)) fl%Lcool(i) = L_table(j)*fact1 + L_table(j+1)*fact2 exit else dL1 = L_table(j+1)-L_table(j) dL2 = L_table(j+2)-L_table(j+1) jump =(max(dabs(dL1),dabs(dL2)) > 2*min(dabs(dL1),dabs(dL2))) end if if( jump ) then fact1 = (fl%tcool(i)-t_table(j+1)) & /(t_table(j)-t_table(j+1)) fact2 = (fl%tcool(i)-t_table(j)) & /(t_table(j+1)-t_table(j)) fl%Lcool(i) = L_table(j)*fact1 + L_table(j+1)*fact2 exit else fact1 = ((fl%tcool(i)-t_table(j+1)) & * (fl%tcool(i)-t_table(j+2))) & / ((t_table(j)-t_table(j+1)) & * (t_table(j)-t_table(j+2))) fact2 = ((fl%tcool(i)-t_table(j)) & * (fl%tcool(i)-t_table(j+2))) & / ((t_table(j+1)-t_table(j)) & * (t_table(j+1)-t_table(j+2))) fact3 = ((fl%tcool(i)-t_table(j)) & * (fl%tcool(i)-t_table(j+1))) & / ((t_table(j+2)-t_table(j)) & * (t_table(j+2)-t_table(j+1))) fl%Lcool(i) = L_table(j)*fact1 + L_table(j+1)*fact2 & + L_table(j+2)*fact3 exit end if end if end do ! end loop to find create one spot on a table end do ! end loop to create one table ! Go from logarithmic to actual values. fl%tcool(1:fl%ncool) = 10.0D0**fl%tcool(1:fl%ncool) fl%Lcool(1:fl%ncool) = 10.0D0**fl%Lcool(1:fl%ncool) ! Change unit of table if SI is used instead of cgs if (si_unit) fl%Lcool(1:fl%ncool) = fl%Lcool(1:fl%ncool) * 10.0d0**(-13) ! Scale both T and Lambda fl%tcool(1:fl%ncool) = fl%tcool(1:fl%ncool) / unit_temperature fl%Lcool(1:fl%ncool) = fl%Lcool(1:fl%ncool) * unit_numberdensity**2 * unit_time / unit_pressure * (1.d0+2.d0*He_abundance_rc) fl%tcoolmin = fl%tcool(1)+smalldouble ! avoid pointless interpolation ! smaller value for lowest temperatures from cooling table and user's choice if (fl%tlow==bigdouble) fl%tlow=fl%tcoolmin fl%tcoolmax = fl%tcool(fl%ncool) fl%lgtcoolmin = dlog10(fl%tcoolmin) fl%lgtcoolmax = dlog10(fl%tcoolmax) fl%lgstep = (fl%lgtcoolmax-fl%lgtcoolmin) * 1.d0 / (fl%ncool-1) deallocate(t_table) deallocate(L_table) if( fl%coolmethod == 'exact' ) then fl%tref = fl%tcoolmax fl%lref = fl%Lcool(fl%ncool) fl%Yc(fl%ncool) = zero do i=fl%ncool-1, 1, -1 fl%Yc(i) = fl%Yc(i+1) do j=1,100 tstep = 1.0d-2*(fl%tcool(i+1)-fl%tcool(i)) call findL(fl%tcool(i+1)-j*tstep, Lstep, fl) fl%Yc(i) = fl%Yc(i) + fl%lref/fl%tref*tstep/Lstep end do end do end if rc_gamma_1=rc_gamma-1.d0 invgam = 1.d0/rc_gamma_1 contains subroutine rc_params_read(fl) use mod_global_parameters, only: unitpar,par_files use mod_constants, only: bigdouble type(rc_fluid), intent(inout) :: fl integer :: n integer :: ncool = 4000 double precision :: cfrac=0.1d0 !> Name of cooling curve character(len=std_len) :: coolcurve='JCcorona' !> Name of cooling method character(len=std_len) :: coolmethod='exact' !> Fixed temperature not lower than tlow logical :: Tfix=.false. !> Lower limit of temperature double precision :: tlow=bigdouble logical :: rad_cut=.false. double precision :: rad_cut_hgt=0.5d0 double precision :: rad_cut_dey=0.15d0 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, Tfix, rad_cut, rad_cut_hgt, rad_cut_dey do n = 1, size(par_files) open(unitpar, file=trim(par_files(n)), status="old") read(unitpar, rc_list, end=111) 111 close(unitpar) end do fl%ncool=ncool fl%coolcurve=coolcurve fl%coolmethod=coolmethod if (coolmethod /= 'exact') call mpistop("Unknown cooling method: "//trim(coolmethod)) fl%tlow=tlow fl%Tfix=Tfix fl%cfrac=cfrac fl%rad_cut=rad_cut fl%rad_cut_hgt=rad_cut_hgt fl%rad_cut_dey=rad_cut_dey end subroutine rc_params_read end subroutine radiative_cooling_init