module_polarfft.F
References to this file elsewhere.
1 MODULE module_polarfft
2
3 USE module_model_constants
4 USE module_wrf_error
5 USE module_positive_definite
6
7 CONTAINS
8
9 SUBROUTINE couple_scalars_for_filter ( field &
10 ,mu,mub &
11 ,ids,ide,jds,jde,kds,kde &
12 ,ims,ime,jms,jme,kms,kme &
13 ,ips,ipe,jps,jpe,kps,kpe )
14 IMPLICIT NONE
15 INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde &
16 ,ims,ime,jms,jme,kms,kme &
17 ,ips,ipe,jps,jpe,kps,kpe
18 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: field
19 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
20
21 INTEGER :: i , j , k
22
23 DO j = jps, MIN(jpe,jde-1)
24 DO k = kps, kpe-1
25 DO i = ips, MIN(ipe,ide-1)
26 field(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
27 END DO
28 END DO
29 END DO
30
31 END SUBROUTINE couple_scalars_for_filter
32
33 SUBROUTINE uncouple_scalars_for_filter ( field &
34 ,mu,mub &
35 ,ids,ide,jds,jde,kds,kde &
36 ,ims,ime,jms,jme,kms,kme &
37 ,ips,ipe,jps,jpe,kps,kpe )
38 IMPLICIT NONE
39 INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde &
40 ,ims,ime,jms,jme,kms,kme &
41 ,ips,ipe,jps,jpe,kps,kpe
42 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: field
43 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
44
45 INTEGER :: i , j , k
46
47 DO j = jps, MIN(jpe,jde-1)
48 DO k = kps, kpe-1
49 DO i = ips, MIN(ipe,ide-1)
50 field(i,k,j)=field(i,k,j)/(mu(i,j)+mub(i,j))
51 END DO
52 END DO
53 END DO
54
55 END SUBROUTINE uncouple_scalars_for_filter
56
57 SUBROUTINE pxft ( grid &
58 ,lineno &
59 ,flag_uv,flag_rurv &
60 ,flag_wph,flag_ww &
61 ,flag_t &
62 ,flag_mu,flag_mut &
63 ,flag_moist &
64 ,flag_chem &
65 ,flag_scalar &
66 ,fft_filter_type, dclat &
67 ,positive_definite &
68 ,moist,chem,scalar &
69 ,ids,ide,jds,jde,kds,kde &
70 ,ims,ime,jms,jme,kms,kme &
71 ,ips,ipe,jps,jpe,kps,kpe &
72 ,imsx,imex,jmsx,jmex,kmsx,kmex &
73 ,ipsx,ipex,jpsx,jpex,kpsx,kpex )
74 USE module_domain
75 USE module_tiles
76 USE module_dm
77 IMPLICIT NONE
78 ! Input data.
79 TYPE(domain) , TARGET :: grid
80 integer, intent(in) :: lineno
81 integer myproc, i, j, k
82 LOGICAL, INTENT(IN) :: positive_definite
83 INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde &
84 ,ims,ime,jms,jme,kms,kme &
85 ,ips,ipe,jps,jpe,kps,kpe &
86 ,imsx,imex,jmsx,jmex,kmsx,kmex &
87 ,ipsx,ipex,jpsx,jpex,kpsx,kpex
88 INTEGER, INTENT(IN) :: fft_filter_type
89 REAL, INTENT(IN) :: dclat
90 INTEGER, INTENT(IN) :: flag_uv &
91 ,flag_rurv &
92 ,flag_ww &
93 ,flag_t,flag_wph &
94 ,flag_mu,flag_mut &
95 ,flag_moist &
96 ,flag_chem &
97 ,flag_scalar
98 REAL, DIMENSION(ims:ime,kms:kme,jms:jme,*) , INTENT(INOUT) :: moist, chem, scalar
99
100 ! Local
101 LOGICAL piggyback_mu, piggyback_mut
102 INTEGER ij, k_end
103 #ifdef DM_PARALLEL
104 #else
105 INTEGER itrace
106 #endif
107
108
109 piggyback_mu = flag_mu .EQ. 1
110 piggyback_mut = flag_mut .EQ. 1
111
112 !<DESCRIPTION>
113 !
114 ! The idea is that this parallel transpose fft routine can be
115 ! called at various points in the solver (solve_em) and it will filter
116 ! the correct fields based on the flag arguments. There are two 2d
117 ! fields mu_2 and mut that may need to be filtered too. Since a two-d
118 ! parallel transpose would be inefficient and because the fields that are
119 ! not staggered in z have an extra layer anyway, these can be
120 ! piggybacked. This is handled using latches to makes sure that *somebody*
121 ! carries one or both of these on their back through the filtering and then
122 ! copies them back afterwards. IMPORTANT NOTE: for simplicity, this routine
123 ! is not completely general. It makes the following assumptions:
124 !
125 ! 1) If both flag_mu and flag_mut are specified then flag_uv is also specified
126 !
127 ! 2) If flag_uv is not specified, then only flag_mu and not flag_mut can be
128 !
129 ! 3) If flag_mu is specified than either flag_uv or flag_t must be
130 !
131 ! This is designed to keep the clutter to a minimum in solve_em.
132 ! This is not intended to be a general abstraction of the polar filtering
133 ! calls in in WRF solvers or if the solve_em algorithms change.
134 ! If the needs of the calling solver change, this logic may have to be
135 ! rewritten.
136 !
137 !</DESCRIPTION>
138 !write(0,*)__FILE__,__LINE__,' short circuit '
139 !return
140
141 !write(0,*)'pxft called from ',lineno
142 call wrf_get_myproc(myproc)
143 !write(20+myproc,*)ipex-ipsx+1,jpex-jpsx+1,' clat_xxx '
144 !do j = jpsx, jpex
145 !do i = ipsx, ipex
146 !write(20+myproc,*)grid%clat_xxx(i,j)
147 !enddo
148 !enddo
149
150 !!!!!!!!!!!!!!!!!!!!!!!
151 ! U & V
152 IF ( flag_uv .EQ. 1 ) THEN
153 IF ( piggyback_mu ) THEN
154 grid%em_u_2(ips:ipe,kde,jps:jpe) = grid%em_mu_2(ips:ipe,jps:jpe)
155 ENDIF
156 #ifdef DM_PARALLEL
157 # include "XPOSE_POLAR_FILTER_V_z2x.inc"
158 CALL polar_filter_3d( grid%v_xxx, grid%clat_xxx, .false., &
159 fft_filter_type, dclat, &
160 ids, ide, jds, jde, kds, kde-1, &
161 imsx, imex, jmsx, jmex, kmsx, kmex, &
162 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex ) )
163 # include "XPOSE_POLAR_FILTER_V_x2z.inc"
164 # include "XPOSE_POLAR_FILTER_U_z2x.inc"
165 k_end = MIN(kde-1,kpex)
166 IF ( piggyback_mu ) k_end = MIN(kde,kpex)
167 CALL polar_filter_3d( grid%u_xxx, grid%clat_xxx, piggyback_mu, &
168 fft_filter_type, 0., &
169 ids, ide, jds, jde, kds, kde, &
170 imsx, imex, jmsx, jmex, kmsx, kmex, &
171 ipsx, ipex, jpsx, jpex, kpsx, k_end )
172 # include "XPOSE_POLAR_FILTER_U_x2z.inc"
173 #else
174 CALL polar_filter_3d( grid%em_v_2, grid%clat, .false., &
175 fft_filter_type, dclat, &
176 ids, ide, jds, jde, kds, kde, &
177 ims, ime, jms, jme, kms, kme, &
178 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
179 k_end = MIN(kde-1,kpe)
180 IF ( piggyback_mu ) k_end = MIN(kde,kpe)
181 CALL polar_filter_3d( grid%em_u_2, grid%clat, piggyback_mu, &
182 fft_filter_type, 0., &
183 ids, ide, jds, jde, kds, kde-1, &
184 ims, ime, jms, jme, kms, kme, &
185 ips, ipe, jps, jpe, kps, k_end )
186 #endif
187
188 IF ( piggyback_mu ) THEN
189 grid%em_mu_2(ips:ipe,jps:jpe) = grid%em_u_2(ips:ipe,kde,jps:jpe)
190 piggyback_mu = .FALSE.
191 ENDIF
192 ENDIF
193
194 !!!!!!!!!!!!!!!!!!!!!!!
195 ! T
196 IF ( flag_t .EQ. 1 ) THEN
197 IF ( piggyback_mu ) THEN
198 grid%em_t_2(ips:ipe,kde,jps:jpe) = grid%em_mu_2(ips:ipe,jps:jpe)
199 ENDIF
200 #ifdef DM_PARALLEL
201 # include "XPOSE_POLAR_FILTER_T_z2x.inc"
202 k_end = MIN(kde-1,kpex)
203 IF ( piggyback_mu ) k_end = MIN(kde,kpex)
204 CALL polar_filter_3d( grid%t_xxx, grid%clat_xxx,piggyback_mu, &
205 fft_filter_type, 0., &
206 ids, ide, jds, jde, kds, kde-1, &
207 imsx, imex, jmsx, jmex, kmsx, kmex, &
208 ipsx, ipex, jpsx, jpex, kpsx, k_end )
209 # include "XPOSE_POLAR_FILTER_T_x2z.inc"
210 #else
211 k_end = MIN(kde-1,kpe)
212 IF ( piggyback_mu ) k_end = MIN(kde,kpe)
213 CALL polar_filter_3d( grid%em_t_2, grid%clat, piggyback_mu, &
214 fft_filter_type, 0., &
215 ids, ide, jds, jde, kds, kde-1, &
216 ims, ime, jms, jme, kms, kme, &
217 ips, ipe, jps, jpe, kps, k_end )
218 #endif
219 IF ( piggyback_mu ) THEN
220 grid%em_mu_2(ips:ipe,jps:jpe) = grid%em_t_2(ips:ipe,kde,jps:jpe)
221 piggyback_mu = .FALSE.
222 ENDIF
223 ENDIF
224
225 !!!!!!!!!!!!!!!!!!!!!!!
226 ! W and PH
227 IF ( flag_wph .EQ. 1 ) THEN
228 ! W AND PH USE ALL LEVELS SO NEVER PIGGYBACK, MU IS OUT OF LUCK HERE
229 #ifdef DM_PARALLEL
230 # include "XPOSE_POLAR_FILTER_W_z2x.inc"
231 CALL polar_filter_3d( grid%w_xxx, grid%clat_xxx, .false., &
232 fft_filter_type, 0., &
233 ids, ide, jds, jde, kds, kde, &
234 imsx, imex, jmsx, jmex, kmsx, kmex, &
235 ipsx, ipex, jpsx, jpex, kpsx, kpex )
236 # include "XPOSE_POLAR_FILTER_W_x2z.inc"
237 # include "XPOSE_POLAR_FILTER_PH_z2x.inc"
238 CALL polar_filter_3d( grid%ph_xxx, grid%clat_xxx, .false., &
239 fft_filter_type, 0., &
240 ids, ide, jds, jde, kds, kde, &
241 imsx, imex, jmsx, jmex, kmsx, kmex, &
242 ipsx, ipex, jpsx, jpex, kpsx, kpex )
243 # include "XPOSE_POLAR_FILTER_PH_x2z.inc"
244 #else
245 CALL polar_filter_3d( grid%em_w_2, grid%clat, .false., &
246 fft_filter_type, 0., &
247 ids, ide, jds, jde, kds, kde, &
248 ims, ime, jms, jme, kms, kme, &
249 ips, ipe, jps, jpe, kps, kpe )
250 CALL polar_filter_3d( grid%em_ph_2, grid%clat, .false., &
251 fft_filter_type, 0., &
252 ids, ide, jds, jde, kds, kde, &
253 ims, ime, jms, jme, kms, kme, &
254 ips, ipe, jps, jpe, kps, kpe )
255 #endif
256 ENDIF
257
258 !!!!!!!!!!!!!!!!!!!!!!!
259 ! WW
260 IF ( flag_ww .EQ. 1 ) THEN
261 ! WW USES ALL LEVELS SO NEVER PIGGYBACK, MU IS OUT OF LUCK HERE
262 #ifdef DM_PARALLEL
263 # include "XPOSE_POLAR_FILTER_WW_z2x.inc"
264 CALL polar_filter_3d( grid%ww_xxx, grid%clat_xxx, .false., &
265 fft_filter_type, 0., &
266 ids, ide, jds, jde, kds, kde, &
267 imsx, imex, jmsx, jmex, kmsx, kmex, &
268 ipsx, ipex, jpsx, jpex, kpsx, kpex )
269 # include "XPOSE_POLAR_FILTER_WW_x2z.inc"
270 #else
271 CALL polar_filter_3d( grid%em_ww_m, grid%clat, .false., &
272 fft_filter_type, 0., &
273 ids, ide, jds, jde, kds, kde, &
274 ims, ime, jms, jme, kms, kme, &
275 ips, ipe, jps, jpe, kps, kpe )
276 #endif
277 ENDIF
278
279 !!!!!!!!!!!!!!!!!!!!!!!
280 ! RU AND RV
281 IF ( flag_rurv .EQ. 1 ) THEN
282 IF ( piggyback_mut ) THEN
283 grid%em_ru_m(ips:ipe,kde,jps:jpe) = grid%em_mut(ips:ipe,jps:jpe)
284 ENDIF
285 #ifdef DM_PARALLEL
286 # include "XPOSE_POLAR_FILTER_RV_z2x.inc"
287 CALL polar_filter_3d( grid%rv_xxx, grid%clat_xxx, .false., &
288 fft_filter_type, dclat, &
289 ids, ide, jds, jde, kds, kde, &
290 imsx, imex, jmsx, jmex, kmsx, kmex, &
291 ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
292 # include "XPOSE_POLAR_FILTER_RV_x2z.inc"
293 # include "XPOSE_POLAR_FILTER_RU_z2x.inc"
294 k_end = MIN(kde-1,kpex)
295 IF ( piggyback_mut ) k_end = MIN(kde,kpex)
296 CALL polar_filter_3d( grid%ru_xxx, grid%clat_xxx, piggyback_mut, &
297 fft_filter_type, 0., &
298 ids, ide, jds, jde, kds, kde, &
299 imsx, imex, jmsx, jmex, kmsx, kmex, &
300 ipsx, ipex, jpsx, jpex, kpsx, k_end )
301 #include "XPOSE_POLAR_FILTER_RU_x2z.inc"
302 #else
303 CALL polar_filter_3d( grid%em_rv_m, grid%clat, .false., &
304 fft_filter_type, dclat, &
305 ids, ide, jds, jde, kds, kde, &
306 ims, ime, jms, jme, kms, kme, &
307 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
308 k_end = MIN(kde-1,kpe)
309 IF ( piggyback_mut ) k_end = MIN(kde,kpe)
310 CALL polar_filter_3d( grid%em_ru_m, grid%clat, piggyback_mut, &
311 fft_filter_type, 0., &
312 ids, ide, jds, jde, kds, kde-1, &
313 ims, ime, jms, jme, kms, kme, &
314 ips, ipe, jps, jpe, kps, k_end )
315 #endif
316 IF ( piggyback_mut ) THEN
317 grid%em_mut(ips:ipe,jps:jpe) = grid%em_ru_m(ips:ipe,kde,jps:jpe)
318 piggyback_mut = .FALSE.
319 ENDIF
320 ENDIF
321
322 !!!!!!!!!!!!!!!!!!!!!!!
323 ! MOIST
324 IF ( flag_moist .GE. PARAM_FIRST_SCALAR ) THEN
325 itrace = flag_moist
326 #ifdef DM_PARALLEL
327 # include "XPOSE_POLAR_FILTER_MOIST_z2x.inc"
328 CALL polar_filter_3d( grid%fourd_xxx, grid%clat_xxx, .false. , &
329 fft_filter_type, 0., &
330 ids, ide, jds, jde, kds, kde, &
331 imsx, imex, jmsx, jmex, kmsx, kmex, &
332 ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1), &
333 positive_definite = positive_definite )
334 # include "XPOSE_POLAR_FILTER_MOIST_x2z.inc"
335 #else
336 CALL polar_filter_3d( moist(ims,kms,jms,itrace), grid%clat, .false., &
337 fft_filter_type, 0., &
338 ids, ide, jds, jde, kds, kde, &
339 ims, ime, jms, jme, kms, kme, &
340 ips, ipe, jps, jpe, kps, MIN(kpe,kde-1), &
341 positive_definite = positive_definite )
342 #endif
343 ENDIF
344
345 !!!!!!!!!!!!!!!!!!!!!!!
346 ! CHEM
347 IF ( flag_chem .GE. PARAM_FIRST_SCALAR ) THEN
348 itrace = flag_chem
349 #ifdef DM_PARALLEL
350 # include "XPOSE_POLAR_FILTER_CHEM_z2x.inc"
351 CALL polar_filter_3d( grid%fourd_xxx, grid%clat_xxx, .false. , &
352 fft_filter_type, 0., &
353 ids, ide, jds, jde, kds, kde, &
354 imsx, imex, jmsx, jmex, kmsx, kmex, &
355 ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1), &
356 positive_definite = positive_definite )
357 # include "XPOSE_POLAR_FILTER_MOIST_x2z.inc"
358 #else
359 CALL polar_filter_3d( chem(ims,kms,jms,itrace), grid%clat, .false. , &
360 fft_filter_type, 0., &
361 ids, ide, jds, jde, kds, kde, &
362 ims, ime, jms, jme, kms, kme, &
363 ips, ipe, jps, jpe, kps, MIN(kpe,kde-1), &
364 positive_definite = positive_definite )
365 #endif
366 ENDIF
367
368 !!!!!!!!!!!!!!!!!!!!!!!
369 ! SCALAR
370 IF ( flag_chem .GE. PARAM_FIRST_SCALAR ) THEN
371 itrace = flag_scalar
372 #ifdef DM_PARALLEL
373 # include "XPOSE_POLAR_FILTER_SCALAR_z2x.inc"
374 CALL polar_filter_3d( grid%fourd_xxx , grid%clat_xxx, .false. , &
375 fft_filter_type, 0., &
376 ids, ide, jds, jde, kds, kde, &
377 imsx, imex, jmsx, jmex, kmsx, kmex, &
378 ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1), &
379 positive_definite = positive_definite )
380 # include "XPOSE_POLAR_FILTER_SCALAR_x2z.inc"
381 #else
382 CALL polar_filter_3d( scalar(ims,kms,jms,itrace) , grid%clat, .false. , &
383 fft_filter_type, 0., &
384 ids, ide, jds, jde, kds, kde, &
385 ims, ime, jms, jme, kms, kme, &
386 ips, ipe, jps, jpe, kps, MIN(kpe,kde-1), &
387 positive_definite = positive_definite )
388 #endif
389 ENDIF
390
391 IF ( flag_mu .EQ. 1 .AND. piggyback_mu ) THEN
392 CALL wrf_error_fatal("mu needed to get piggybacked on a transpose and did not")
393 ENDIF
394 IF ( flag_mut .EQ. 1 .AND. piggyback_mut ) THEN
395 CALL wrf_error_fatal("mut needed to get piggybacked on a transpose and did not")
396 ENDIF
397
398 !write(0,*)'pxft back to ',lineno
399 RETURN
400 END SUBROUTINE pxft
401
402 SUBROUTINE polar_filter_3d( f, xlat, piggyback, fft_filter_type, dvlat, &
403 ids, ide, jds, jde, kds, kde, &
404 ims, ime, jms, jme, kms, kme, &
405 its, ite, jts, jte, kts, kte, &
406 positive_definite )
407
408 IMPLICIT NONE
409
410 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
411 ims, ime, jms, jme, kms, kme, &
412 its, ite, jts, jte, kts, kte
413 INTEGER , INTENT(IN ) :: fft_filter_type
414
415 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: f
416 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: xlat
417 REAL , INTENT(IN) :: dvlat
418 LOGICAL , INTENT(IN), OPTIONAL :: positive_definite
419 LOGICAL , INTENT(IN) :: piggyback
420
421 REAL , DIMENSION(1:ide-ids,1:kte-kts+1) :: sheet
422 REAL , DIMENSION(1:kte-kts+1) :: sheet_total
423 REAL :: lat, avg, rnboxw
424 INTEGER :: ig, jg, i, j, j_end, nx, ny, nmax, kw
425 INTEGER :: k, nboxw, nbox2, istart, iend, overlap
426 INTEGER, DIMENSION(6) :: wavenumber = (/ 1, 3, 7, 10, 13, 16 /)
427 REAL, PARAMETER :: filter_latitude = 45.
428
429 ! Variables will stay in domain form since this routine is meaningless
430 ! unless tile extent is the same as domain extent in E/W direction, i.e.,
431 ! the processor has access to all grid points in E/W direction.
432 ! There may be other ways of doing FFTs, but we haven't learned them yet...
433
434 ! Check to make sure we have full access to all E/W points
435 IF ((its /= ids) .OR. (ite /= ide)) THEN
436 WRITE ( wrf_err_message , * ) 'module_polarfft: 3d: (its /= ids) or (ite /= ide)',its,ids,ite,ide
437 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
438 END IF
439
440
441 nx = ide-ids ! "U" stagger variables will be repeated by periodic BCs
442 ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
443 lat = 0.
444 j_end = MIN(jte, jde-1)
445 IF (dvlat /= 0. .and. j_end == jde-1) j_end = jde
446 DO j = jts, j_end
447 ! jg is the J index in the global (domain) span of the array.
448 jg = j-jds+1
449
450 ! determine whether or not to filter the data
451
452 nmax = 0
453 lat = xlat(ids,j)-dvlat
454 if(abs(lat) >= filter_latitude) nmax = -1
455 IF (nmax /= 0) THEN
456 DO k=kts,kte
457 DO i=ids,ide-1
458 sheet(i-ids+1,k-kts+1) = f(i,k,j)
459 END DO
460 END DO
461
462 IF (PRESENT(positive_definite)) THEN
463 IF (positive_definite) THEN
464 DO k=1,ny
465 sheet_total(k) = SUM(sheet(:,k))
466 END DO
467 END IF
468 END IF
469
470 CALL polar_filter_fft_2d_ncar(nx,ny,sheet,nmax,lat,filter_latitude,piggyback)
471
472 IF (PRESENT(positive_definite)) THEN
473 IF (positive_definite) THEN
474 CALL positive_definite_sheet(sheet,sheet_total,nx,ny)
475 END IF
476 END IF
477
478 DO k=kts,kte
479 DO i=ids,ide-1
480 f(i,k,j) = sheet(i-ids+1,k-kts+1)
481 END DO
482 ! setting up ims-ime with x periodicity:
483 ! enforce periodicity as in set_physical_bc3d
484 DO i=1,ids-ims
485 f(ids-i,k,j)=f(ide-i,k,j)
486 END DO
487 DO i=1,ime-ide+1
488 f(ide+i-1,k,j)=f(ids+i-1,k,j)
489 END DO
490 END DO
491 END IF
492 END DO
493
494 END SUBROUTINE polar_filter_3d
495
496 !------------------------------------------------------------------------------
497
498 SUBROUTINE polar_filter_fft_2d_ncar(nx,ny,fin,nmax,lat,filter_latitude,piggyback)
499 IMPLICIT NONE
500 INTEGER , INTENT(IN) :: nx, ny, nmax
501 REAL , DIMENSION(nx,ny), INTENT(INOUT) :: fin
502 REAL , INTENT(IN) :: lat, filter_latitude
503 LOGICAL, INTENT(IN) :: piggyback
504
505 REAL :: pi, rcosref, freq, c, cf
506 INTEGER :: i, j
507 REAL, dimension(nx,ny) :: fp
508
509 INTEGER :: lensave, ier, nh, n1
510 INTEGER :: lot, jump, n, inc, lenr, lensav, lenwrk
511 REAL, DIMENSION(nx+15) :: wsave
512 REAL, DIMENSION(nx,ny) :: work
513 REAL, PARAMETER :: alpha = 0.0
514 REAL :: factor_k
515
516 INTEGER :: ntop
517
518 ! *** only transform if nmax < 0 ***
519
520 IF( nmax < 0 ) THEN
521
522 pi = ACOS(-1.)
523 rcosref = 1./COS(filter_latitude*pi/180.)
524
525 ! we are following the naming convention of the fftpack5 routines
526
527 n = nx
528 lot = ny
529 lensav = n+15
530 inc = 1
531 lenr = nx*ny
532 jump = nx
533 lenwrk = lenr
534 ntop = ny
535 IF(piggyback) ntop = ny-1
536
537 ! forward transform
538 ! initialize coefficients, place in wsave
539 ! (should place this in init and save wsave at program start)
540
541 call rfftmi(n,wsave,lensav,ier)
542 IF(ier /= 0) THEN
543 write(0,*) ' error in rfftmi ',ier
544 END IF
545
546 ! do the forward transform
547
548 call rfftmf( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
549 IF(ier /= 0) THEN
550 write(0,*) ' error in rfftmf ',ier
551 END IF
552
553 if(MOD(n,2) == 0) then
554 nh = n/2 - 1
555 else
556 nh = (n-1)/2
557 end if
558
559 DO j=1,ny
560 fp(1,j) = 1.
561 ENDDO
562
563 DO i=2,nh+1
564 freq=REAL(i-1)/REAL(n)
565 c = (rcosref*COS(lat*pi/180.)/SIN(freq*pi))**2
566 ! c = MAX(0.,MIN(1.,c))
567 do j=1,ntop
568 factor_k = (1.-alpha)+alpha*min(1.,float(ntop - j)/10.)
569 cf = c*factor_k*factor_k
570 cf = MAX(0.,MIN(1.,cf))
571 fp(2*(i-1),j) = cf
572 fp(2*(i-1)+1,j) = cf
573 enddo
574 if(piggyback) then
575 cf = MAX(0.,MIN(1.,c))
576 fp(2*(i-1),ny) = cf
577 fp(2*(i-1)+1,ny) = cf
578 endif
579 END DO
580
581 IF(MOD(n,2) == 0) THEN
582 c = (rcosref*COS(lat*pi/180.))**2
583 ! c = MAX(0.,MIN(1.,c))
584 do j=1,ntop
585 factor_k = (1.-alpha)+alpha*min(1.,float(ntop - j)/10.)
586 cf = c*factor_k*factor_k
587 cf = MAX(0.,MIN(1.,cf))
588 fp(n,j) = cf
589 enddo
590 if(piggyback) then
591 cf = MAX(0.,MIN(1.,c))
592 fp(n,ny) = cf
593 endif
594 END IF
595
596 DO j=1,ny
597 DO i=1,nx
598 fin(i,j) = fp(i,j)*fin(i,j)
599 ENDDO
600 ENDDO
601
602 ! do the backward transform
603
604 call rfftmb( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
605 IF(ier /= 0) THEN
606 write(0,*) ' error in rfftmb ',ier
607 END IF
608
609 END IF
610
611 END SUBROUTINE polar_filter_fft_2d_ncar
612
613 !------------------------------------------------------------------------------
614
615 END MODULE module_polarfft
616