-
Notifications
You must be signed in to change notification settings - Fork 2
/
setperturb.f90
182 lines (147 loc) · 4 KB
/
setperturb.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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
subroutine setperturb
! Random noise
use vars
use params
use microphysics, only: micro_field, index_water_vapor
use sgs, only: setperturb_sgs
implicit none
integer i,j,k,ptype,it,jt,n,ii,jj
real rrr,ranf_
real xxx,yyy,zzz,sss
call ranset_(3+nensemble)
ptype = perturb_type
call setperturb_sgs(ptype) ! set sgs fields
call task_rank_to_index(rank,it,jt)
select case (ptype)
case(-1) ! no perturbation is set
case(0)
sss = 0.
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
sss = sss + rrr
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(k.le.5) then
t(i,j,k)=t(i,j,k)+0.02*rrr*(6-k)
endif
end if
end do
end do
end do
if(masterproc) print*,'>>>>>> sss=',sss
case(1)
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(q0(k).gt.6.e-3) then
t(i,j,k)=t(i,j,k)+0.1*rrr
endif
end if
end do
end do
end do
case(2) ! warm bubble
if(masterproc) then
print*, 'initialize with warm bubble:'
print*, 'bubble_x0=',bubble_x0
print*, 'bubble_y0=',bubble_y0
print*, 'bubble_z0=',bubble_z0
print*, 'bubble_radius_hor=',bubble_radius_hor
print*, 'bubble_radius_ver=',bubble_radius_ver
print*, 'bubble_dtemp=',bubble_dtemp
print*, 'bubble_dq=',bubble_dq
end if
do k=1,nzm
zzz = z(k)
do j=1,ny
yyy = dy*(j+jt)
do i=1,nx
xxx = dx*(i+it)
if((xxx-bubble_x0)**2+YES3D*(yyy-bubble_y0)**2.lt.bubble_radius_hor**2 &
.and.(zzz-bubble_z0)**2.lt.bubble_radius_ver**2) then
rrr = cos(pi/2.*(xxx-bubble_x0)/bubble_radius_hor)**2 &
*cos(pi/2.*(yyy-bubble_y0)/bubble_radius_hor)**2 &
*cos(pi/2.*(zzz-bubble_z0)/bubble_radius_ver)**2
t(i,j,k) = t(i,j,k) + bubble_dtemp*rrr
micro_field(i,j,k,index_water_vapor) = &
micro_field(i,j,k,index_water_vapor) + bubble_dq*rrr
end if
end do
end do
end do
case(3) ! gcss wg1 smoke-cloud case
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(q0(k).gt.0.5e-3) then
t(i,j,k)=t(i,j,k)+0.1*rrr
endif
end if
end do
end do
end do
case(4) ! gcss wg1 arm case
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(z(k).le.200.) then
t(i,j,k)=t(i,j,k)+0.1*rrr*(1.-z(k)/200.)
endif
end if
end do
end do
end do
case(5) ! gcss wg1 BOMEX case
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(z(k).le.1600.) then
t(i,j,k)=t(i,j,k)+0.1*rrr
micro_field(i,j,k,index_water_vapor)= &
micro_field(i,j,k,index_water_vapor)+0.025e-3*rrr
endif
end if
end do
end do
end do
case(6) ! GCSS Lagragngian ASTEX
do k=1,nzm
do jj=1,ny_gl
j=jj-jt
do ii=1,nx_gl
i=ii-it
rrr=1.-2.*ranf_()
if(i.ge.1.and.i.le.nx.and.j.ge.1.and.j.le.ny) then
if(q0(k).gt.6.e-3) then
t(i,j,k)=t(i,j,k)+0.1*rrr
micro_field(i,j,k,index_water_vapor)= &
micro_field(i,j,k,index_water_vapor)+2.5e-5*rrr
endif
end if
end do
end do
end do
case default
if(masterproc) print*,'perturb_type is not defined in setperturb(). Exitting...'
call task_abort()
end select
end