Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_polarfft.F @ 3094

Last change on this file since 3094 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

Rev	Line
[2759]	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