Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_polarfft.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 22.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format