Line data Source code
1 : module filter_module
2 : use shr_kind_mod,only: r8 => shr_kind_r8
3 : use cam_logfile, only: iulog
4 : use cam_abortutils,only: endrun
5 : use edyn_geogrid,only: nlon, nlat
6 :
7 : implicit none
8 : private
9 : public :: ntrigs, ifax
10 : public :: filter_init
11 : public :: filter1, filter2, ringfilter
12 : public :: kut1, kut2
13 : !
14 : ! Coefficients and factors for fft. Sub setfft is called once per run from edyn_init
15 : !
16 : integer :: ntrigs ! = 3*nlon/2+1
17 : real(r8), allocatable :: trigs(:)
18 : integer :: ifax(13)
19 : !--------------------------------------------------------------------------
20 : !
21 : ! For filter1:
22 : !
23 : ! This is used by TIEGCM for basic filtering (t,u,v, et.al.),
24 : ! when nlat=72 (2.5 deg res):
25 : !
26 : ! integer,parameter :: kut(nlat) =
27 : ! | (/1 ,1 ,2 ,2 ,4 ,4 ,8 ,8 ,10 ,10 ,12 ,12,
28 : ! | 15 ,15 ,18 ,18 ,22 ,22 ,26 ,26 ,30 ,30 ,32 ,32,
29 : ! | 34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34,
30 : ! | 34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34 ,34,
31 : ! | 32 ,32 ,30 ,30 ,26 ,26 ,22 ,22 ,18 ,18 ,15 ,15,
32 : ! | 12 ,12 ,10 ,10 ,8 ,8 ,4 ,4 ,2 ,2 ,1 ,1/)
33 :
34 : integer, allocatable, protected :: kut1(:)
35 : integer, allocatable, protected :: kut2(:)
36 : integer, allocatable, protected :: nn(:)
37 :
38 : integer, parameter :: nxlat = 96
39 : integer,parameter :: kut1_x(nxlat) = &
40 : (/0 ,0 ,0 ,0 ,1 ,1 ,1 ,1 ,2 ,2 ,2 ,2 , &
41 : 3 ,3 ,3 ,3 ,4 ,4 ,4 ,4 ,6 ,6 ,6 ,6 , &
42 : 8 ,8 ,8 ,8 ,10 ,10 ,10 ,10 ,12 ,12 ,12 ,12, &
43 : 16 ,16 ,18 ,18 ,20 ,20 ,22 ,22 ,24 ,24 ,26 ,26, &
44 : 26 ,26 ,24 ,24 ,22 ,22 ,20 ,20 ,18 ,18 ,16 ,16, &
45 : 12 ,12 ,12 ,12 ,10 ,10 ,10 ,10 ,8 ,8 ,8 ,8 , &
46 : 6 ,6 ,6 ,6 ,4 ,4 ,4 ,4 ,3 ,3 ,3 ,3 , &
47 : 2 ,2 ,2 ,2 ,1 ,1 ,1 ,1 ,0 ,0 ,0 ,0 /)
48 : !--------------------------------------------------------------------------
49 : !
50 : ! For filter2:
51 : !
52 : ! This is used by TIEGCM for O+ filtering when nlat=72 (2.5 deg res):
53 : !
54 : ! kut2=(/0, 0, 1, 2, 4, 4, 6, 6, 8, 8,10,10,12,12,15,15,18,18,
55 : ! | 20,20,20,20,18,18,15,12, 8, 8, 4, 4, 4, 4, 2, 2, 1, 1,
56 : ! | 1, 1, 2, 2, 4, 4, 4, 4, 8, 8,12,15,18,18,20,20,20,20,
57 : ! | 18,18,15,15,12,12,10,10, 8, 8, 6, 6, 4, 4, 2, 1, 0, 0/) ! 2.5 deg
58 : !
59 : ! nn=(/90,90,40,40,22,22,14,14,10,10, 8, 8, 6, 6, 4, 4, 2, 2,
60 : ! | 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
61 : ! | 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
62 : ! | 2, 2, 4, 4, 6, 6, 8, 8,10,10,14,14,22,22,40,40,90,90/) ! 2.5 deg
63 : !
64 : ! At 1.9 deg resolution, nlat==96
65 : !
66 : integer, parameter :: kut2_x(nxlat) = &
67 : (/ 0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 6, 7, &
68 : 8, 8, 9, 10, 11, 11, 13, 14, 15, 16, 17, 18, &
69 : 19, 19, 20, 20, 19, 19, 18, 17, 15, 13, 10, 8, &
70 : 7, 5, 4, 4, 4, 3, 3, 2, 1, 1, 1, 1, &
71 : 1, 1, 1, 1, 2, 3, 3, 4, 4, 4, 5, 7, &
72 : 8, 10, 13, 15, 17, 18, 19, 19, 20, 20, 19, 19, &
73 : 18, 17, 16, 15, 14, 13, 11, 11, 10, 9, 8, 8, &
74 : 7, 6, 5, 5, 4, 3, 2, 1, 0, 0, 0, 0 /)
75 :
76 : integer, parameter :: nn_x(nxlat) = &
77 : (/255, 171, 104, 60, 42, 32, 26, 20, 17, 15, 12, 11, &
78 : 9, 9, 8, 7, 6, 6, 5, 4, 3, 3, 2, 1, &
79 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
80 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
81 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
82 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
83 : 1, 2, 3, 3, 4, 5, 6, 6, 7, 8, 9, 9, &
84 : 11, 12, 15, 17, 20, 26, 32, 42, 60, 104, 171, 255 /)
85 :
86 : ! for ring filter
87 : integer :: nlat_filter
88 : integer, allocatable :: chunk_array(:)
89 :
90 : contains
91 : !-----------------------------------------------------------------------
92 768 : subroutine filter_init(use_ringfilter)
93 : use interpolate_data, only : lininterp
94 :
95 : logical, intent(in) :: use_ringfilter
96 :
97 1536 : real(r8) :: xlats(nxlat), lats(nlat), kut1_out(nlat), kut2_out(nlat)
98 1536 : real(r8) :: nn_out(nlat)
99 : integer :: i
100 : character(len=128) :: errmsg
101 :
102 768 : if (use_ringfilter) then
103 :
104 0 : select case (nlon)
105 : case (72)
106 : ! 5 deg
107 0 : nlat_filter = 6
108 0 : allocate(chunk_array(nlat_filter))
109 : chunk_array(:nlat_filter) = &
110 0 : (/9,18,36,36,72,72/)
111 : case (144)
112 : ! 2.5deg
113 768 : nlat_filter = 16
114 768 : allocate(chunk_array(nlat_filter))
115 : chunk_array(:nlat_filter) = &
116 13056 : (/9,9,9,18,18,18,36,36,36,36,36,72,72,72,72,72/)
117 : case (288)
118 : ! 1.25 deg
119 0 : nlat_filter = 26
120 0 : allocate(chunk_array(nlat_filter))
121 : chunk_array(:nlat_filter) = &
122 : (/9,9,9,18,18,18,36,36,36,36,36,36,72,72,72,72,72,72, &
123 0 : 144,144,144,144,144,144,144,144/)
124 : case (576)
125 : ! 0.625 deg
126 0 : nlat_filter = 45
127 0 : allocate(chunk_array(nlat_filter))
128 : chunk_array(:nlat_filter) = &
129 : (/9,9,9,9,9,18,18,18,18,18,36,36,36,36,36, &
130 : 72,72,72,72,72,72,72,72,72,72, &
131 : 144,144,144,144,144,144,144,144,144,144, &
132 0 : 288,288,288,288,288,288,288,288,288,288/)
133 : case (1152)
134 : ! 0.3125 deg
135 0 : nlat_filter = 111
136 0 : allocate(chunk_array(nlat_filter))
137 : chunk_array(:nlat_filter) = &
138 : (/9,9,9,9,9,9,9,9,9, &
139 : 18,18,18,18,18,18,18,18,18, &
140 : 36,36,36,36,36,36,36,36,36, &
141 : 72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72, &
142 : 144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144,144, &
143 : 288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288,288, &
144 0 : 576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576,576/)
145 : case default
146 0 : write (errmsg,'(a,i6)') 'ringfilter -- resolution is not recognized: number of global longitudes = ',nlon
147 768 : call endrun(errmsg)
148 : end select
149 :
150 : else
151 :
152 0 : ntrigs = 3*nlon/2+1
153 0 : allocate(trigs(ntrigs),kut1(nlat),kut2(nlat),nn(nlat))
154 :
155 0 : call set99(trigs,ifax,nlon) ! initialize fft for O+ polar filtering
156 :
157 0 : do i = 1,nlat
158 0 : lats(i) = -90._r8 + (i-1)*(180._r8/(nlat-1))
159 : enddo
160 :
161 0 : do i = 1,nxlat
162 0 : xlats(i) = -90._r8 + (i-1)*(180._r8/(nxlat-1))
163 : end do
164 :
165 0 : call lininterp( dble(kut1_x), xlats, nxlat, kut1_out, lats, nlat )
166 0 : call lininterp( dble(kut2_x), xlats, nxlat, kut2_out, lats, nlat )
167 :
168 0 : kut1 = int( kut1_out )
169 0 : kut2 = int( kut2_out )
170 :
171 0 : call lininterp( dble(nn_x), xlats, nxlat, nn_out, lats, nlat )
172 0 : nn = int( nn_out )
173 :
174 : endif
175 :
176 768 : end subroutine filter_init
177 : !-----------------------------------------------------------------------
178 0 : subroutine filter1(f,lev0,lev1,lat)
179 : !
180 : ! Remove longitudinal waves of prognostic variables with global fft.
181 : ! Remove wave numbers greater than kut(nlat). This is called after
182 : ! mp_gatherlons, and only by tasks with mytidi==0. On entry, task must
183 : ! have global longitude data defined (mp_gatherlons).
184 : !
185 : ! Args:
186 : integer,intent(in) :: lev0,lev1,lat
187 : real(r8),intent(inout) :: f(nlon,lev0:lev1)
188 : !
189 : ! Local:
190 : integer :: n1,n2,k,i,nlevs
191 0 : real(r8) :: fx(nlon+2,lev1-lev0+1), wfft(nlon+1,lev1-lev0+1)
192 :
193 0 : nlevs = lev1-lev0+1
194 0 : n1 = 2*kut1(lat)+3 ! nyquist freq (?)
195 0 : n2 = nlon+2
196 0 : if (n1 > n2) then
197 : write(iulog,"('filter1: lat=',i2,' kutj=',i2,' n1,2=',2i3,' n1 > n2')") &
198 0 : lat,kut1(lat),n1,n2
199 : return
200 : endif
201 : !
202 : ! Load fx from f for the fft:
203 0 : fx(:,:) = 0._r8
204 0 : do k=lev0,lev1
205 0 : do i=1,nlon
206 0 : fx(i,k) = f(i,k)
207 : enddo
208 : enddo
209 :
210 0 : call fft991(fx, wfft, trigs, ifax,1,nlon+2,nlon,nlevs,-1)
211 : !
212 : ! Remove wave numbers greater than kut(lat)
213 0 : do k = 1,nlevs
214 0 : do i=n1,n2
215 0 : fx(i,k) = 0.0_r8
216 : enddo
217 : enddo
218 : !
219 : ! Inverse transform fourier back to gridpoint:
220 : !
221 0 : call fft991(fx, wfft, trigs, ifax,1,nlon+2,nlon,nlevs,1)
222 : !
223 : ! Redefine f from fx:
224 0 : do k=lev0,lev1
225 0 : do i=1,nlon
226 0 : f(i,k) = fx(i,k)
227 : enddo
228 : enddo
229 : end subroutine filter1
230 : !-----------------------------------------------------------------------
231 0 : subroutine filter2(f,lev0,lev1,lat)
232 : use edyn_geogrid,only : dlamda
233 : !
234 : ! Remove longitudinal waves of prognostic variables with global fft.
235 : ! Remove wave numbers greater than kut2(nlat). This is called after
236 : ! mp_gatherlons, and only by tasks with mytidi==0. On entry, task must
237 : ! have global longitude data defined (mp_gatherlons).
238 : !
239 : ! Args:
240 : integer,intent(in) :: lev0,lev1,lat
241 : real(r8),intent(inout) :: f(nlon,lev0:lev1)
242 : !
243 : ! Local:
244 : integer :: n1,k,i,nlevs
245 0 : real(r8) :: fx(nlon+2,lev1-lev0+1), wfft(nlon+1,lev1-lev0+1)
246 : real(r8) :: smoothfunc,coslon
247 : !
248 0 : nlevs = lev1-lev0+1
249 : !
250 : ! Load local fx from inout f subdomain for the fft:
251 : !
252 0 : fx(:,:) = 0._r8
253 0 : do k=lev0,lev1
254 0 : do i=1,nlon
255 0 : fx(i,k) = f(i,k)
256 : enddo
257 : enddo
258 :
259 0 : call fft991(fx, wfft, trigs, ifax,1,nlon+2,nlon,nlevs,-1)
260 : !
261 : ! Wenbin's comments from TIEGCM:
262 : ! Change filters so that it does not over filtering at high latitudes, it will be the
263 : ! same as filter for low wavenumer, but wrapping up smoothly for large wavenumers, not a
264 : ! sharp transition, so there is still filtering effect in the lower latitudes
265 : ! Wenbin Wang 06/11/13
266 0 : n1=2*(kut2(lat)-1)+3
267 : !
268 : ! Multiply by smoothing function:
269 : ! Test coslon to avoid underflow in smoothfunc at the poles
270 : !
271 :
272 0 : do k=lev0,lev1
273 0 : do i=n1,nlon
274 0 : coslon = cos(((i-n1)/2._r8)*dlamda/2._r8)
275 0 : if ((coslon<0.1_r8) .or. (coslon<0.9_r8 .and. nn(lat)>50)) then
276 0 : fx(i,k) = 0._r8
277 : else
278 0 : smoothfunc = coslon**(2*nn(lat))
279 0 : fx(i,k) = fx(i,k)*smoothfunc
280 : endif
281 : enddo ! i=1,nlon
282 : enddo ! k=lev0,lev1
283 : !
284 : ! Inverse transform fourier back to gridpoint:
285 : !
286 0 : call fft991(fx, wfft, trigs, ifax,1,nlon+2,nlon,nlevs,1)
287 : !
288 : ! Redefine f from fx:
289 0 : do k=lev0,lev1
290 0 : do i=1,nlon
291 0 : f(i,k) = fx(i,k)
292 : enddo
293 : enddo
294 0 : end subroutine filter2
295 : !-----------------------------------------------------------------------
296 :
297 9120 : subroutine ringfilter(lev0,lev1,lat,f)
298 : !
299 : ! Ringfilter for the second order of FFT
300 : ! keep first and second order of fourier series, and filter orders
301 : ! coded by Dang, 2017
302 : ! Args:
303 : integer,intent(in) :: lat,lev0,lev1
304 : real(r8),intent(inout) :: f(nlon,lev0:lev1)
305 : !
306 : ! Local:
307 18240 : real(r8) :: fx(nlon), f_out(nlon)
308 9120 : real(r8), allocatable :: average(:)
309 18240 : real(r8) :: w(nlon),wm1(nlon),a0,a1,b1,theta,dtheta
310 : real(r8) :: fL,fR,fm2,fm1,ff,fp1,fp2,a,b,c
311 : integer :: i,k,points,n,chunk,j_index,m
312 :
313 9120 : near_pole: if(lat .LE. nlat_filter .OR. lat .GE. (nlat-nlat_filter+1) ) then
314 :
315 57760 : allocate(average(maxval(chunk_array)))
316 :
317 3040 : dtheta=2._r8*3.14159_r8/real(nlon)
318 :
319 142880 : do k=lev0,lev1
320 :
321 : ! Load field data into w
322 : ! Fourier expansion: f(x)=a0+a1*cos(x)+b1*sin(x)+others
323 : a1=0._r8
324 : b1=0._r8
325 20276800 : do i=1,nlon
326 20136960 : w(i) = f(i,k)
327 20136960 : theta=dtheta*i
328 20136960 : a1=a1+w(i)*cos(theta)
329 20276800 : b1=b1+w(i)*sin(theta)
330 : enddo
331 139840 : a1=a1*2._r8/real(nlon)
332 139840 : b1=b1*2._r8/real(nlon)
333 20276800 : a0=sum(w)/real(nlon)
334 :
335 : ! Chunk numbers in this latitude
336 139840 : if(lat .LE. nlat_filter) chunk=chunk_array(lat)
337 139840 : if(lat .GE. (nlat-nlat_filter+1)) chunk=chunk_array(nlat-lat+1)
338 :
339 : ! w(i)=wm1(i)+fx(i), then filter fx(i)
340 20276800 : do i=1,nlon
341 20136960 : theta=dtheta*i
342 20136960 : wm1(i)=a0+a1*cos(theta)+b1*sin(theta)
343 : ! fx(i)=w(i)-wm1(i) ! This option is to apply the ring filter to perturbations with wavenumber larger than 1
344 20276800 : fx(i)=w(i) ! This is to apply the ring filter to the full field
345 : enddo
346 :
347 : ! Start the ring average filtering
348 :
349 : ! Grid points in each chunk
350 139840 : points=nlon/chunk
351 139840 : n=points
352 :
353 : ! Calculate the average value in each chunk
354 5567380 : do i=1,chunk ! i is the chunk number in each ring
355 25704340 : average(i)=sum(fx((i-1)*points+1:i*points))/real(points)
356 : enddo
357 :
358 : ! Then do the linear interpolation between each fL, fR
359 5567380 : do i=1,chunk ! i is the chunk number in each ring
360 :
361 : ! Calculate f,fL,fR
362 5427540 : if(i == 1) then
363 :
364 139840 : fm2 = average(chunk-1)
365 139840 : fm1 = average(chunk)
366 139840 : ff = average(i)
367 139840 : fp1 = average(i+1)
368 139840 : fp2 = average(i+2)
369 :
370 5287700 : else if(i == 2) then
371 :
372 139840 : fm2 = average(chunk)
373 139840 : fm1 = average(i-1)
374 139840 : ff = average(i)
375 139840 : fp1 = average(i+1)
376 139840 : fp2 = average(i+2)
377 :
378 5147860 : else if(i == chunk-1) then
379 :
380 139840 : fm2 = average(i-2)
381 139840 : fm1 = average(i-1)
382 139840 : ff = average(i)
383 139840 : fp1 = average(i+1)
384 139840 : fp2 = average(1)
385 :
386 5008020 : else if(i == chunk) then
387 :
388 139840 : fm2 = average(i-2)
389 139840 : fm1 = average(i-1)
390 139840 : ff = average(i)
391 139840 : fp1 = average(1)
392 139840 : fp2 = average(2)
393 :
394 : else
395 :
396 4868180 : fm2 = average(i-2)
397 4868180 : fm1 = average(i-1)
398 4868180 : ff = average(i)
399 4868180 : fp1 = average(i+1)
400 4868180 : fp2 = average(i+2)
401 :
402 : endif
403 :
404 5427540 : fL = (-fm2+7._r8*fm1+7._r8*ff-fp1)/12._r8
405 5427540 : fR = (-fm1+7._r8*ff+7._r8*fp1-fp2)/12._r8
406 :
407 5427540 : a = 3._r8*(fL + fR - 2._r8*ff)
408 5427540 : b = 2._r8*(3._r8*ff - fR - 2._r8*fL)
409 5427540 : c = fL
410 :
411 : ! Calculate the filtered data at j_index
412 25704340 : do m=1,n
413 20136960 : j_index=m+(i-1)*points
414 20136960 : f_out(j_index)=(a/3.0_r8)*(3*m*m-3*m+1)/(n*n) &
415 45701460 : + 0.5_r8*b*(2*m-1)/n + c
416 : enddo
417 :
418 : enddo ! i=1,chunk
419 :
420 20276800 : fx(:)=f_out(:)
421 :
422 : ! Save filtered field:
423 20279840 : do i=1,nlon
424 20276800 : f(i,k) = fx(i) ! This corresponds to the option applying ring filter to the full field.
425 : enddo ! i=1,nlon
426 :
427 : enddo ! k=lev0,lev1
428 :
429 3040 : deallocate(average)
430 :
431 : endif near_pole
432 :
433 9120 : end subroutine ringfilter
434 : !-----------------------------------------------------------------------
435 : end module filter_module
|
