COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cgeomag.f
Go to the documentation of this file.
1 
4 !c test crdGeomag
5 ! real*8 year
6 ! read(*,*) year
7 ! call crdGeomag('../../Data/Geomag/wmm', year)
8 ! end
9 ! **************************************************************
10 ! *
11 ! * cgeomag: geomagnetic filed strength is obtained
12 ! *
13 ! **************************************************************
14 ! /usage/ call cgeomag(year, llh, h, icon)
15 ! year: real*8. input. such as 1990.5
16 ! llh: /coord/ input. position around the earth.
17 ! in 'llh' form is better. if not 'llh'
18 ! conversion is done here.
19 ! h: /magfield/. output. magnetic field is set in
20 ! the form of 'ned' (north,
21 ! east-down). The unit is T.
22 ! icon: output. integer*4 0---> o.k
23 ! 1---> input parameter wrong.
24 !
25 
26 !
32  subroutine cgeomag(yearin, llh, h, icon)
33  implicit none
34 #include "Zglobalc.h"
35 #include "Zcoord.h"
36 #include "Zmagfield.h"
37 #include "Zearth.h"
38  real*8 yearin
39  type(coord)::llh
40  type(magfield)::h
41  integer icon
42 
43 ! ************************************
44  integer nmax
45  parameter(nmax = 13)
46  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
47  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
48  integer nmx
49  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
50 ! ************************************
51 
52 !
53 !
54  real*4 r, sumn, sume, sumd, t, cost, sint, x, tlonr, gmnc,
55  * cosml, sinml, hmnc, temp
56  real*8 gn, ge, gd
57  real*4 ssumd, ssumn, ssume
58  integer m, n
59  type(coord)::cdata
60  real*4 kdpmnxsinn, kpmnxn, kpmnxisinn
61 
62 !
63 ! check data type
64  if(llh%sys .eq. 'llh') then
65  cdata = llh
66  else ! convet to llh
67  call ctranscoord2('llh', llh, cdata)
68  endif
69 !
70  icon = 0
71  if( abs(yearin-year) .gt. 7.) then
72  icon=1
73  call cerrormsg(
74  * ' Year spec. for geomag data is invalid.', 1)
75 ! elseif(abs(cdata.lat) .gt. 90.) then
76  elseif(abs(cdata%r(1)) .gt. 90.5d0) then
77  icon=2
78 ! elseif(abs(cdata.long) .gt. 360.) then
79  elseif(abs(cdata%r(2)) .gt. 360.5d0) then
80  icon=2
81 ! elseif(cdata.h .gt. 5000.d3) then
82  elseif(cdata%r(3) .gt. 100000.d3) then
83  icon=1
84 ! elseif(cdata.h .lt. -3000.d3) then
85  elseif(cdata%r(3) .lt. -3000.d3) then
86  icon=1
87  endif
88  if(icon .ne. 2) then
89 ! r=1./( 1.+cdata.h/Eradius )
90  r=1./( 1.+cdata%r(3)/eradius )
91  sumn=0.
92  sume=0.
93  sumd=0.
94 ! t=(90.-cdata.lat)*Torad
95  t=(90.-cdata%r(1))*torad
96  cost=cos(t)
97  sint=sin(t)
98  x=cost
99 ! tlonr=cdata.long*Torad
100  tlonr=cdata%r(2)*torad
101  do n=1, nmx
102  ssumn = 0.
103  ssume = 0.
104  ssumd = 0.
105  do m=0, n
106  gmnc= gnm(n, m)
107  cosml=cos(m*tlonr)
108  sinml=sin(m*tlonr)
109  if(m .gt. 0) then
110  hmnc=hnm(n,m)
111  ssumn = ssumn+ (gmnc*cosml+hmnc*sinml)*
112  * kdpmnxsinn(m, n, sint, x)
113  ssume = ssume+ m *(-gmnc*sinml + hmnc*cosml)
114  * * kpmnxisinn(m, n, sint, x)
115  ssumd = ssumd + (gmnc*cosml+ hmnc*sinml)
116  * * kpmnxn(m, n, sint, x)
117  else
118  ssumn = ssumn + gmnc*kdpmnxsinn(m, n, sint, x)
119  ssumd = ssumd + gmnc*kpmnxn(m, n, sint, x)
120  endif
121  enddo
122  temp = r**(n+2)
123  sumn = sumn + ssumn*temp
124  sume = sume + ssume*temp
125  sumd = sumd + ssumd*temp *(n + 1)
126  enddo
127 ! original formula gives data in nT.
128 ! north component
129  gn = -sumn /1.e9 ! to T.
130 ! east component
131  ge = -sume /1.e9
132 ! down
133  gd = -sumd /1.e9
134  call csetmagfield('ned', gn, ge, gd, h)
135  endif
136  if(icon .eq. 2) then
137  call cerrormsg('Geometrical input data wrong', 0)
138  endif
139  end
140 
143  subroutine crdgeomag(filepath, yearin)
144  implicit none
145 ! #include "Zmanagerp.h"
146  character*(*) filepath
147  real*8 yearin
148  integer klena
149  character*65 msg
150 
151 
152 ! ************************************
153  integer nmax
154  parameter(nmax = 13)
155  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
156  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
157  integer nmx
158  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
159 ! ************************************
160 
161  character*128 path
162  character*1 NULL
163  integer devn, icon, leng, kgetenv
164  save devn
165  data devn /12/
166  null = char(0)
167 
168  if(filepath .eq. ' ') then
169 ! assume it is Cosmos/Data/Geomag/igrf
170 ! get $COSMOSTOP
171  leng = kgetenv("COSMOSTOP"//null, path)
172  if(leng .eq. 0) then
173  call cerrormsg(
174  * 'Env. variable "COSMOSTOP" cannot be found', 0)
175  endif
176  path = path(1:leng)//"/Data/Geomag/igrf"
177  else
178  path = filepath
179  endif
180  call copenf(devn, path(1:klena(path)), icon)
181  if(icon .eq. 1) then
182  call cerrormsg(path, 1)
183  call cerrormsg(
184  * 'above file cannot be opened for geomagnetic data', 0)
185  endif
186  call cgmgigrf(devn, path, yearin)
187  end
188 ! *************************
194  subroutine cgmgigrf(devn, filepath, yearin)
195  implicit none
196 ! assume at least two year data is included.
197  integer devn ! input file dev. no.
198  character*(*) filepath
199  real*8 yearin ! input. year of geomag.
200 
201  integer i
202  character*1024 msg
203 
204 
205 ! ************************************
206  integer nmax
207  parameter(nmax = 13)
208  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
209  real*4 yearmin, yearmax
210  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
211  integer nmx
212  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
213 ! ************************************
214 ! format from vesion Cosmos 7.55 is diff. from
215 ! old one; the data format for old one cannot be
216 ! found in the web site of igrf. The format found
217 ! in the web site is
218 ! g/h n m 1900.0 1905.0 ... 2005 2010 SV
219 ! g 1 0 -31543 -3146 ... -29554.63 -29496.5 11.4
220 ! g 1 1 -2298 -2298 ... -1669.05 -1585.9 16.7
221 ! h 1 1 5922 5909 ... 5077.99 4945.1 -28.8
222 ! g 2 0 -677 -728 ... -2337.24 -2396.6 -11.3
223 ! ....
224 ! g 13 13 0 0 ... -0.18 -0.3 0.0
225 ! h 13 13 0 0 .. -0.82 -0.8 0.0
226 !
227  integer klena, msglen
228 
229  integer nf ! # of fields in the data
230  integer n1, n2
231  integer:: n11, n12, n21, n22, n, m
232  integer:: targetf1, targetf2
233 
234  real(4):: year1, year2, sv, dy, yearx
235 
236  rewind devn
237 
238 
239 ! first find the appropriate year
240  year = yearin
241  msg = ' '
242 ! format changed from igrf12: skip if first col. is # or c/s
243  do
244  read(devn, '(a)' ) msg
245  if( msg(1:1) /= "#" .and. msg(1:3) /= "c/s") exit
246  enddo
247 ! chagne for igrf12 ended. igrf10 can be read also.
248 
249 
250 ! now msg should start with "g/h"
251  call kcountfields(msg, nf)
252  if(nf < 5) then
253  write(0,*) ' input igrf data is invalid'
254  write(0,*) ' the first line is ', msg
255  write(0,*) ' detected in cgmgIgrf '
256  stop 9999
257  endif
258  call kcgetafield(msg, 4, n1, n2)
259  if(n1 <= 0) then
260  write(0,*) ' strange in cgmgIgrf'
261  stop 9999
262  endif
263  read( msg(n1:n2), *) year1
264  yearmin = year1
265 
266  call kcgetafield(msg, nf-1, n1, n2)
267  read( msg(n1:n2), *) yearmax
268  targetf2 = 5
269  if(year < year1 ) then
270  targetf1 = 4
271  targetf2 = 5
272  dy = year - year1
273  else
274  do i = 5, nf-1
275  call kcgetafield(msg, i, n1, n2) ! 2nd year to last year
276  read(msg(n1:n2), *) year2
277  if( year < year2 .and. year >= year1 ) then
278  targetf1 = i - 1
279  targetf2 = i
280  dy = year - year1
281  exit
282  else
283  year1 = year2
284  if( i == nf-1 ) then
285  targetf1 = nf -1
286  targetf2 = nf -1
287  dy = year - year1
288  year2 = year1+5
289  endif
290  endif
291  enddo
292  endif
293  if( abs(dy) > 11. ) then
294  write(0,*) 'yearin to cgmgIgrf is out of range'
295  write(0,*) 'yearin=', yearin
296  stop 8888
297  endif
298 ! read coeff.
299  nmx = 0
300  do while (.true.)
301  read(devn, '(a)', end=1000) msg
302  call kcgetafield(msg, 2, n1, n2)
303  read( msg(n1:n2), *) n
304  call kcgetafield(msg, 3, n1, n2)
305  read( msg(n1:n2), *) m
306  if( n > nmax .or. m > nmax ) then
307  write(0,*) ' too large m or n; m=', m, ' n=',n
308  write(0,*) ' max should be ', nmax
309  write(0,*) ' cgmgIgrf '
310  stop 1111
311  endif
312  nmx = max(n,m,nmx)
313  call kcgetafield(msg, 1, n1, n2) ! g/h
314  call kcgetafield(msg, targetf1, n11, n12) ! coef. for year1
315  call kcgetafield(msg, targetf2, n21, n22) ! // 2
316 
317 
318  if( msg(n1:n2) == 'g' ) then
319  read( msg(n11:n12), *) gnm(n,m)
320  read( msg(n21:n22), *) gnmd(n,m)
321  if( year2 <= yearmax) then
322  gnmd(n,m) = ( gnmd(n,m) - gnm(n,m) )/ (year2-year1)
323  else
324  call kcgetafield(msg, nf, n1, n2)
325  read(msg(n1:n2), *) sv
326  gnmd(n,m) = sv
327  endif
328  elseif(msg(n1:n2) == 'h' ) then
329  read( msg(n11:n12), *) hnm(n,m)
330  read( msg(n21:n22), *) hnmd(n,m)
331  if( year2 <= yearmax ) then
332  hnmd(n,m) = ( hnmd(n,m) - hnm(n,m) )/ (year2-year1)
333  else
334  call kcgetafield(msg, nf, n1, n2)
335  read(msg(n1:n2), *) sv
336  hnmd(n,m) = sv
337  endif
338  else
339  write(0,*) ' igff data strange the line is'
340  write(0,*) msg
341  stop 12345
342  endif
343  enddo
344  1000 continue
345  do n = 1, nmx
346  do m=0, n
347  gnm(n, m) = gnm(n, m) + dy*gnmd(n, m)
348  if(m .gt. 0) then
349  hnm(n, m) = hnm(n,m) + dy*hnmd(n, m)
350  endif
351  enddo
352  enddo
353 
354 
355  write(0, *) 'Geomagnetic data has been read from', filepath
356  write(0, *) ' year=',yearin, ' # of expansion terms=', nmx
357  end
358 
359 
360 ! *************************
362  subroutine cprgeomag
363  implicit none
364 ! print current geomag coeff.
365 
366 
367  integer m, n
368 
369 
370 ! ************************************
371  integer nmax
372  parameter(nmax = 13)
373  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
374  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
375  integer nmx
376  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
377 ! ************************************
378 
379  write(*,*) 'year=',year, ' nmax=', nmx
380  write(*,'(" m n gnm hnm dgnm dhnm")')
381  do n = 1, nmx
382  do m =0, n
383  if(m .eq. 0) then
384  write(*,'(2i3,2f7.0,2f7.1)') m, n,
385  * gnm(n, m), 0., gnmd(n, m), 0.
386  else
387  write(*,'(2i3,2f7.0,2f7.1)') m, n,
388  * gnm(n, m), hnm(n,m), gnmd(n,m), hnmd(n,m)
389  endif
390  enddo
391  enddo
392  end
393 
394 
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
cprgeomag
subroutine cprgeomag
Definition: cgeomag.f:477
cos
! for length to thickness conversion or v v ! integer maxnodes real Hinf ! This is used when making table for dim simulation ! The slant length for vertical height to km is ! divided by LenStep m steps ! It can cover the slant length of about km for cos
Definition: Zatmos.h:8
crdgeomag
subroutine crdgeomag(filepath, yearin)
Definition: cgeomag.f:133
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
cgeomag
subroutine cgeomag(yearin, llh, h, icon)
Definition: cgeomag.f:25
ctranscoord2
subroutine ctranscoord2(sys, a, b)
Definition: ctransCoord2.f:2
cgmgigrf
subroutine cgmgigrf(devn, filepath, yearin)
Definition: cgeomag.f:202
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
csetmagfield
subroutine csetmagfield(sys, b1, b2, b3, b)
brief set Calculated magnetic field to /magfield/ b
Definition: csetMagField.f:6
coord
Definition: Zcoord.h:43