source: trunk/MagicSoft/Simulation/Corsika/Mmcs/rhof.f

Last change on this file was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 4.2 KB
Line 
1C=======================================================================
2
3 DOUBLE PRECISION FUNCTION RHOF( ARG )
4
5C-----------------------------------------------------------------------
6C RHO (DENSITY) F(UNCTION)
7C
8C CALCULATES DENSITY (G/CM**3) OF ATMOSPHERE DEPENDING ON HEIGHT (CM)
9C (US STANDARD ATMOSPHERE)
10C THIS FUNCTION IS CALLED FROM ININKG, UPDATE, CERENE, CERENH
11C ARGUMENT:
12C ARG = HEIGHT IN CM
13C-----------------------------------------------------------------------
14
15 IMPLICIT NONE
16
17*KEEP,ATMOS.
18 COMMON /ATMOS/ AATM,BATM,CATM,DATM
19 DOUBLE PRECISION AATM(5),BATM(5),CATM(5),DATM(5)
20*KEEP,RUNPAR.
21 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
22 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
23 * MONIOU,MDEBUG,NUCNUC,
24 * CETAPE,
25 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
26 * N1STTR,MDBASE,
27 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
28 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
29 * ,GHEISH,GHESIG
30 COMMON /RUNPAC/ DSN,HOST,USER
31 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
32 REAL STEPFC
33 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
34 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
35 * N1STTR,MDBASE
36 INTEGER CETAPE
37 CHARACTER*79 DSN
38 CHARACTER*20 HOST,USER
39
40 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
41 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
42 * ,GHEISH,GHESIG
43c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
44c Try
45c------------------------------------------------------------
46*KEEP,PARPAR.
47 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
48 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
49 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
50 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
51 INTEGER ITYPE,LEVL
52c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
53*KEND.
54
55C*******************************************************************
56C Modificado por Aitor (5-febrero-98)
57
58 common /aitor/ aitoth
59 double precision aitoth
60C*******************************************************************
61
62
63 DOUBLE PRECISION ARG,H,RT
64 PARAMETER (RT=6348.0D5)
65
66C-----------------------------------------------------------------------
67
68CC IF ( DEBUG ) WRITE(MDEBUG,*) 'RHOF : ARG=',SNGL(ARG)
69
70c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
71c>> Modification (HZA trick) cancelled >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
72c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
73c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
74 IF ( ARG .LT. 4.D5 ) THEN
75 RHOF = BATM(1) * DATM(1) * EXP ( -ARG * DATM(1) )
76 ELSEIF ( ARG .LT. 1.D6 ) THEN
77 RHOF = BATM(2) * DATM(2) * EXP ( -ARG * DATM(2) )
78 ELSEIF ( ARG .LT. 4.D6 ) THEN
79 RHOF = BATM(3) * DATM(3) * EXP ( -ARG * DATM(3) )
80 ELSEIF ( ARG .LT. 1.D7 ) THEN
81 RHOF = BATM(4) * DATM(4) * EXP ( -ARG * DATM(4) )
82 ELSE
83 RHOF = CATM(5)
84 ENDIF
85c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
86cC*******************************************************************
87cC Modificado por Aitor (5-febrero-98)
88c
89c H = -RT + SQRT(RT**2 + (ARG/COS(aitoth))**2 +(2.0D0*RT*ARG))
90cC*******************************************************************
91c
92cC R = SQRT(CURPAR(7)**2+CURPAR(8)**2)
93cC H = SQRT((RT+ARG)**2+R**2)-RT
94cc print *,'RHOF>>',arg,r,h,curpar(7),curpar(8)
95c
96c IF ( H .LT. 4.D5 ) THEN
97c RHOF = BATM(1) * DATM(1) * EXP ( -H * DATM(1) )
98c ELSEIF ( H .LT. 1.D6 ) THEN
99c RHOF = BATM(2) * DATM(2) * EXP ( -H * DATM(2) )
100c ELSEIF ( H .LT. 4.D6 ) THEN
101c RHOF = BATM(3) * DATM(3) * EXP ( -H * DATM(3) )
102c ELSEIF ( H .LT. 1.D7 ) THEN
103c RHOF = BATM(4) * DATM(4) * EXP ( -H * DATM(4) )
104c ELSE
105c RHOF = CATM(5)
106c ENDIF
107c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
108
109 RETURN
110 END
Note: See TracBrowser for help on using the repository browser.