DISPATCH
cooling_mod.f90
1 !==============================================================================
2 !==============================================================================
3 MODULE cooling_mod
4  USE io_mod
5  USE trace_mod
6  USE gpatch_mod
7  USE scaling_mod
8  implicit none
9  integer, parameter:: dp=4
10  real, parameter:: t_mc=10.0
11  private
12  type, public:: cooling_t
13  real:: dtime
14  procedure(cooling_casiana), pointer:: function => null()
15  real, dimension(:,:,:), pointer:: tt => null()
16  real, dimension(:,:,:), pointer:: qq => null()
17  real, dimension(:,:,:), pointer:: e => null()
18  contains
19  procedure:: init
20  procedure:: pre_update
21  procedure:: courant_condition
22  end type
23  real:: courant=0.5
24  logical:: on=.false.
25  integer:: type=0, verbose=0
26  logical, save:: first_time=.true.
27  type(cooling_t), public:: cooling
28 CONTAINS
29 
30 !==============================================================================
31 !==============================================================================
32 SUBROUTINE init (self)
33  class(cooling_t):: self
34  integer:: iostat
35  namelist /cooling_params/ on, verbose, type, courant
36  !------------------------------------------------------------------------------
37  call trace%begin ('simple/cooling_t%init')
38  !$omp critical (input_cr)
39  if (first_time) then
40  first_time = .false.
41  rewind(io%input)
42  read (io%input, cooling_params, iostat=iostat)
43  write (io%output, cooling_params)
44  end if
45  !$omp end critical (input_cr)
46  select case (type)
47  case (1)
48  self%function => cooling_casiana
49  case (2)
50  self%function => cooling_dalgarno
51  case (3)
52  self%function => cooling_rosen_etal
53  case (4)
54  self%function => cooling_shure
55  case default
56  self%function => cooling_gnedin
57  end select
58  call trace%end
59 END SUBROUTINE init
60 
61 !==============================================================================
62 !> Implement actual cooling
63 !==============================================================================
64 SUBROUTINE pre_update (self, patch)
65  class(cooling_t):: self
66  class(gpatch_t), pointer:: patch
67  integer:: i1, i2, i3, l(3), u(3)
68  real:: scaling_constant, q1, q2, eth, tt, qq, d, e, px, py, pz
69  real:: t_min, t_max, q_min, q_max, e_min, e_max
70  !-----------------------------------------------------------------------------
71  call trace%begin ('simple/cooling_t%pre_update')
72  if (on) then
73  if (.not.allocated(patch%heating_per_unit_volume)) then
74  u = patch%gn
75  allocate (patch%heating_per_unit_volume(u(1),u(2),u(3)))
76  patch%heating_per_unit_volume = 0.0
77  end if
78  self%dtime=huge(1.)
79  scaling_constant = scaling%t/scaling%p
80  t_min = huge(1.)
81  t_max = 0.0
82  q_min = huge(1.)
83  q_max = 0.0
84  e_min = huge(1.)
85  e_max = 0.0
86  l = patch%mesh%lb
87  u = patch%mesh%ub
88  do i3=l(3),u(3)
89  do i2=l(2),u(2)
90  do i1=l(1),u(1)
91  !-------------------------------------------------------------------------
92  ! Subtract kinetic energy and convert to temperature
93  !-------------------------------------------------------------------------
94  d = patch%mem(i1,i2,i3,patch%idx%d ,patch%it,1)
95  e = patch%mem(i1,i2,i3,patch%idx%e ,patch%it,1)
96  px = patch%mem(i1,i2,i3,patch%idx%px,patch%it,1)
97  py = patch%mem(i1,i2,i3,patch%idx%py,patch%it,1)
98  pz = patch%mem(i1,i2,i3,patch%idx%pz,patch%it,1)
99  eth = e - 0.5*(px**2+py**2+pz**2)/d
100  tt = real(scaling%temp*(patch%gamma-1d0))*eth/d
101  !-------------------------------------------------------------------------
102  ! Cool above T_MC, heat below
103  !-------------------------------------------------------------------------
104  q1 = self%function(tt)*d*scaling_constant
105  q2 = (t_mc/tt-1.0)*eth/patch%dtime*courant
106  qq = merge(q1,q2,tt>t_mc)
107  !patch%heating_per_unit_volume(i1,i2,i3) = &
108  !patch%heating_per_unit_volume(i1,i2,i3) + qq
109  self%dtime = min(self%dtime, courant*eth/qq)
110  t_min = min(t_min,tt)
111  t_max = max(t_max,tt)
112  q_min = min(q_min,qq)
113  q_max = max(q_max,qq)
114  e_min = min(e_min,eth)
115  e_max = max(e_max,eth)
116  end do
117  end do
118  end do
119  if (patch%istep < 2) then
120  patch%heating_per_unit_volume = 0.0
121  end if
122  if (verbose > 0) then
123  write(stdout,*) 'simple/cooling_t%pre_update: T_min,max =', t_min, t_max
124  write(stdout,*) 'simple/cooling_t%pre_update: Q_min,max =', q_min, q_max
125  end if
126  end if
127  call trace%end
128 END SUBROUTINE pre_update
129 
130 !==============================================================================
131 !> Implement actual cooling
132 !==============================================================================
133 SUBROUTINE courant_condition (self, patch)
134  class(cooling_t):: self
135  class(gpatch_t), pointer:: patch
136  integer:: i1, i2, i3, l(3), u(3)
137  real, dimension(:,:,:), pointer:: e
138  !----------------------------------------------------------------------------
139  call trace%begin ('simple/cooling_t%courant_condition')
140  if (on .and. self%dtime < patch%dtime) then
141  if (verbose > 0) then
142  write (stdout,'(a,i6,1p,2e12.3)') &
143  'simple/cooling_t%courant_condition: WARNING', &
144  patch%id, self%dtime, patch%dtime
145  end if
146  !!!$omp atomic write
147  !!patch%dtime = self%dtime
148  end if
149  call trace%end
150 END SUBROUTINE courant_condition
151 
152 !==============================================================================
153 !==============================================================================
154 include "cooling_Casiana.f90"
155 include "cooling_Dalgarno.f90"
156 include "cooling_Gnedin.f90"
157 include "cooling_Rosen_etal.f90"
158 include "cooling_Shure.f90"
159 
160 END MODULE cooling_mod
The gpath_t layer now essentially only handles restarts.
Definition: gpatch_mod.f90:4
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Definition: io_mod.f90:4