Added heat source

.gitignore

@@ -43,6 +43,7 @@ import_eqdsk
 jorek2_main
 jorek2_target2vtk
 jorek2_powers
+jorek2_radiative
 jorek2_import_perturbation
 jorek2_import_perturbation_new_flags
 jorek2_import_perturbation_new_flags

namelists/lithiumtest.nml

@@ -4,18 +4,18 @@
   ADASname   = '96_li'                      ! Prefix of ADAS files
 
   tstep    = 0.01d0                         ! Timestep
-  ftime    = 300.d0                         ! Final time
+  ftime    = 10000.d0                       ! Final time
   nout     = 100                            ! Output every
 
   temperature    = 0.04d0                   ! Initial electron temperature (everything in JOREK-units, central_density=1.d0)
   density_plasma = 1.d0                     ! Plasma ion density
   density_imp    = 1.d-1                    ! Impurity ion density
 
-  t_rad  = 0.d0                             ! Turn on radiative cooling term after
+  t_rad  = 0.d00                            ! Turn on radiative cooling term after
   t_equi = 1.d99                            ! Assume coronal equilibrium after
 
-  C_cooling(1)  = 0.d0                      ! Linear Cooling rate dT/dt
-  C_cooling(2)  = 1.d99                     ! Start cooling after
-  C_cooling(3)  = 2.d-5                     ! Stop cooling, if Te < 1 eV
+  C_cooling(1)  = 9.d-6                     ! Linear Cooling rate dT/dt
+  C_cooling(2)  = 1.d3                      ! Start cooling after ...
+  C_cooling(3)  = 1.d-2                     ! Stop cooling, if Te < ...
 
 /

radiative/mod_non_coronal.f90

@@ -23,6 +23,7 @@ type non_coronal
   real*8                         :: density_plasma    !< constant density of ions of the background plasma
   real*8                         :: density_imp       !< constant density of the imuprity ions
   real*8, dimension(3)           :: C_cooling         !< Coefficient for linear aritifical cooling C(1): dT/dt C(2): when starts cooling, C(3): Tmin
+  real*8                         :: t_rad, t_equi     !< when starts radiative cooling, when go into coronal equilbrium
   real*8, allocatable            :: fractions(:)      !< Fractional charge states of impurites
 
   ! Constants
@@ -86,7 +87,9 @@ cor%density_plasma = density_plasma
 cor%density_imp    = density_imp
 cor%fractions      = 0.0d0
 cor%fractions(1)   = 1.0d0
-
+cor%C_cooling      = C_cooling
+cor%t_rad          = t_rad
+cor%t_equi         = t_equi
 
 call cor%calculate_GCRs()
 
@@ -173,7 +176,8 @@ ddensity_dt        = 0.d0
 dfractions_dt      = 0.d0
 
 ! ---- non coronal
-do iz=1,cor%n_Z0
+if ( cor%time .lt. cor%t_equi ) then
+  do iz=1,cor%n_Z0
     if ( iz .gt. 1) then
       dfractions_dt(iz) = dfractions_dt(iz)                                                  &
                           + cor%fractions(iz-1) * cor%density * cor%ion_rate(iz-1)   & ! Gain by Ionization    from (iz-1) to iz
@@ -184,24 +188,26 @@ do iz=1,cor%n_Z0
                           + cor%fractions(iz+1) * cor%density * cor%rec_rate(iz)     & ! Gain by Recombination from (iz+1) to iz
                           - cor%fractions(iz)   * cor%density * cor%ion_rate(iz)       ! Loss by Ionization    from iz     to (iz+1)
      end if
-end do
+  end do
 
-cor%fractions     = fractions0   + dfractions_dt *dtime
-ddensity_dt       = cor%density_imp * dot_product(dfractions_dt,cor%Z)
+  cor%fractions     = fractions0   + dfractions_dt *dtime
+  ddensity_dt       = cor%density_imp * dot_product(dfractions_dt,cor%Z)
 
 
 ! --- coronal equilibrium
-!cor%fractions = cor%equilibrium()
-!ddensity_dt   = cor%density_imp * ( dot_product(cor%fractions,cor%Z) - dot_product(fractions0,cor%Z) ) / dtime
-
-
+else
+  cor%fractions = cor%equilibrium()
+  ddensity_dt   = cor%density_imp * ( dot_product(cor%fractions,cor%Z) - dot_product(fractions0,cor%Z) ) / dtime
+end if
 
 !Diffusion and fueling
 !dfractions_dt    = cor%fractions/tau
 !dfractions_dt(1) = sum(cor%dfractions_dt(2:cor%n_Z0))
 
-
-dtemperature_dt   = - 2.d0 / 3.d0 *   dnT_rad_dt / density_tot0 &
+if ( cor%time .gt. cor%t_rad) then
+  dtemperature_dt   = - 2.d0 / 3.d0 *   dnT_rad_dt / density_tot0
+end if
+  dtemperature_dt   = dtemperature_dt  &
                       - dT_cool_dt     &
                       - ddensity_dt * temperature0 / density_tot0
 
@@ -320,7 +326,11 @@ function calculate_cooling(cor)
 class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
 real*8                         :: calculate_cooling
 
-calculate_cooling = 0.d0
+if ( cor%time .gt. cor%C_cooling(2) .and. cor%temperature .gt. cor%C_cooling(3) ) then
+  calculate_cooling = cor%C_cooling(1)
+else
+  calculate_cooling = 0.d0
+end if
 
 end function calculate_cooling