get_GOES_flux_grid Subroutine

    subroutine get_GOES_flux_grid(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,dV,xboxmin1,&
       xboxmin2,xboxmin3,xboxmax1,xboxmax2,xboxmax3,xbmin1,xbmin2,xbmin3,&
       xbmax1,xbmax2,xbmax3,fl,eflux_grid)
      use mod_global_parameters

      integer, intent(in)           :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
         ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
      double precision, intent(in)  :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
         ixImin3:ixImax3,1:ndim),dV(ixImin1:ixImax1,ixImin2:ixImax2,&
         ixImin3:ixImax3)
      double precision, intent(in)  :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
         ixImin3:ixImax3,nw)
      double precision, intent(in)  :: xboxmin1,xboxmin2,xboxmin3,xboxmax1,&
         xboxmax2,xboxmax3,xbmin1,xbmin2,xbmin3,xbmax1,xbmax2,xbmax3
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: eflux_grid

      integer :: ix1,ix2,ix3,ixO1,ixO2,ixO3,ixbmin1,ixbmin2,ixbmin3,ixbmax1,&
         ixbmax2,ixbmax3
      integer :: iE,numE,j,inbox
      double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
      double precision :: pth(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
         Te(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
         kbT(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3)
      double precision :: Ne(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
         EM(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3)
      double precision :: gff,fi,erg_SI

      ! check whether the grid is inside given box
      inbox=0
      if (xbmin1<xboxmax1 .and. xbmax1>xboxmin1) inbox=inbox+1
      if (xbmin2<xboxmax2 .and. xbmax2>xboxmin2) inbox=inbox+1
      if (xbmin3<xboxmax3 .and. xbmax3>xboxmin3) inbox=inbox+1

      if (inbox==ndim) then
        ! indexes for cells inside given box
        ixbmin1=ixOmin1;ixbmin2=ixOmin2;ixbmin3=ixOmin3;
        ixbmax1=ixOmax1;ixbmax2=ixOmax2;ixbmax3=ixOmax3;
        if (xbmax1>xboxmax1) ixbmax1=ixOmax1-ceiling((xbmax1-&
           xboxmax1)/dxlevel(1))
        if (xbmax2>xboxmax2) ixbmax2=ixOmax2-ceiling((xbmax2-&
           xboxmax2)/dxlevel(2))
        if (xbmax3>xboxmax3) ixbmax3=ixOmax3-ceiling((xbmax3-&
           xboxmax3)/dxlevel(3))
        if (xbmin1<xboxmin1) ixbmin1=ceiling((xboxmin1-xbmin1)/dxlevel(1))+&
           ixOmin1
        if (xbmin2<xboxmin2) ixbmin2=ceiling((xboxmin2-xbmin2)/dxlevel(2))+&
           ixOmin2
        if (xbmin3<xboxmin3) ixbmin3=ceiling((xboxmin3-xbmin3)/dxlevel(3))+&
           ixOmin3

        I0=1.07d-38 !photon flux index for observed at 1AU [photon cm3,cm3,cm3 m-2,m-2,m-2 s-1,s-1,s-1 keV-1],keV-1],keV-1]
        kb=const_kb
        keV=1.0d3*const_ev
        erg_SI=1.d-7
        A_cgs=1.d-8 ! Angstrom
        El=const_h*const_c/(8.d0*A_cgs)/keV ! 8 A
        Eu=const_h*const_c/(1.d0*A_cgs)/keV ! 1 A
        dE=0.1  ! keV
        numE=floor((Eu-El)/dE)
        call fl%get_pthermal(w,x,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
           ixImax3,ixbmin1,ixbmin2,ixbmin3,ixbmax1,ixbmax2,ixbmax3,pth)
        call fl%get_rho(w,x,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
           ixbmin1,ixbmin2,ixbmin3,ixbmax1,ixbmax2,ixbmax3,Ne)
        call fl%get_var_Rfactor(w,x,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
           ixImax3,ixbmin1,ixbmin2,ixbmin3,ixbmax1,ixbmax2,ixbmax3,Te)
        Te(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3)=pth(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3)/(Ne(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3)*Te(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3))*unit_temperature
        if (SI_unit) then
          Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)*unit_numberdensity/1.d6 !m-3,m-3,m-3 -> cm-3
          EM(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3)=(I0*(Ne(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3))**2)*dV(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3)*(unit_length*1.d2)**3 !cm-3,cm-3,cm-3
        else
          Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)*unit_numberdensity
          EM(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3)=(I0*(Ne(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3))**2)*dV(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
             ixbmin3:ixbmax3)*unit_length**3
        endif
        kbT(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3)=kb*Te(ixbmin1:ixbmax1,ixbmin2:ixbmax2,&
           ixbmin3:ixbmax3)/keV
        eflux_grid=0.0d0

        do iE=0,numE-1
          Ei=dE*iE+El
          do ix3=ixbmin3,ixbmax3
          do ix2=ixbmin2,ixbmax2
          do ix1=ixbmin1,ixbmax1
            if (kbT(ix1,ix2,ix3)>1.d-2*Ei) then
              if(kbT(ix1,ix2,ix3)<Ei) then
                gff=(kbT(ix1,ix2,ix3)/Ei)**0.4
              else
                gff=1.d0
              endif
              fi=(EM(ix1,ix2,ix3)*gff)*dexp(-Ei/(kbT(ix1,ix2,&
                 ix3)))/(Ei*dsqrt(kbT(ix1,ix2,ix3)))
              eflux_grid=eflux_grid+fi*dE*Ei
            endif
          enddo
          enddo
          enddo
        enddo
        eflux_grid=eflux_grid*keV*erg_SI
      endif

    end subroutine get_GOES_flux_grid