cool_exact Subroutine

public subroutine cool_exact(qdt, wCT, wCTprim, wnew, x, fl)

Arguments

Type IntentOptional Attributes Name
double precision, intent(in) :: qdt
double precision, intent(in) :: wCT(nw_phys)
double precision, intent(in) :: wCTprim(nw_phys)
double precision, intent(inout) :: wnew(nw_phys)
double precision, intent(in) :: x(1:ndim)
type(rc_fluid), intent(in) :: fl

Calls

proc~~cool_exact~~CallsGraph proc~cool_exact cool_exact proc~calc_l_extended calc_l_extended proc~cool_exact->proc~calc_l_extended proc~findt findT proc~cool_exact->proc~findt proc~findy findY proc~cool_exact->proc~findy proc~get_pthermal get_pthermal proc~cool_exact->proc~get_pthermal proc~get_rfactor get_Rfactor proc~cool_exact->proc~get_rfactor proc~get_rho get_rho proc~cool_exact->proc~get_rho

Called by

proc~~cool_exact~~CalledByGraph proc~cool_exact cool_exact proc~radiative_cooling_add_source radiative_cooling_add_source proc~radiative_cooling_add_source->proc~cool_exact

Source Code

    subroutine cool_exact(qdt,wCT,wCTprim,wnew,x,fl)
      
    !  Cooling routine using exact integration method from Townsend 2009
      double precision, intent(in)    :: qdt, wCT(nw_phys), wCTprim(nw_phys), x(1:ndim)
      double precision, intent(inout) :: wnew(nw_phys)
      type(rc_fluid), intent(in) :: fl
      double precision :: Y1, Y2
      double precision :: L1, Tlocal2, pnew
      double precision :: rho, Te, rhonew, Rfactor
      double precision :: emin, Lmax, fact
      double precision :: de, emax

      rho = get_rho(wCT,x)
      Rfactor = get_Rfactor()
      Te=wCTprim(iw_e)/(rho*Rfactor)
      pnew = get_pthermal(wCT,x)
      rhonew = get_rho(wnew,x)

      fact = fl%lref*qdt/fl%tref

      emin = rhonew*fl%tlow*Rfactor*invgam
      Lmax = max(zero,pnew*invgam-emin)/qdt
      emax = max(zero,pnew*invgam-emin)
      !  Determine explicit cooling
      !  If temperature is below floor level, no cooling.
      !  Stop wasting time and go to next gridpoint.
      !  If the temperature is higher than the maximum,
      !  assume Bremsstrahlung
      if( Te<=fl%tcoolmin ) then
        L1 = zero
      else if( Te>=fl%tcoolmax )then
        call calc_l_extended(Te, L1, fl)
        L1 = L1*rho**2

        L1 = min(L1,Lmax)
        if(slab_uniform .and. fl%rad_cut .and. x(ndim) .le. fl%rad_cut_hgt) then
          L1 = L1*exp(-(x(ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
        end if
        wnew(iw_e) = wnew(iw_e)-L1*qdt
      else
        call findY(Te,Y1,fl)
        Y2 = Y1 + fact*rho*rc_gamma_1
        call findT(Tlocal2,Y2,fl)
        if(Tlocal2<=fl%tcoolmin) then
          de = emax
        else
          de = (Te-Tlocal2)*rho*Rfactor*invgam
        end if
        de = min(de,emax)
        if(slab_uniform .and. fl%rad_cut .and. x(ndim) .le. fl%rad_cut_hgt) then
          de = de*exp(-(x(ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
        end if
        wnew(iw_e) = wnew(iw_e)-de
      end if
    end subroutine cool_exact