-
Notifications
You must be signed in to change notification settings - Fork 10
/
sw.f90
executable file
·128 lines (102 loc) · 3.3 KB
/
sw.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
MODULE sw
IMPLICIT NONE
INTEGER,PARAMETER :: r8 = SELECTED_REAL_KIND(12) ! 8 byte real
REAL(r8), PARAMETER :: &
alpha = 0. ,&! angle of axis rotation about the poles
!
!
g = 9.80616d0 ,&! gravitational acceleration (m/s^2)
a = 6371229.d0,&! Earth's radius in m
omega = 7.29212d-5,&! angular velocity 1/s
pi = 3.14159265358979323846_R8,& ! pi
deg2rad = pi/180.d0
CONTAINS
!
!**************************************************************************
!
! Surface geopotential
!
!**************************************************************************
!
REAL(r8) FUNCTION surface_geopotential(lon,lat,test_case)
IMPLICIT NONE
REAL(r8), INTENT(IN) :: lon, lat
INTEGER , INTENT(IN) :: test_case
REAL(r8) :: r, r0, p1(2), p2(2)
select case (test_case)
case(1)
surface_geopotential = 0.0
case(2)
surface_geopotential = 0.0
case(5)
r0 = pi/9.
p1(1) = pi/2.
p1(2) = pi/6.
p2(1) = lon
p2(2) = lat
r = SQRT(MIN(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2))))
surface_geopotential = 2000.0*g*(1.0-(r/r0))
end select
END FUNCTION surface_geopotential
!
!**************************************************************************
!
! heights
!
!**************************************************************************
!
REAL(r8) FUNCTION height(lon,lat,test_case)
IMPLICIT NONE
REAL(r8), INTENT(IN) :: lon, lat
INTEGER , INTENT(IN) :: test_case
REAL(r8) :: Ubar
REAL(r8) :: r, r0, p1(2), p2(2)
select case (test_case)
case(1)
r0 = a/3.
p1(1) = pi/2.
p1(2) = 0.0
p2(1) = lon
p2(2) = lat
r = a*acos( sin(p1(2))*sin(p2(2)) + cos(p1(2))*cos(p2(2))*cos(p2(1)-p1(1)) )
if (r < r0) then
height = (1000.0/2.0)*(1.0+cos(pi*r/r0))
else
height = 0.0
endif
case(2)
Ubar = (2.0*pi*a)/(12.0*86400.0)
height = 2.94e4 - (a*omega*Ubar + (Ubar*Ubar)/2.) * &
( -1.*cos(lon)*cos(lat)*sin(alpha) + &
sin(lat)*cos(alpha) ) ** 2
case(5)
Ubar = (2.0*pi*a)/(12.0*86400.0)
height = 5400.*g - (a*omega*Ubar + (Ubar*Ubar)/2.) * &
( -1.*cos(lon)*cos(lat)*sin(alpha) + &
sin(lat)*cos(alpha) ) ** 2
end select
END FUNCTION height
!
!********************************************************************
!
! wind components
!
!********************************************************************
!
REAL(r8) FUNCTION u_wind(lon,lat,test_case)
IMPLICIT NONE
REAL(r8), INTENT(IN) :: lon,lat
INTEGER, INTENT(IN) :: test_case
REAL(r8) :: Ubar
Ubar = (2.0*pi*a)/(12.0*86400.0)
u_wind = Ubar*( cos(lat)*cos(alpha) + cos(lon)*sin(lat)*sin(alpha))
END FUNCTION u_wind
REAL(r8) FUNCTION v_wind(lon,lat,test_case)
IMPLICIT NONE
REAL(r8), INTENT(IN) :: lon,lat
INTEGER, INTENT(IN) :: test_case
REAL(r8) :: Ubar
Ubar = (2.0*pi*a)/(12.0*86400.0)
v_wind = -Ubar*sin(lon)*sin(alpha)
END FUNCTION v_wind
END MODULE sw