-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathinit_module.f90
156 lines (123 loc) · 4.12 KB
/
init_module.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
module init_module
use kind_module, only: i4b, dp
use math_module, only: pi=>math_pi
implicit none
private
integer, parameter, private :: nc = 2
real(kind=dp), dimension(nc), private :: &
lonc = (/5.0_dp*pi/6.0_dp, 7.0_dp*pi/6.0_dp/), &
latc = (/0.0_dp, 0.0_dp/)
public :: init_ghill, init_ghill2, init_cbell2, init_scyli2, init_ccbel2
contains
! Ritche 1987
subroutine init_ghill(lon,lat,gphi)
use planet_module, only: a=>planet_radius
use math_module, only: deg2rad=>math_deg2rad
use sphere_module, only: orthodrome
implicit none
real(kind=dp), dimension(:), intent(in) :: lon, lat
real(kind=dp), dimension(:,:), intent(inout) :: gphi
real(kind=dp), parameter :: &
phimax = 100.0_dp, xi = 0.0_dp, yi = 0.0_dp, L = 2.5e6_dp ! Gaussian hill
integer(kind=i4b) :: i, j, nx, ny
real(kind=dp) :: r, loni, lati
nx = size(lon)
ny = size(lat)
loni = xi*deg2rad
lati = yi*deg2rad
do j=1, ny
do i=1, nx
r = a*orthodrome(lon(i),lat(j),loni,lati)
gphi(i,j) = phimax * exp(-(pi*r/L)**2)
end do
end do
end subroutine init_ghill
subroutine init_ghill2(lon,lat,gphi)
use sphere_module, only: lonlat2xyz
implicit none
real(kind=dp), dimension(:), intent(in) :: lon, lat
real(kind=dp), dimension(:,:), intent(inout) :: gphi
real(kind=dp), parameter :: hmax = 0.95_dp, b0 = 5.0_dp
integer(kind=i4b) :: i, j, k, nx, ny
! real(kind=dp) :: x0, y0, z0, x, y, z
nx = size(lon)
ny = size(lat)
gphi(:,:) = 0.0_dp
do k=1, nc
! call lonlat2xyz(lonc(k),latc(k),x0,y0,z0)
do j=1, ny
do i=1, nx
! call lonlat2xyz(lon(i),lat(j),x,y,z)
! gphi(i,j) = gphi(i,j) + hmax*exp(-b0*((x-x0)**2+(y-y0)**2+(z-z0)**2))
gphi(i,j) = gphi(i,j) + hmax*exp(-2*b0*(1.0_dp-cos(lat(j))*cos(lon(i)-lonc(k))))
end do
end do
end do
end subroutine init_ghill2
subroutine init_cbell2(lon,lat,gphi)
use sphere_module, only: orthodrome
implicit none
real(kind=dp), dimension(:), intent(in) :: lon, lat
real(kind=dp), dimension(:,:), intent(inout) :: gphi
real(kind=dp), parameter :: &
hmax = 1.0_dp, r0 = 0.5_dp, b = 0.1_dp, c = 0.9_dp
integer(kind=i4b) :: i, j, k, nx, ny
real(kind=dp) :: hch, pir0r, r
nx = size(lon)
ny = size(lat)
gphi(:,:) = b
pir0r = pi/r0
hch = c*0.5_dp*hmax
do k=1, nc
do j=1, ny
do i=1, nx
r = orthodrome(lon(i),lat(j),lonc(k),latc(k))
if (r<r0) then
gphi(i,j) = gphi(i,j) + hch*(1.0_dp+cos(pir0r*r))
end if
end do
end do
end do
end subroutine init_cbell2
subroutine init_scyli2(lon,lat,gphi)
use sphere_module, only: orthodrome
implicit none
real(kind=dp), dimension(:), intent(in) :: lon, lat
real(kind=dp), dimension(:,:), intent(inout) :: gphi
real(kind=dp), parameter :: &
hmax = 1.0_dp, r0 = 0.5_dp, b = 0.1_dp, c = 1.0_dp
integer(kind=i4b) :: i, j, nx, ny
real(kind=dp) :: &
r1, r2, dlon0, dlon1, dlon2, dlat0, dlat1, dlat2
nx = size(lon)
ny = size(lat)
gphi(:,:) = b
dlon0 = r0/6.0_dp
dlat0 = 5.0_dp/12.0_dp*r0
do j=1, ny
do i=1, nx
r1 = orthodrome(lon(i),lat(j),lonc(1),latc(1))
r2 = orthodrome(lon(i),lat(j),lonc(2),latc(2))
dlon1 = abs(lon(i)-lonc(1))
dlon2 = abs(lon(i)-lonc(2))
dlat1 = lat(j)-latc(1)
dlat2 = lat(j)-latc(2)
if (((r1<=r0).and.(dlon1>=dlon0)).or. &
((r2<=r0).and.(dlon2>=dlon0)).or. &
((r1<=r0).and.(dlon1<dlon0).and.dlat1<-dlat0).or. &
((r2<=r0).and.(dlon2<dlon0).and.dlat2>dlat0)) then
gphi(i,j) = c
end if
end do
end do
end subroutine init_scyli2
subroutine init_ccbel2(lon,lat,gphi)
use sphere_module, only: orthodrome
implicit none
real(kind=dp), dimension(:), intent(in) :: lon, lat
real(kind=dp), dimension(:,:), intent(inout) :: gphi
real(kind=dp), parameter :: a = -0.8, b = 0.9
call init_cbell2(lon,lat,gphi)
gphi(:,:) = a*gphi(:,:)**2 + b
end subroutine init_ccbel2
end module init_module