COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cstdatmos0-multi-seg.f
Go to the documentation of this file.
1 ! Programs contained in this file manage the atmosphere in
2 ! the vertical direction.
3 ! stdatmos2.d is used to describe the atmosphere. > 11km
4 !
5 ! cstdatmos0: to read data for the atmospeher
6 ! cvh2den : vertical height to denstiy
7 ! cvh2thick : vertical height to thickness of air
8 ! cvthick2h : vertical thicknes to height.
9 ! cvh2denp : d rho/ dz
10 ! cvh2den2p : d (d rho/ dz)/dz
11 !
12 
13 ! They follow the folloing formulas.
14 !
15 ! rho = rho0 * (1+ a(z-z0)/H(z0))**(-1-1/a) (a != 0)
16 ! = rho0 * exp(- (z-z0)/H) (a =0; hence H is const)
17 
18 !
19 ! The gramage between given heights, z1 and z2 is by
20 !
21 ! d = d0 *(fd(z1) - fd(z2)) where
22 !
23 ! fd(z) = (1+ a(z-z0)/H(z0))**(-1/a) (a != 0)
24 !
25 ! = exp(-(z-z0)/H ) (a = 0)
26 !
27 ! where d0 = rho0*H(z0)
28 !
29 ! Other formula is used at < 11 km. I.e.,
30 !
31 ! rho = rho00*( (ha - z)/hl )** pw
32 !
33 !
34 
35 
36 !
37  subroutine cstdatmos0
38  implicit none
39 #include "Zmanager.h"
40 #include "Zearth.h"
41 #include "Zstdatmos.h"
42 
43  integer icon, ios, i
44  external blkstdatmos
45 c---- include 'Zstdatmosf.h' ! rho=rho(h) function
46 #include "Zstdatmosf.h"
47 
48  if(first) then
49  call copenf(tempdev, "stdatmos2%d", icon)
50  if(icon .ne. 0) stop 9999
51  call cskipcomment(tempdev, icon)
52  if(icon .ne. 0) stop 9999
53 
54  nodes = 0
55  do while( .true. )
56  read(tempdev, *, iostat=ios)
57  * znode(nodes+1), tnode(nodes+1), pnode(nodes+1),
58  * rho(nodes+1), alfa(nodes+1), d0(nodes+1),
59  * dsum(nodes+1), scaleh(nodes+1)
60  if(ios .ne. 0) goto 10
61  if(nodes .ge. maxnodes) then
62  write(*,*) 'numbr of nodes for atmosphere > ', maxnodes
63  stop 9999
64  endif
65  nodes = nodes + 1
66 ! write(0,*) znode(nodes), alfa(nodes)
67  enddo
68  10 continue
69 ! write(*, *) " # of nodal points =", nodes
70 
71  do i=2, nodes
72  if(alfa(i-1) .ne. 0.) then
73  fd1i(i) = fd1(znode(i), alfa(i-1), znode(i-1),
74  * scaleh(i-1))
75  rhop(i) = rho(i-1) *(-1.d0 -1.d0/alfa(i-1)) *
76  * alfa(i-1)/scaleh(i-1)
77  pwp(i) =-2.d0-1.d0/alfa(i-1)
78  else
79  fd0i(i) = fd0(znode(i), znode(i-1), scaleh(i-1) )
80  endif
81  enddo
82  fd3 = ((ha-znode(2))/hl)**(pw+1.d0)
83  first = .false.
84  endif
85  end
86 ! ************************* to be updated
87  block data blkstdatmos
88 ! *************************
89  implicit none
90 c---- include 'Zearth.h'
91 #include "Zearth.h"
92 c---- include 'Zstdatmos.h'
93 #include "Zstdatmos.h"
94  data first/.true./
95  data ha/45.3425d3/, hl/11.865d3/,
96  * pw /4.36815d0/, rho00/3.50618d-3/
97  end
98 ! *********************************
99  real*8 function cvh2den(z)
100 ! *********************************
101 ! gives density of air at a given vertical height.
102 ! h: vertical height in m
103 ! function value: density in kg/m3.
104 !
105  implicit none
106 c---- include 'Zearth.h'
107 #include "Zearth.h"
108 c---- include 'Zstdatmos.h'
109 #include "Zstdatmos.h"
110 !
111  integer i
112 
113  real*8 z
114 
115  if(first) call cstdatmos0
116 
117 ! check if z < 11.2km
118  if(z .lt. znode(2)) then
119  cvh2den = rho00 * ( (ha - z)/hl )**pw
120  else
121  do i = 3, nodes
122  if(z .lt. znode(i) .or. i .eq. nodes ) then
123  if(alfa(i-1) .ne. 0.)then
124  if(z .lt. znode(i-1) ) then
125  write(0,*) ' z=',z, ' znode =', znode(i-1)
126  endif
127  cvh2den =
128  * rho(i-1) *
129  * (1.d0 + alfa(i-1)*
130  * (z-znode(i-1))/scaleh(i-1) )**(-1.0-1.d0/alfa(i-1))
131  goto 10
132  else
133  cvh2den =
134  * rho(i-1) * exp(- (z-znode(i-1))/scaleh(i-1))
135  goto 10
136  endif
137  endif
138  enddo
139  10 continue
140  endif
141  end
142 ! *************************
143  real*8 function cvh2thick(z)
144 ! vertical height to thickness of air above that point.
145 ! z: input. vertical height in m
146 ! function value. output. thickness in kg/m2.
147 !
148 !
149  implicit none
150 c---- include 'Zearth.h'
151 #include "Zearth.h"
152 c---- include 'Zstdatmos.h'
153 #include "Zstdatmos.h"
154 !
155  integer i
156 
157  real*8 z
158 
159 c---- include 'Zstdatmosf.h'
160 #include "Zstdatmosf.h"
161 
162 
163  if(first) call cstdatmos0
164 
165 ! check if z < 11.2km
166  if(z .lt. znode(2)) then
167  cvh2thick =hl* rho00 /(pw+1.d0)* (
168  * ( (ha-z)/hl)**(pw+1.d0) -
169  * fd3 )
170  * + dsum(2)
171 ! where fd3 is
172 ! * ((ha-znode(2))/hl)**(pw+1.d0) )
173  else
174 ! The gramage between given heights, z1 and z2 is by
175 !
176 ! d = d0 *(fd(z1) - fd(z2)) where
177 !
178 ! fd(z) = (1+ a(z-z0)/H(z0))**(-1/a) (a != 0)
179 !
180 ! = exp(-(z-z0)/H ) (a = 0)
181  do i = 3, nodes-1
182  if(z .lt. znode(i) .or. i .eq. nodes-1 ) then
183  if(alfa(i-1) .ne. 0.)then
184  cvh2thick = d0(i-1) * (
185  * fd1( z, alfa(i-1), znode(i-1), scaleh(i-1) )
186  * - fd1i(i)
187  * ) + dsum(i)
188 
189 ! where fd1i(i) is
190 ! * fd1(znode(i), alfa(i-1),
191 ! * znode(i-1), scaleh(i-1) )
192 
193  goto 10
194  else
195  cvh2thick = d0(i-1) * (
196  * fd0(z, znode(i-1), scaleh(i-1) )
197  * - fd0i(i)
198  * ) + dsum(i)
199 ! where fd0i is
200 ! * fd0(znode(i), znode(i-1), scaleh(i-1) )
201 ! * )
202  goto 10
203  endif
204  endif
205  enddo
206  10 continue
207  endif
208  end
209 ! *************************
210  real*8 function cvthick2h(t)
211 ! vertical thickness to vertical height conversion
212 ! t: input. vertical tikness in kg/m2
213 ! function value. output. height in m.
214 !
215 !
216  implicit none
217 c---- include 'Zearth.h'
218 #include "Zearth.h"
219 c---- include 'Zstdatmos.h'
220 #include "Zstdatmos.h"
221 !
222  integer i
223 
224  real*8 t, temp
225 
226 
227  if(first) call cstdatmos0
228 
229 ! check if z < 11.2km
230  if(t .gt. dsum(2)) then
231 ! t =hl* rho00 /(pw+1.d0)* (
232 ! * ( (ha-z)/hl)**(pw+1.d0) -
233 ! * ((ha-znode(2))/hl)**(pw+1.d0) )
234 ! * + dsum(2)
235 ! solve above eq.
236  temp =( t - dsum(2))/hl/rho00 * (pw+1.d0)
237  * + fd3
238  cvthick2h =ha - hl* temp**(1.d0/(pw+1.d0))
239  else
240 ! The gramage between given heights, z1 and z2 is by
241 !
242 ! d = d0 *(fd(z1) - fd(z2)) where
243 !
244 ! fd(z) = (1+ a(z-z0)/H(z0))**(-1/a) (a != 0)
245 !
246 ! = exp(-(z-z0)/H ) (a = 0)
247  do i = 3, nodes-1
248  if(t .gt. dsum(i) .or. i .eq. nodes-1 ) then
249  if(alfa(i-1) .ne. 0.)then
250 ! cvh2thick = d0(i-1) * (
251 ! * fd1( z, alfa(i-1), znode(i-1), scaleh(i-1) )
252 ! * - fd1(znode(i), alfa(i-1),
253 ! * znode(i-1), scaleh(i-1) )
254 ! * ) + dsum(i)
255 ! solve above eq.
256  temp =( t -dsum(i))/d0(i-1)
257  * + fd1i(i)
258  cvthick2h = (temp**(-alfa(i-1)) -1.d0)*
259  * scaleh(i-1)/alfa(i-1) + znode(i-1)
260  goto 10
261  else
262 ! cvh2thick = d0(i-1) * (
263 ! * fd0(z, znode(i-1), scaleh(i-1) )
264 ! * -fd0(znode(i), znode(i-1), scaleh(i-1) )
265 ! * ) + dsum(i)
266 ! solve above eq.
267  temp = (t -dsum(i)) /d0(i-1) +
268  * + fd0i(i)
269 ! temp = exp(-(z-z0)/H )
270  cvthick2h =znode(i-1)- log(temp)*scaleh(i-1)
271  goto 10
272  endif
273  endif
274  enddo
275  10 continue
276  endif
277  end
278 ! *********************************
279  real*8 function cvh2denp(z)
280 ! *********************************
281 ! gives derivative of the density of
282 ! air at a given vertical height.
283 ! h: vertical height in m
284 ! function value: drho/dz density in kg/m4
285 !
286  implicit none
287 c---- include 'Zearth.h'
288 #include "Zearth.h"
289 c---- include 'Zstdatmos.h'
290 #include "Zstdatmos.h"
291 !
292  integer i
293 
294  real*8 z
295 
296  if(first) call cstdatmos0
297 
298 ! check if z < 11.2km
299  if(z .lt. znode(2)) then
300  cvh2denp =- rho00 * pw/hl* ( (ha - z)/hl )**(pw-1.d0)
301  else
302  do i = 3, nodes
303  if(z .lt. znode(i) .or. i .eq. nodes ) then
304  if(alfa(i-1) .ne. 0.)then
305  cvh2denp =
306 ! * rho(i-1) *(-1.d0 -1.d0/alfa(i-1)) *
307 ! * alfa(i-1)/scaleh(i-1) * = rhop(i)
308  * rhop(i) *
309  * (1.d0 + alfa(i-1)*
310 ! * (z-znode(i-1))/scaleh(i-1) )**(-2.d0-1.d0/alfa(i-1))
311  * (z-znode(i-1))/scaleh(i-1) )**pwp(i)
312  goto 10
313  else
314  cvh2denp =
315  * -rho(i-1) * exp(- (z-znode(i-1))/scaleh(i-1))
316  * /scaleh(i-1)
317  goto 10
318  endif
319  endif
320  enddo
321  10 continue
322  endif
323  end
324 ! *********************************
325  real*8 function cvh2den2p(z)
326 ! *********************************
327 ! gives double derivative of the density of
328 ! air at a given vertical height.
329 ! h: vertical height in m
330 ! function value: d (drho/dz)/dz density in kg/m5
331 !
332  implicit none
333 c---- include 'Zearth.h'
334 #include "Zearth.h"
335 c---- include 'Zstdatmos.h'
336 #include "Zstdatmos.h"
337 !
338  integer i
339 
340  real*8 z
341 
342  if(first) call cstdatmos0
343 
344 ! check if z < 11.2km
345  if(z .lt. znode(2)) then
346  cvh2den2p = rho00 * pw*(pw-1.d0)/hl/hl
347  * * ( (ha - z)/hl )**(pw-2.d0)
348  else
349  do i = 3, nodes
350  if(z .lt. znode(i) .or. i .eq. nodes ) then
351  if(alfa(i-1) .ne. 0.)then
352  cvh2den2p =
353  * rhop(i) * pwp(i) *alfa(i-1)/scaleh(i-1)*
354  * (1.d0 + alfa(i-1)*
355  * (z-znode(i-1))/scaleh(i-1) )**(pwp(i)-1.d0)
356  goto 10
357  else
358  cvh2den2p =
359  * rho(i-1) * exp(- (z-znode(i-1))/scaleh(i-1))
360  * /scaleh(i-1)**2
361  goto 10
362  endif
363  endif
364  enddo
365  10 continue
366  endif
367  end
368 
cvthick2h
real *8 function cvthick2h(t)
Definition: ciniSegAtoms.f:148
cvh2denp
real *8 function cvh2denp(z)
Definition: ciniSegAtoms.f:201
d3
block data cblkIncident data *Za1ry *HeightOfInj d3
Definition: cblkIncident.h:7
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
hl
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hl
Definition: Zstdatmos.h:7
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
copenf
subroutine copenf(io, fnin, icon)
Definition: copenf.f:8
cvh2den2p
real *8 function cvh2den2p(z)
Definition: ciniSegAtoms.f:251
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
ha
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hlhmi common comstdatm ha
Definition: Zstdatmos.h:7
cvh2den
real *8 function cvh2den(z)
Definition: ciniSegAtoms.f:54
false
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
cstdatmos0
subroutine cstdatmos0
Definition: cstdatmos0-multi-seg.f:38
fd1
! to be included just before the execution code ! density as a function of height real * fd1
Definition: Zstdatmosf.h:3
cskipcomment
subroutine cskipcomment(io, icon)
Definition: cskipComment.f:19
cvh2thick
real *8 function cvh2thick(z)
Definition: ciniSegAtoms.f:93