-
Notifications
You must be signed in to change notification settings - Fork 0
/
modules.f90
168 lines (121 loc) · 4.56 KB
/
modules.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
! Last change: MG 20 Mar 2003 2:44 pm
! *******************************************************
! ** MODULE randmod **
! *******************************************************
MODULE randmod
! Purpose:
! --------
! To generate a random number between 0 and 1.
!
! Description
! -----------
! This random number generator combines:
! (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period 2^32.
! (2) A 3-shift shift-register generator, period 2^32-1,
! (3) Two 16-bit multiply-with-carry generators, period 597273182964842497>2^59
! Overall period>2^123; Default seeds x,y,z,w.
!
! Record and revisions
! --------------------
! Date Programmer Description of change
! ==== ========== =====================
! 1993 George Marsaglia Original code in F77
! 2000 Mike Metcalf Converted code to F90
! 1/10/2003 Manuel Gimond Modified code to ouptut
! values between 0 and 1
!
IMPLICIT NONE
PRIVATE
PUBLIC :: rand
CONTAINS
FUNCTION rand ()
USE rand_global
REAL :: rand
! This random number generator combines:
! (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period 2^32.
! (2) A 3-shift shift-register generator, period 2^32-1,
! (3) Two 16-bit multiply-with-carry generators, period 597273182964842497>2^59
! Overall period>2^123; Default seeds x,y,z,w.
! Set your own seeds with statement i=kisset(ix,iy,iz,iw).
!
xseed = 69069 * xseed + 1327217885
yseed = m (m (m (yseed, 13), - 17), 5)
zseed = 18000 * IAND (zseed, 65535) + ISHFT (zseed, - 16)
wseed = 30903 * IAND (wseed, 65535) + ISHFT (wseed, - 16)
rand = ABS(xseed + yseed + ISHFT (zseed, 16) + wseed)/2147483648. !modified by MG
CONTAINS
FUNCTION m(k, n)
INTEGER :: m, k, n
m = IEOR (k, ISHFT (k, n) )
END FUNCTION m
END FUNCTION rand
END MODULE randmod
! *******************************************************
! ** MODULE anglesum **
! *******************************************************
MODULE anglesum_mod
! Purpose:
! --------
! To calculate the resulting angle between [0,PI] or [0,2PI[
! from the summation of two angles.
!
! Description
! -----------
! If the angles to be summed are polar, the output will result
! in an angle ranging between 0 (downward) and PI (upward)
! Else, the resulting output will range between 0 and 2PI.
! Input angles : radian
! Output angle : radian
!
!
! Record and revisions
! --------------------
! Date Programmer Description of change
! ==== ========== =====================
! 12/27/2002 Manuel Gimond Original code
IMPLICIT NONE
PRIVATE
PUBLIC :: anglesum
CONTAINS
FUNCTION anglesum(switch,orig,add)
REAL :: anglesum
INTEGER,INTENT(IN) :: switch ! 0=> poleward , 1=>horizontal
REAL,INTENT(IN) :: orig ! original angle(in radian)
REAL,INTENT(IN) :: add ! angle to add (in radian)
REAL, PARAMETER :: PI = 3.14159265
INTEGER :: tip
SELECT CASE (switch)
CASE (0) !Polar angle (range [0,PI])
anglesum = ACOS(COS(orig + add))
CASE DEFAULT ! (range [0,2PI[
IF ( ( ABS((2. * PI) - (orig+add)) < 0.0000001 ) .AND. ((orig+add) >= 0.)) THEN
anglesum = orig + add
tip=1
ELSE IF ( ABS( COS(orig + add) - 1) < 0.0000001 ) THEN
anglesum = 0.0
tip=2
ELSE IF ( SIN(orig + add) < 0) THEN
anglesum = 2 * PI - ACOS( COS(orig + add) )
tip=3
ELSE IF ( SIN(orig + add) > 0 ) THEN
anglesum = ACOS( COS(orig + add) )
tip=4
ELSE IF ( ABS( COS(orig + add) + 1) < 0.001 ) THEN
anglesum = PI
tip=5
ELSE
anglesum = 0.
tip=6
END IF
IF( (anglesum - 2 * PI) > 0.0) THEN
! If sum of angle in the azimuth direction is greater than 2PI,
! then let the user know. Angle should alwayd be less than 2PI.
WRITE(*,*) 'BAD SUMMATION IN anglesum ROUTINE ', 'orig = ',orig, 'new = ', anglesum, tip
anglesum = anglesum - ABS(anglesum - 2 * PI) !! angle must no be greater than 2PI
END IF
IF( ABS(anglesum - 2 * PI) < 0.0001) THEN
anglesum = 0.0
END IF
END SELECT
END FUNCTION anglesum
END MODULE anglesum_mod