DISPATCH
newton_mod.f90
1 !===============================================================================
2 !> Add Newton cooling to any solver
3 !===============================================================================
4 MODULE newton_mod
5  USE link_mod
6  USE gpatch_mod
7  USE mesh_mod
8  USE io_unit_mod
9  USE trace_mod
10  implicit none
11  private
12  type, public:: newton_t
13  class(gpatch_t), pointer:: patch
14  contains
15  procedure:: init
16  procedure:: pre_update
17  end type
18  logical:: on=.false.
19  integer:: verbose=0, axis=3
20  ! -- the default values are reasonable for the solar photosphere
21  real:: position=-0.5, scale=0.1, ee0=5.3, ee1=0.5, time=0.1
22 CONTAINS
23 
24 !===============================================================================
25 !> Initialize Newton cooling
26 !===============================================================================
27 SUBROUTINE init (self, link)
28  class(newton_t):: self
29  class(link_t), pointer :: link
30  ! ............................................................................
31  integer:: iostat
32  logical, save:: first_time=.true.
33  namelist /newton_params/ on, axis, position, scale, ee0, ee1, time
34  !-----------------------------------------------------------------------------
35  call trace%begin ('newton_t%init')
36  !$omp critical (input_cr)
37  if (first_time) then
38  first_time = .false.
39  rewind(io_unit%input)
40  read (io_unit%input, newton_params, iostat=iostat)
41  write (io_unit%output, newton_params)
42  end if
43  !$omp end critical (input_cr)
44  !-----------------------------------------------------------------------------
45  associate(task => link%task)
46  select type (task)
47  class is (gpatch_t)
48  self%patch => task
49  if (.not.allocated(task%heating_per_unit_volume)) &
50  allocate (task%heating_per_unit_volume(task%gn(1),task%gn(2),task%gn(3)))
51  end select
52  end associate
53  call trace%end()
54 END SUBROUTINE init
55 
56 !===============================================================================
57 !> Apply Newton cooling
58 !===============================================================================
59 SUBROUTINE pre_update (self)
60  class(newton_t) :: self
61  !.............................................................................
62  class(mesh_t), pointer:: m1, m2, m3
63  real, dimension(:,:,:), pointer:: d, e
64  real:: f, r, ee, s
65  integer:: ix, iy, iz
66  ! ----------------------------------------------------------------------------
67  d => self%patch%mem(:,:,:,self%patch%idx%d,self%patch%it,1)
68  e => self%patch%mem(:,:,:,self%patch%idx%e,self%patch%it,1)
69  if (on) then
70  m1 => self%patch%mesh(1)
71  m2 => self%patch%mesh(2)
72  m3 => self%patch%mesh(3)
73  if (axis==1) then
74  do iz = m3%li, m3%ui
75  do iy = m2%li, m2%ui
76  do ix = m1%li, m1%ui
77  s = m1%p + m1%r(ix)-position
78  f = exp(-s/scale)
79  r = f/((1.0+f)*time)
80  ee = ee0 + ee1*s
81  self%patch%heating_per_unit_volume(ix,iy,iz) = &
82  self%patch%heating_per_unit_volume(ix,iy,iz) + &
83  (ee-e(ix,iy,iz)/d(ix,iy,iz))*d(ix,iy,iz)*r
84  end do
85  end do
86  end do
87  else
88  do iz = m3%li, m3%ui
89  if (axis==3) then
90  s = m3%p + m3%r(iz)-position
91  f = exp(-s/scale)/time
92  r = f/((1.0+f)*time)
93  ee = ee0 + ee1*s
94  end if
95  do iy = m2%li, m2%ui
96  if (axis==2) then
97  s = m2%p + m2%r(iy)-position
98  f = exp(-s/scale)
99  r = f/((1.0+f)*time)
100  ee = ee0 + ee1*s
101  end if
102  do ix = m1%li, m1%ui
103  self%patch%heating_per_unit_volume(ix,iy,iz) = &
104  self%patch%heating_per_unit_volume(ix,iy,iz) + &
105  (ee-e(ix,iy,iz)/d(ix,iy,iz))*d(ix,iy,iz)*r
106  end do
107  end do
108  end do
109  end if
110  end if
111 END SUBROUTINE pre_update
112 
113 END MODULE newton_mod
The gpath_t layer now essentially only handles restarts.
Definition: gpatch_mod.f90:4
Add Newton cooling to any solver.
Definition: newton_mod.f90:4
Template module for mesh.
Definition: mesh_mod.f90:4