Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_polarfft.F @ 1313

Last change on this file since 1313 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format