Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dfi.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: 51.8 KB

Rev	Line
[2759]	1	SUBROUTINE dfi_accumulate( grid )
	2
	3	USE module_domain
	4	USE module_configure
	5	USE module_driver_constants
	6	USE module_machine
	7	USE module_dm
	8	USE module_model_constants
	9	USE module_state_description
	10
	11	IMPLICIT NONE
	12
	13	REAL hn
	14
	15	! Input data.
	16
	17	TYPE(domain) , POINTER :: grid
	18
	19	#if (EM_CORE == 1)
	20
	21	IF ( grid%dfi_opt .EQ. DFI_NODFI .OR. grid%dfi_stage .EQ. DFI_FST ) RETURN
	22
	23	hn = grid%hcoeff(grid%itimestep+1)
	24
	25	! accumulate dynamic variables
	26	grid%dfi_mu(:,:) = grid%dfi_mu(:,:) + grid%mu_2(:,:) * hn
	27	grid%dfi_u(:,:,:) = grid%dfi_u(:,:,:) + grid%u_2(:,:,:) * hn
	28	grid%dfi_v(:,:,:) = grid%dfi_v(:,:,:) + grid%v_2(:,:,:) * hn
	29	grid%dfi_w(:,:,:) = grid%dfi_w(:,:,:) + grid%w_2(:,:,:) * hn
	30	grid%dfi_ww(:,:,:) = grid%dfi_ww(:,:,:) + grid%ww(:,:,:) * hn
	31	grid%dfi_t(:,:,:) = grid%dfi_t(:,:,:) + grid%t_2(:,:,:) * hn
	32	grid%dfi_phb(:,:,:) = grid%dfi_phb(:,:,:) + grid%phb(:,:,:) * hn
	33	grid%dfi_ph0(:,:,:) = grid%dfi_ph0(:,:,:) + grid%ph0(:,:,:) * hn
	34	grid%dfi_php(:,:,:) = grid%dfi_php(:,:,:) + grid%php(:,:,:) * hn
	35	grid%dfi_p(:,:,:) = grid%dfi_p(:,:,:) + grid%p(:,:,:) * hn
	36	grid%dfi_ph(:,:,:) = grid%dfi_ph(:,:,:) + grid%ph_2(:,:,:) * hn
	37	grid%dfi_tke(:,:,:) = grid%dfi_tke(:,:,:) + grid%tke_2(:,:,:) * hn
	38	grid%dfi_al(:,:,:) = grid%dfi_al(:,:,:) + grid%al(:,:,:) * hn
	39	grid%dfi_alt(:,:,:) = grid%dfi_alt(:,:,:) + grid%alt(:,:,:) * hn
	40	grid%dfi_pb(:,:,:) = grid%dfi_pb(:,:,:) + grid%pb(:,:,:) * hn
	41	grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
	42	grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
	43
	44	! accumulate DFI coefficient
	45	grid%hcoeff_tot = grid%hcoeff_tot + hn
	46	#endif
	47
	48	END SUBROUTINE dfi_accumulate
	49
	50
	51	#if (EM_CORE == 1)
	52	SUBROUTINE wrf_dfi_bck_init ( )
	53
	54	USE module_domain
	55	USE module_utility
	56
	57	IMPLICIT NONE
	58
	59	INTEGER rc
	60
	61	INTERFACE
	62	SUBROUTINE Setup_Timekeeping(grid)
	63	USE module_domain
	64	TYPE (domain), POINTER :: grid
	65	END SUBROUTINE Setup_Timekeeping
	66
	67	SUBROUTINE dfi_save_arrays(grid)
	68	USE module_domain
	69	TYPE (domain), POINTER :: grid
	70	END SUBROUTINE dfi_save_arrays
	71
	72	SUBROUTINE dfi_clear_accumulation(grid)
	73	USE module_domain
	74	TYPE (domain), POINTER :: grid
	75	END SUBROUTINE dfi_clear_accumulation
	76
	77	SUBROUTINE optfil_driver(grid)
	78	USE module_domain
	79	TYPE (domain), POINTER :: grid
	80	END SUBROUTINE optfil_driver
	81
	82	SUBROUTINE start_domain(grid, allowed_to_read)
	83	USE module_domain
	84	TYPE (domain) :: grid
	85	LOGICAL, INTENT(IN) :: allowed_to_read
	86	END SUBROUTINE start_domain
	87	END INTERFACE
	88
	89	head_grid%dfi_stage = DFI_BCK
	90
	91	! Negate time step
	92	CALL nl_set_time_step ( 1, -head_grid%time_step )
	93
	94	CALL Setup_Timekeeping (head_grid)
	95
	96	! set physics options to zero
	97	CALL nl_set_mp_physics( 1, 0 )
	98	CALL nl_set_ra_lw_physics( 1, 0 )
	99	CALL nl_set_ra_sw_physics( 1, 0 )
	100	CALL nl_set_sf_surface_physics( 1, 0 )
	101	CALL nl_set_sf_sfclay_physics( 1, 0 )
	102	CALL nl_set_bl_pbl_physics( 1, 0 )
	103	CALL nl_set_cu_physics( 1, 0 )
	104
	105	! set diffusion to zero for backward integration
	106	CALL nl_set_km_opt( 1, head_grid%km_opt_dfi)
	107	CALL nl_set_pd_moist( 1, head_grid%pd_moist_dfi)
	108
	109	head_grid%start_subtime = domain_get_start_time ( head_grid )
	110	head_grid%stop_subtime = domain_get_stop_time ( head_grid )
	111
	112	CALL WRFU_ClockSet(head_grid%domain_clock, currTime=head_grid%start_subtime, rc=rc)
	113
	114	CALL dfi_save_arrays ( head_grid )
	115	CALL dfi_clear_accumulation( head_grid )
	116	CALL optfil_driver(head_grid)
	117
	118	!tgs need to call start_domain here to reset bc initialization for negative dt
	119	CALL start_domain ( head_grid , .TRUE. )
	120
	121	END SUBROUTINE wrf_dfi_bck_init
	122
	123
	124	SUBROUTINE wrf_dfi_fwd_init ( )
	125
	126	USE module_domain
	127	USE module_utility
	128
	129	IMPLICIT NONE
	130
	131	INTEGER rc
	132
	133	INTERFACE
	134	SUBROUTINE Setup_Timekeeping(grid)
	135	USE module_domain
	136	TYPE (domain), POINTER :: grid
	137	END SUBROUTINE Setup_Timekeeping
	138
	139	SUBROUTINE dfi_save_arrays(grid)
	140	USE module_domain
	141	TYPE (domain), POINTER :: grid
	142	END SUBROUTINE dfi_save_arrays
	143
	144	SUBROUTINE dfi_clear_accumulation(grid)
	145	USE module_domain
	146	TYPE (domain), POINTER :: grid
	147	END SUBROUTINE dfi_clear_accumulation
	148
	149	SUBROUTINE optfil_driver(grid)
	150	USE module_domain
	151	TYPE (domain), POINTER :: grid
	152	END SUBROUTINE optfil_driver
	153
	154	SUBROUTINE start_domain(grid, allowed_to_read)
	155	USE module_domain
	156	TYPE (domain) :: grid
	157	LOGICAL, INTENT(IN) :: allowed_to_read
	158	END SUBROUTINE start_domain
	159	END INTERFACE
	160
	161	head_grid%dfi_stage = DFI_FWD
	162
	163	! get the negative time step from the namelist and store
	164	! it as positive again.
	165	! for Setup_Timekeeping to use when setting the clock
	166	! note that this ignores fractional parts of time step
	167	CALL nl_get_time_step( 1, head_grid%time_step )
	168	head_grid%time_step = abs(head_grid%time_step)
	169	CALL nl_set_time_step( 1, head_grid%time_step )
	170
	171	head_grid%itimestep=0
	172	head_grid%xtime=0.
	173
	174	! reset physics options to normal
	175	CALL nl_set_mp_physics( 1, head_grid%mp_physics)
	176	CALL nl_set_ra_lw_physics( 1, head_grid%ra_lw_physics)
	177	CALL nl_set_ra_sw_physics( 1, head_grid%ra_sw_physics)
	178	CALL nl_set_sf_surface_physics( 1, head_grid%sf_surface_physics)
	179	CALL nl_set_sf_sfclay_physics( 1, head_grid%sf_sfclay_physics)
	180	CALL nl_set_bl_pbl_physics( 1, head_grid%bl_pbl_physics)
	181	CALL nl_set_cu_physics( 1, head_grid%cu_physics)
	182
	183	! reset km_opt to norlmal
	184	CALL nl_set_km_opt( 1, head_grid%km_opt)
	185	CALL nl_set_pd_moist( 1, head_grid%pd_moist)
	186
	187	CALL Setup_Timekeeping (head_grid)
	188	head_grid%start_subtime = domain_get_start_time ( head_grid )
	189	head_grid%stop_subtime = domain_get_stop_time ( head_grid )
	190
	191	CALL WRFU_ClockSet(head_grid%domain_clock, currTime=head_grid%start_subtime, rc=rc)
	192
	193	IF ( head_grid%dfi_opt .EQ. DFI_DFL ) THEN
	194	CALL dfi_save_arrays ( head_grid )
	195	END IF
	196	CALL dfi_clear_accumulation( head_grid )
	197	CALL optfil_driver(head_grid)
	198
	199	!tgs need to call it here to reset bc initialization for positive time_step
	200	CALL start_domain ( head_grid , .TRUE. )
	201
	202	END SUBROUTINE wrf_dfi_fwd_init
	203
	204
	205	SUBROUTINE wrf_dfi_fst_init ( )
	206
	207	USE module_domain
	208
	209	IMPLICIT NONE
	210
	211	CHARACTER (LEN=80) :: wrf_error_message
	212
	213	INTERFACE
	214	SUBROUTINE Setup_Timekeeping(grid)
	215	USE module_domain
	216	TYPE (domain), POINTER :: grid
	217	END SUBROUTINE Setup_Timekeeping
	218
	219	SUBROUTINE dfi_save_arrays(grid)
	220	USE module_domain
	221	TYPE (domain), POINTER :: grid
	222	END SUBROUTINE dfi_save_arrays
	223
	224	SUBROUTINE dfi_clear_accumulation(grid)
	225	USE module_domain
	226	TYPE (domain), POINTER :: grid
	227	END SUBROUTINE dfi_clear_accumulation
	228
	229	SUBROUTINE optfil_driver(grid)
	230	USE module_domain
	231	TYPE (domain), POINTER :: grid
	232	END SUBROUTINE optfil_driver
	233
	234	SUBROUTINE start_domain(grid, allowed_to_read)
	235	USE module_domain
	236	TYPE (domain) :: grid
	237	LOGICAL, INTENT(IN) :: allowed_to_read
	238	END SUBROUTINE start_domain
	239	END INTERFACE
	240
	241	head_grid%dfi_stage = DFI_FST
	242
	243	head_grid%itimestep=0
	244	head_grid%xtime=0. ! BUG: This will probably not work for all DFI options
	245
	246	CALL Setup_Timekeeping (head_grid)
	247	head_grid%start_subtime = domain_get_start_time ( head_grid )
	248	head_grid%stop_subtime = domain_get_stop_time ( head_grid )
	249
	250	CALL start_domain ( head_grid , .TRUE. )
	251
	252	END SUBROUTINE wrf_dfi_fst_init
	253
	254
	255	SUBROUTINE wrf_dfi_write_initialized_state( )
	256
	257	! Driver layer
	258	USE module_domain
	259	USE module_io_domain
	260	USE module_configure
	261
	262	IMPLICIT NONE
	263
	264	INTEGER :: fid, ierr
	265	CHARACTER (LEN=80) :: wrf_error_message
	266	CHARACTER (LEN=80) :: rstname
	267	CHARACTER (LEN=132) :: message
	268
	269	TYPE (grid_config_rec_type) :: config_flags
	270
	271	CALL model_to_grid_config_rec ( head_grid%id , model_config_rec , config_flags )
	272
	273	WRITE (wrf_err_message,'(A,I4)') 'Writing out initialized model state'
	274	CALL wrf_message(TRIM(wrf_err_message))
	275
	276	rstname = 'wrfinput_initialized_d01'
	277	CALL open_w_dataset ( fid, TRIM(rstname), head_grid, config_flags, output_model_input, "DATASET=INPUT", ierr )
	278	IF ( ierr .NE. 0 ) THEN
	279	WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
	280	CALL WRF_ERROR_FATAL ( message )
	281	END IF
	282	CALL output_model_input ( fid, head_grid, config_flags, ierr )
	283	CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
	284
	285	END SUBROUTINE wrf_dfi_write_initialized_state
	286
	287
	288	SUBROUTINE wrf_dfi_array_reset ( )
	289
	290	USE module_domain
	291
	292	IMPLICIT NONE
	293
	294	INTERFACE
	295	SUBROUTINE dfi_array_reset(grid)
	296	USE module_domain
	297	TYPE (domain), POINTER :: grid
	298	END SUBROUTINE dfi_array_reset
	299	END INTERFACE
	300
	301	! Copy filtered arrays back into state arrays in grid structure, and
	302	! restore original land surface fields
	303	CALL dfi_array_reset( head_grid )
	304
	305	END SUBROUTINE wrf_dfi_array_reset
	306
	307
	308	SUBROUTINE optfil_driver( grid )
	309
	310	USE module_domain
	311	USE module_wrf_error
	312	USE module_timing
	313	USE module_date_time
	314	USE module_configure
	315	USE module_state_description
	316	USE module_model_constants
	317
	318	IMPLICIT NONE
	319
	320	TYPE (domain) , POINTER :: grid
	321
	322	! Local variables
	323	integer :: nstep2, nstepmax, rundfi, i, rc
	324	real :: timestep, tauc
	325	TYPE(WRFU_TimeInterval) :: run_interval
	326
	327	timestep=abs(grid%dt)
	328	run_interval = grid%stop_subtime - grid%start_subtime
	329
	330	CALL WRFU_TimeIntervalGet( run_interval, S=rundfi, rc=rc )
	331	rundfi = abs(rundfi)
	332
	333	nstep2= ceiling((1.0 + real(rundfi)/timestep) / 2.0)
	334
	335	! nstep2 is equal to a half of timesteps per initialization period,
	336	! should not exceed nstepmax
	337
	338	tauc = real(grid%dfi_cutoff_seconds)
	339
	340	! Get DFI coefficient
	341	grid%hcoeff(:) = 0.0
	342	IF ( grid%dfi_nfilter < 0 .OR. grid%dfi_nfilter > 8 ) THEN
	343	write(0,*) 'Invalid filter specified in namelist.'
	344	ELSE
	345	call dfcoef(nstep2-1, grid%dt, tauc, grid%dfi_nfilter, grid%hcoeff)
	346	END IF
	347
	348	IF ( MOD(int(1.0 + real(rundfi)/timestep),2) /= 0 ) THEN
	349	DO i=1,nstep2-1
	350	grid%hcoeff(2*nstep2-i) = grid%hcoeff(i)
	351	END DO
	352	ELSE
	353	DO i=1,nstep2
	354	grid%hcoeff(2*nstep2-i+1) = grid%hcoeff(i)
	355	END DO
	356	END IF
	357
	358	END SUBROUTINE optfil_driver
	359
	360
	361	SUBROUTINE dfi_clear_accumulation( grid )
	362
	363	USE module_domain
	364	USE module_configure
	365	USE module_driver_constants
	366	USE module_machine
	367	USE module_dm
	368	USE module_model_constants
	369	USE module_state_description
	370
	371	IMPLICIT NONE
	372
	373	! Input data.
	374	TYPE(domain) , POINTER :: grid
	375
	376	grid%dfi_mu(:,:) = 0.
	377	grid%dfi_u(:,:,:) = 0.
	378	grid%dfi_v(:,:,:) = 0.
	379	grid%dfi_w(:,:,:) = 0.
	380	grid%dfi_ww(:,:,:) = 0.
	381	grid%dfi_t(:,:,:) = 0.
	382	grid%dfi_phb(:,:,:) = 0.
	383	grid%dfi_ph0(:,:,:) = 0.
	384	grid%dfi_php(:,:,:) = 0.
	385	grid%dfi_p(:,:,:) = 0.
	386	grid%dfi_ph(:,:,:) = 0.
	387	grid%dfi_tke(:,:,:) = 0.
	388	grid%dfi_al(:,:,:) = 0.
	389	grid%dfi_alt(:,:,:) = 0.
	390	grid%dfi_pb(:,:,:) = 0.
	391	grid%dfi_moist(:,:,:,:) = 0.
	392	grid%dfi_scalar(:,:,:,:) = 0.
	393
	394	grid%hcoeff_tot = 0.0
	395
	396	END SUBROUTINE dfi_clear_accumulation
	397
	398
	399	SUBROUTINE dfi_save_arrays( grid )
	400
	401	USE module_domain
	402	USE module_configure
	403	USE module_driver_constants
	404	USE module_machine
	405	USE module_dm
	406	USE module_model_constants
	407	USE module_state_description
	408
	409	IMPLICIT NONE
	410
	411	! Input data.
	412	TYPE(domain) , POINTER :: grid
	413
	414	! save surface 2-D fields
	415	grid%dfi_SNOW(:,:) = grid%SNOW(:,:)
	416	grid%dfi_SNOWH(:,:) = grid%SNOWH(:,:)
	417	grid%dfi_SNOWC(:,:) = grid%SNOWC(:,:)
	418	grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
	419	grid%dfi_TSK(:,:) = grid%TSK(:,:)
	420	grid%dfi_QVG(:,:) = grid%QVG(:,:)
	421	grid%dfi_SOILT1(:,:) = grid%SOILT1(:,:)
	422	grid%dfi_TSNAV(:,:) = grid%TSNAV(:,:)
	423
	424	! save soil fields
	425	grid%dfi_TSLB(:,:,:) = grid%TSLB(:,:,:)
	426	grid%dfi_SMOIS(:,:,:) = grid%SMOIS(:,:,:)
	427	grid%dfi_SMFR3D(:,:,:) = grid%SMFR3D(:,:,:)
	428	grid%dfi_KEEPFR3DFLAG(:,:,:) = grid%KEEPFR3DFLAG(:,:,:)
	429
	430	END SUBROUTINE dfi_save_arrays
	431
	432
	433	SUBROUTINE dfi_array_reset( grid )
	434
	435	USE module_domain
	436	USE module_configure
	437	USE module_driver_constants
	438	USE module_machine
	439	USE module_dm
	440	USE module_model_constants
	441	USE module_state_description
	442
	443	IMPLICIT NONE
	444
	445	! Input data.
	446	TYPE(domain) , POINTER :: grid
	447
	448	IF ( grid%dfi_opt .EQ. DFI_NODFI ) RETURN
	449
	450
	451	! Set dynamic variables
	452	! divide by total DFI coefficient
	453	grid%mu_2(:,:) = grid%dfi_mu(:,:) / grid%hcoeff_tot
	454	grid%u_2(:,:,:) = grid%dfi_u(:,:,:) / grid%hcoeff_tot
	455	grid%v_2(:,:,:) = grid%dfi_v(:,:,:) / grid%hcoeff_tot
	456	grid%w_2(:,:,:) = grid%dfi_w(:,:,:) / grid%hcoeff_tot
	457	grid%ww(:,:,:) = grid%dfi_ww(:,:,:) / grid%hcoeff_tot
	458	grid%t_2(:,:,:) = grid%dfi_t(:,:,:) / grid%hcoeff_tot
	459	grid%phb(:,:,:) = grid%dfi_phb(:,:,:) / grid%hcoeff_tot
	460	grid%ph0(:,:,:) = grid%dfi_ph0(:,:,:) / grid%hcoeff_tot
	461	grid%php(:,:,:) = grid%dfi_php(:,:,:) / grid%hcoeff_tot
	462	grid%p(:,:,:) = grid%dfi_p(:,:,:) / grid%hcoeff_tot
	463	grid%ph_2(:,:,:) = grid%dfi_ph(:,:,:) / grid%hcoeff_tot
	464	grid%tke_2(:,:,:) = grid%dfi_tke(:,:,:) / grid%hcoeff_tot
	465	grid%al(:,:,:) = grid%dfi_al(:,:,:) / grid%hcoeff_tot
	466	grid%alt(:,:,:) = grid%dfi_alt(:,:,:) / grid%hcoeff_tot
	467	grid%pb(:,:,:) = grid%dfi_pb(:,:,:) / grid%hcoeff_tot
	468	grid%moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot
	469	grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot
	470
	471	! restore initial fields
	472	grid%SNOW (:,:) = grid%dfi_SNOW (:,:)
	473	grid%SNOWH (:,:) = grid%dfi_SNOWH (:,:)
	474	grid%SNOWC (:,:) = grid%dfi_SNOWC (:,:)
	475	grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
	476	grid%TSK (:,:) = grid%dfi_TSK (:,:)
	477	grid%QVG (:,:) = grid%dfi_QVG (:,:)
	478	grid%SOILT1(:,:) = grid%dfi_SOILT1(:,:)
	479	grid%TSNAV (:,:) = grid%dfi_TSNAV (:,:)
	480
	481	grid%TSLB (:,:,:) = grid%dfi_TSLB (:,:,:)
	482	grid%SMOIS (:,:,:) = grid%dfi_SMOIS (:,:,:)
	483	grid%SMFR3D(:,:,:) = grid%dfi_SMFR3D (:,:,:)
	484	grid%KEEPFR3DFLAG(:,:,:) = grid%dfi_KEEPFR3DFLAG(:,:,:)
	485
	486	END SUBROUTINE dfi_array_reset
	487
	488
	489	SUBROUTINE dfcoef (NSTEPS,DT,TAUC,NFILT,H)
	490	!
	491	! calculate filter weights with selected window.
	492	!
	493	! peter lynch and xiang-yu huang
	494	!
	495	! ref: see hamming, r.w., 1989: digital filters,
	496	! prentice-hall international. 3rd edition.
	497	!
	498	! input: nsteps - number of timesteps
	499	! forward or backward.
	500	! dt - time step in seconds.
	501	! tauc - cut-off period in seconds.
	502	! nfilt - indicator for selected filter.
	503	!
	504	! output: h - array(0:nsteps) with the
	505	! required filter weights
	506	!
	507	!------------------------------------------------------------
	508
	509	implicit none
	510
	511	integer, intent(in) :: nsteps, nfilt
	512	real , intent(in) :: dt, tauc
	513	real, intent(out) :: h(1:nsteps+1)
	514
	515	! Local data
	516
	517	integer :: n
	518	real :: pi, omegac, x, window, deltat
	519	real :: hh(0:nsteps)
	520
	521	pi=4*ATAN(1.)
	522	deltat=ABS(dt)
	523
	524	! windows are defined by a call and held in hh.
	525
	526	if ( nfilt .eq. -1) call debug (nsteps,h)
	527
	528	IF ( NFILT .EQ. 0 ) CALL UNIFORM (NSTEPS,HH)
	529	IF ( NFILT .EQ. 1 ) CALL LANCZOS (NSTEPS,HH)
	530	IF ( NFILT .EQ. 2 ) CALL HAMMING (NSTEPS,HH)
	531	IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
	532	IF ( NFILT .EQ. 4 ) CALL KAISER (NSTEPS,HH)
	533	IF ( NFILT .EQ. 5 ) CALL POTTER2 (NSTEPS,HH)
	534	IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
	535
	536	IF ( NFILT .LE. 6 ) THEN ! sinc-windowed filters
	537
	538	! calculate the cutoff frequency
	539	OMEGAC = 2.*PI/TAUC
	540
	541	DO N=0,NSTEPS
	542	WINDOW = HH(N)
	543	IF ( N .EQ. 0 ) THEN
	544	X = (OMEGAC*DELTAT/PI)
	545	ELSE
	546	X = SIN(NOMEGACDELTAT)/(N*PI)
	547	END IF
	548	HH(N) = X*WINDOW
	549	END DO
	550
	551	! normalize the sums to be unity
	552	CALL NORMLZ(HH,NSTEPS)
	553
	554	DO N=0,NSTEPS
	555	H(N+1) = HH(NSTEPS-N)
	556	END DO
	557
	558	ELSE IF ( NFILT .EQ. 7 ) THEN ! dolph filter
	559
	560	CALL DOLPH(DT,TAUC,NSTEPS,H)
	561
	562	ELSE IF ( NFILT .EQ. 8 ) THEN ! 2nd order, 2nd type quick start filter (RHO)
	563
	564	CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS)
	565
	566	END IF
	567
	568	RETURN
	569
	570	END SUBROUTINE dfcoef
	571
	572
	573	SUBROUTINE NORMLZ(HH,NMAX)
	574
	575	! normalize the sum of hh to be unity
	576
	577	implicit none
	578
	579	integer, intent(in) :: nmax
	580	real , dimension(0:nmax), intent(out) :: hh
	581
	582	! local data
	583	real :: sumhh
	584	integer :: n
	585
	586	SUMHH = HH(0)
	587	DO N=1,NMAX
	588	SUMHH = SUMHH + 2*HH(N)
	589	ENDDO
	590	DO N=0,NMAX
	591	HH(N) = HH(N)/SUMHH
	592	ENDDO
	593
	594	RETURN
	595
	596	END subroutine normlz
	597
	598
	599	subroutine debug(nsteps, ww)
	600
	601	implicit none
	602
	603	integer, intent(in) :: nsteps
	604	real , dimension(0:nsteps), intent(out) :: ww
	605	integer :: n
	606
	607	do n=0,nsteps
	608	ww(n)=0
	609	end do
	610
	611	ww(int(nsteps/2))=1
	612
	613	return
	614
	615	end subroutine debug
	616
	617
	618	SUBROUTINE UNIFORM(NSTEPS,WW)
	619
	620	! define uniform or rectangular window function.
	621
	622	implicit none
	623
	624	integer, intent(in) :: nsteps
	625	real , dimension(0:nsteps), intent(out) :: ww
	626
	627	integer :: n
	628
	629	DO N=0,NSTEPS
	630	WW(N) = 1.
	631	ENDDO
	632
	633	RETURN
	634
	635	END subroutine uniform
	636
	637
	638	SUBROUTINE LANCZOS(NSTEPS,WW)
	639
	640	! define (genaralised) lanczos window function.
	641
	642	implicit none
	643
	644	integer, parameter :: nmax = 1000
	645	integer, intent(in) :: nsteps
	646	real , dimension(0:nmax), intent(out) :: ww
	647	integer :: n
	648	real :: power, pi, w
	649
	650	! (for the usual lanczos window, power = 1 )
	651	POWER = 1
	652
	653	PI=4*ATAN(1.)
	654	DO N=0,NSTEPS
	655	IF ( N .EQ. 0 ) THEN
	656	W = 1.0
	657	ELSE
	658	W = SIN(NPI/(NSTEPS+1)) / ( NPI/(NSTEPS+1))
	659	ENDIF
	660	WW(N) = W**POWER
	661	ENDDO
	662
	663	RETURN
	664
	665	END SUBROUTINE lanczos
	666
	667
	668	SUBROUTINE HAMMING(NSTEPS,WW)
	669
	670	! define (genaralised) hamming window function.
	671
	672	implicit none
	673
	674	integer, intent(in) :: nsteps
	675	real, dimension(0:nsteps) :: ww
	676	integer :: n
	677	real :: alpha, pi, w
	678
	679	! (for the usual hamming window, alpha=0.54,
	680	! for the hann window, alpha=0.50).
	681	ALPHA=0.54
	682
	683	PI=4*ATAN(1.)
	684	DO N=0,NSTEPS
	685	IF ( N .EQ. 0 ) THEN
	686	W = 1.0
	687	ELSE
	688	W = ALPHA + (1-ALPHA)COS(NPI/(NSTEPS))
	689	ENDIF
	690	WW(N) = W
	691	ENDDO
	692
	693	RETURN
	694
	695	END SUBROUTINE hamming
	696
	697
	698	SUBROUTINE BLACKMAN(NSTEPS,WW)
	699
	700	! define blackman window function.
	701
	702	implicit none
	703
	704	integer, intent(in) :: nsteps
	705	real, dimension(0:nsteps) :: ww
	706	integer :: n
	707
	708	real :: pi, w
	709
	710	PI=4*ATAN(1.)
	711	DO N=0,NSTEPS
	712	IF ( N .EQ. 0 ) THEN
	713	W = 1.0
	714	ELSE
	715	W = 0.42 + 0.50COS( NPI/(NSTEPS)) &
	716	+ 0.08COS(2N*PI/(NSTEPS))
	717	ENDIF
	718	WW(N) = W
	719	ENDDO
	720
	721	RETURN
	722
	723	END SUBROUTINE blackman
	724
	725
	726	SUBROUTINE KAISER(NSTEPS,WW)
	727
	728	! define kaiser window function.
	729
	730	implicit none
	731
	732	real, external :: bessi0
	733
	734	integer, intent(in) :: nsteps
	735	real, dimension(0:nsteps) :: ww
	736	integer :: n
	737	real :: alpha, xi0a, xn, as
	738
	739	ALPHA = 1
	740
	741	XI0A = BESSI0(ALPHA)
	742	DO N=0,NSTEPS
	743	XN = N
	744	AS = ALPHASQRT(1.-(XN/NSTEPS)*2)
	745	WW(N) = BESSI0(AS) / XI0A
	746	ENDDO
	747
	748	RETURN
	749
	750	END SUBROUTINE kaiser
	751
	752
	753	REAL FUNCTION BESSI0(X)
	754
	755	! from numerical recipes (press, et al.)
	756
	757	implicit none
	758
	759	real(8) :: Y
	760	real(8) :: P1 = 1.0d0
	761	real(8) :: P2 = 3.5156230D0
	762	real(8) :: P3 = 3.0899424D0
	763	real(8) :: P4 = 1.2067492D0
	764	real(8) :: P5 = 0.2659732D0
	765	real(8) :: P6 = 0.360768D-1
	766	real(8) :: P7 = 0.45813D-2
	767
	768	real*8 :: Q1 = 0.39894228D0
	769	real*8 :: Q2 = 0.1328592D-1
	770	real*8 :: Q3 = 0.225319D-2
	771	real*8 :: Q4 = -0.157565D-2
	772	real*8 :: Q5 = 0.916281D-2
	773	real*8 :: Q6 = -0.2057706D-1
	774	real*8 :: Q7 = 0.2635537D-1
	775	real*8 :: Q8 = -0.1647633D-1
	776	real*8 :: Q9 = 0.392377D-2
	777
	778	real :: x, ax
	779
	780
	781	IF (ABS(X).LT.3.75) THEN
	782	Y=(X/3.75)**2
	783	BESSI0=P1+Y(P2+Y(P3+Y(P4+Y(P5+Y(P6+YP7)))))
	784	ELSE
	785	AX=ABS(X)
	786	Y=3.75/AX
	787	BESSI0=(EXP(AX)/SQRT(AX))(Q1+Y(Q2+Y(Q3+Y(Q4 &
	788	+Y(Q5+Y(Q6+Y(Q7+Y(Q8+Y*Q9))))))))
	789	ENDIF
	790	RETURN
	791
	792	END FUNCTION bessi0
	793
	794
	795	SUBROUTINE POTTER2(NSTEPS,WW)
	796
	797	! define potter window function.
	798	! modified to fall off over twice the range.
	799
	800	implicit none
	801
	802	integer, intent(in) :: nsteps
	803	real, dimension(0:nsteps),intent(out) :: ww
	804	integer :: n
	805	real :: ck, sum, arg
	806
	807	! local data
	808	real, dimension(0:3) :: d
	809	real :: pi
	810	integer :: ip
	811
	812	d(0) = 0.35577019
	813	d(1) = 0.2436983
	814	d(2) = 0.07211497
	815	d(3) = 0.00630165
	816
	817	PI=4*ATAN(1.)
	818
	819	CK = 1.0
	820	DO N=0,NSTEPS
	821	IF (N.EQ.NSTEPS) CK = 0.5
	822	ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
	823	!min--- modification in next statement
	824	ARG = ARG/2.
	825	!min--- end of modification
	826	SUM = D(0)
	827	DO IP=1,3
	828	SUM = SUM + 2.D(IP)COS(ARG*FLOAT(IP))
	829	END DO
	830	WW(N) = CK*SUM
	831	END DO
	832
	833	RETURN
	834
	835	END SUBROUTINE potter2
	836
	837
	838	SUBROUTINE dolphwin(m, window)
	839
	840	! calculation of dolph-chebyshev window or, for short,
	841	! dolph window, using the expression in the reference:
	842	!
	843	! antoniou, andreas, 1993: digital filters: analysis,
	844	! design and applications. mcgraw-hill, inc., 689pp.
	845	!
	846	! the dolph window is optimal in the following sense:
	847	! for a given main-lobe width, the stop-band attenuation
	848	! is minimal; for a given stop-band level, the main-lobe
	849	! width is minimal.
	850	!
	851	! it is possible to specify either the ripple-ratio r
	852	! or the stop-band edge thetas.
	853
	854	IMPLICIT NONE
	855
	856	! Arguments
	857	INTEGER, INTENT(IN) :: m
	858	REAL, DIMENSION(0:M), INTENT(OUT) :: window
	859
	860	! local data
	861	REAL, DIMENSION(0:2*M) :: t
	862	REAL, DIMENSION(0:M) :: w, time
	863	REAL :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
	864	INTEGER :: n, nm1, nt, i
	865
	866	PI = 4*ATAN(1.D0)
	867	THETAS = 2*PI/M
	868
	869	N = 2*M+1
	870	NM1 = N-1
	871	X0 = 1/COS(THETAS/2)
	872
	873	TERM1 = (X0 + SQRT(X02-1))(FLOAT(N-1))
	874	TERM2 = (X0 - SQRT(X02-1))(FLOAT(N-1))
	875	RR = 0.5*(TERM1+TERM2)
	876	R = 1/RR
	877	DB = 20*LOG10(R)
	878	WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
	879	WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
	880	WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
	881
	882	DO NT=0,M
	883	SUM = RR
	884	DO I=1,M
	885	ARG = X0COS(IPI/N)
	886	CALL CHEBY(T,NM1,ARG)
	887	TERM1 = T(NM1)
	888	TERM2 = COS(2NTPI*I/N)
	889	SUM = SUM + 2TERM1TERM2
	890	ENDDO
	891	W(NT) = SUM/N
	892	TIME(NT) = NT
	893	ENDDO
	894
	895	! fill up the array for return
	896	DO NT=0,M
	897	WINDOW(NT) = W(NT)
	898	ENDDO
	899
	900	RETURN
	901
	902	END SUBROUTINE dolphwin
	903
	904
	905	SUBROUTINE dolph(deltat, taus, m, window)
	906
	907	! calculation of dolph-chebyshev window or, for short,
	908	! dolph window, using the expression in the reference:
	909	!
	910	! antoniou, andreas, 1993: digital filters: analysis,
	911	! design and applications. mcgraw-hill, inc., 689pp.
	912	!
	913	! the dolph window is optimal in the following sense:
	914	! for a given main-lobe width, the stop-band attenuation
	915	! is minimal; for a given stop-band level, the main-lobe
	916	! width is minimal.
	917
	918	IMPLICIT NONE
	919
	920	! Arguments
	921	INTEGER, INTENT(IN) :: m
	922	REAL, DIMENSION(0:M), INTENT(OUT) :: window
	923	REAL, INTENT(IN) :: deltat, taus
	924
	925	! local data
	926	integer, PARAMETER :: NMAX = 5000
	927	REAL, dimension(0:NMAX) :: t, w, time
	928	real, dimension(0:2*nmax) :: w2
	929	INTEGER :: NPRPE=0 ! no of pe
	930	CHARACTER*80 :: MES
	931
	932	real :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
	933	integer :: n, nm1, i, nt
	934
	935	PI = 4*ATAN(1.D0)
	936
	937	print *, 'in dfcoef, deltat = ', deltat, 'taus=',taus
	938
	939	N = 2*M+1
	940	NM1 = N-1
	941
	942	THETAS = 2PIABS(DELTAT/TAUS)
	943	X0 = 1/COS(THETAS/2)
	944	TERM1 = (X0 + SQRT(X02-1))(FLOAT(N-1))
	945	TERM2 = (X0 - SQRT(X02-1))(FLOAT(N-1))
	946	RR = 0.5*(TERM1+TERM2)
	947	R = 1/RR
	948	DB = 20*LOG10(R)
	949
	950
	951	WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
	952	WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
	953	WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
	954
	955	DO NT=0,M
	956	SUM = 1
	957	DO I=1,M
	958	ARG = X0COS(IPI/N)
	959	CALL CHEBY(T,NM1,ARG)
	960	TERM1 = T(NM1)
	961	TERM2 = COS(2NTPI*I/N)
	962	SUM = SUM + R2TERM1*TERM2
	963	ENDDO
	964	W(NT) = SUM/N
	965	TIME(NT) = NT
	966	WRITE(*,'(1X,''DOLPH: TIME, W='',F10.6,2X,E17.7)') &
	967	TIME(NT), W(NT)
	968	ENDDO
	969	! fill in the negative-time values by symmetry.
	970	DO NT=0,M
	971	W2(M+NT) = W(NT)
	972	W2(M-NT) = W(NT)
	973	ENDDO
	974
	975	! fill up the array for return
	976	SUMW = 0.
	977	DO NT=0,2*M
	978	SUMW = SUMW + W2(NT)
	979	ENDDO
	980	WRITE(*,'(1X,''DOLPH: SUM OF WEIGHTS W2='',F10.4)')SUMW
	981
	982	DO NT=0,2*M
	983	WINDOW(NT) = W2(NT)
	984	ENDDO
	985
	986	RETURN
	987
	988	END SUBROUTINE dolph
	989
	990
	991	SUBROUTINE cheby(t, n, x)
	992
	993	! calculate all chebyshev polynomials up to order n
	994	! for the argument value x.
	995
	996	! reference: numerical recipes, page 184, recurrence
	997	! t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) , n>=2.
	998
	999	IMPLICIT NONE
	1000
	1001	! Arguments
	1002	INTEGER, INTENT(IN) :: n
	1003	REAL, INTENT(IN) :: x
	1004	REAL, DIMENSION(0:N) :: t
	1005
	1006	integer :: nn
	1007
	1008	T(0) = 1
	1009	T(1) = X
	1010	IF(N.LT.2) RETURN
	1011	DO NN=2,N
	1012	T(NN) = 2XT(NN-1) - T(NN-2)
	1013	ENDDO
	1014
	1015	RETURN
	1016
	1017	END SUBROUTINE cheby
	1018
	1019
	1020	SUBROUTINE rhofil(dt, tauc, norder, nstep, ictype, frow, nosize)
	1021
	1022	! RHO = recurssive high order.
	1023	!
	1024	! This routine calculates and returns the
	1025	! Last Row, FROW, of the FILTER matrix.
	1026	!
	1027	! Input Parameters:
	1028	! DT : Time Step in seconds
	1029	! TAUC : Cut-off period (hours)
	1030	! NORDER : Order of QS Filter
	1031	! NSTEP : Number of step/Size of row.
	1032	! ICTYPE : Initial Conditions
	1033	! NOSIZE : Max. side of FROW.
	1034	!
	1035	! Working Fields:
	1036	! ACOEF : X-coefficients of filter
	1037	! BCOEF : Y-coefficients of filter
	1038	! FILTER : Filter Matrix.
	1039	!
	1040	! Output Parameters:
	1041	! FROW : Last Row of Filter Matrix.
	1042	!
	1043	! Note: Two types of initial conditions are permitted.
	1044	! ICTYPE = 1 : Order increasing each row to NORDER.
	1045	! ICTYPE = 2 : Order fixed at NORDER throughout.
	1046	!
	1047	! DOUBLE PRECISION USED THROUGHOUT.
	1048
	1049	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	1050
	1051	DOUBLE PRECISION MUC
	1052
	1053	! N.B. Single Precision for List Parameters.
	1054	REAL, intent(in) :: DT,TAUC
	1055
	1056	! Space for the last row of FILTER.
	1057	integer, intent(in) :: norder, nstep, ictype, nosize
	1058	REAL , dimension(0:nosize), intent(out):: FROW
	1059
	1060	! Arrays for rho filter.
	1061	integer, PARAMETER :: NOMAX=100
	1062	real , dimension(0:NOMAX) :: acoef, bcoef
	1063	real , dimension(0:NOMAX,0:NOMAX) :: filter
	1064	! Working space.
	1065	real , dimension(0:NOMAX) :: alpha, beta
	1066
	1067	real :: DTT
	1068
	1069	DTT = ABS(DT)
	1070	PI = 2*DASIN(1.D0)
	1071	IOTA = CMPLX(0.,1.)
	1072
	1073	! Filtering Parameters (derived).
	1074	THETAC = 2PIDTT/(TAUC)
	1075	MUC = tan(THETAC/2)
	1076	FC = THETAC/(2*PI)
	1077
	1078	! Clear the arrays.
	1079	DO NC=0,NOMAX
	1080	ACOEF(NC) = 0.
	1081	BCOEF(NC) = 0.
	1082	ALPHA(NC) = 0.
	1083	BETA (NC) = 0.
	1084	FROW (NC) = 0.
	1085	DO NR=0,NOMAX
	1086	FILTER(NR,NC) = 0.
	1087	ENDDO
	1088	ENDDO
	1089
	1090	! Fill up the Filter Matrix.
	1091	FILTER(0,0) = 1.
	1092
	1093	! Get the coefficients of the Filter.
	1094	IF ( ICTYPE.eq.2 ) THEN
	1095	CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
	1096	ENDIF
	1097
	1098	DO 100 NROW=1,NSTEP
	1099
	1100	IF ( ICTYPE.eq.1 ) THEN
	1101	NORD = MIN(NROW,NORDER)
	1102	IF ( NORD.le.NORDER) THEN
	1103	CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
	1104	ENDIF
	1105	ENDIF
	1106
	1107	DO K=0,NROW
	1108	ALPHA(K) = ACOEF(NROW-K)
	1109	IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
	1110	ENDDO
	1111
	1112	! Correction for terms of negative index.
	1113	IF ( ICTYPE.eq.2 ) THEN
	1114	IF ( NROW.lt.NORDER ) THEN
	1115	CN = 0.
	1116	DO NN=NROW+1,NORDER
	1117	CN = CN + (ACOEF(NN)+BCOEF(NN))
	1118	ENDDO
	1119	ALPHA(0) = ALPHA(0) + CN
	1120	ENDIF
	1121	ENDIF
	1122
	1123	! Check sum of ALPHAs and BETAs = 1
	1124	SUMAB = 0.
	1125	DO NN=0,NROW
	1126	SUMAB = SUMAB + ALPHA(NN)
	1127	IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
	1128	ENDDO
	1129
	1130	DO KK=0,NROW-1
	1131	SUMBF = 0.
	1132	DO LL=0,NROW-1
	1133	SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
	1134	ENDDO
	1135	FILTER(NROW,KK) = ALPHA(KK)+SUMBF
	1136	ENDDO
	1137	FILTER(NROW,NROW) = ALPHA(NROW)
	1138
	1139	! Check sum of row elements = 1
	1140	SUMROW = 0.
	1141	DO NN=0,NROW
	1142	SUMROW = SUMROW + FILTER(NROW,NN)
	1143	ENDDO
	1144
	1145	100 CONTINUE
	1146
	1147	DO NC=0,NSTEP
	1148	FROW(NC) = FILTER(NSTEP,NC)
	1149	ENDDO
	1150
	1151	RETURN
	1152
	1153	END SUBROUTINE rhofil
	1154
	1155
	1156	SUBROUTINE rhocof(nord, nomax, muc, ca, cb)
	1157
	1158	! Get the coefficients of the RHO Filter
	1159
	1160	! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	1161	IMPLICIT NONE
	1162
	1163	! Arguments
	1164	integer, intent(in) :: nord, nomax
	1165	real, dimension(0:nomax) :: ca, cb
	1166
	1167	! Functions
	1168	double precision, external :: cnr
	1169
	1170	! Local variables
	1171	INTEGER :: nn
	1172	COMPLEX :: IOTA
	1173	DOUBLE PRECISION :: MUC, ZN
	1174	DOUBLE PRECISION :: pi, root2, rn, sigma, gain, sumcof
	1175
	1176	PI = 2*ASIN(1.)
	1177	ROOT2 = SQRT(2.)
	1178	IOTA = (0.,1.)
	1179
	1180	RN = 1./FLOAT(NORD)
	1181	SIGMA = 1./( SQRT(2.**RN-1.) )
	1182
	1183	GAIN = (MUCSIGMA/(1+MUCSIGMA))**NORD
	1184	ZN = (1-MUCSIGMA)/(1+MUCSIGMA)
	1185
	1186	DO NN=0,NORD
	1187	CA(NN) = CNR(NORD,NN)*GAIN
	1188	IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)(-ZN)*NN
	1189	ENDDO
	1190
	1191	! Check sum of coefficients = 1
	1192	SUMCOF = 0.
	1193	DO NN=0,NORD
	1194	SUMCOF = SUMCOF + CA(NN)
	1195	IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
	1196	ENDDO
	1197
	1198	RETURN
	1199
	1200	END SUBROUTINE RHOCOF
	1201
	1202
	1203	DOUBLE PRECISION FUNCTION cnr(n,r)
	1204
	1205	! Binomial Coefficient (n,r).
	1206
	1207	! IMPLICIT DOUBLE PRECISION(C,X)
	1208	IMPLICIT NONE
	1209
	1210	! Arguments
	1211	INTEGER , intent(in) :: n, R
	1212
	1213	! Local variables
	1214	INTEGER :: k
	1215	DOUBLE PRECISION :: coeff, xn, xr, xk
	1216
	1217	IF ( R.eq.0 ) THEN
	1218	CNR = 1.0
	1219	RETURN
	1220	ENDIF
	1221	Coeff = 1.0
	1222	XN = DFLOAT(N)
	1223	XR = DFLOAT(R)
	1224	DO K=1,R
	1225	XK = DFLOAT(K)
	1226	COEFF = COEFF * ( (XN-XR+XK)/XK )
	1227	ENDDO
	1228	CNR = COEFF
	1229
	1230	RETURN
	1231
	1232	END FUNCTION cnr
	1233
	1234
	1235	SUBROUTINE optfil (grid,NH,DELTAT,NHMAX)
	1236	!----------------------------------------------------------------------
	1237
	1238	! SUBROUTINE optfil (NH,DELTAT,TAUP,TAUS,LPRINT, &
	1239	! H,NHMAX)
	1240	!
	1241	! - Huang and Lynch optimal filter
	1242	! Monthly Weather Review, Feb 1993
	1243	!----------------------------------------------------------
	1244	! Input Parameters in List:
	1245	! NH : Half-length of the Filter
	1246	! DELTAT : Time-step (in seconds).
	1247	! TAUP : Period of pass-band edge (hours).
	1248	! TAUS : Period of stop-band edge (hours).
	1249	! LPRINT : Logical switch for messages.
	1250	! NHMAX : Maximum permitted Half-length.
	1251	!
	1252	! Output Parameters in List:
	1253	! H : Impulse Response.
	1254	! DP : Deviation in pass-band (db)
	1255	! DS : Deviation in stop-band (db)
	1256	!----------------------------------------------------------
	1257	!
	1258	USE module_domain
	1259
	1260	TYPE(domain) , POINTER :: grid
	1261
	1262	REAL,DIMENSION( 20) :: EDGE
	1263	REAL,DIMENSION( 10) :: FX, WTX, DEVIAT
	1264	REAL,DIMENSION(2*NHMAX+1) :: H
	1265	logical LPRINT
	1266	REAL, INTENT (IN) :: DELTAT
	1267	INTEGER, INTENT (IN) :: NH, NHMAX
	1268	!
	1269	TAUP = 3.
	1270	TAUS = 1.5
	1271	LPRINT = .true.
	1272	!initialize H array
	1273
	1274	NL=2*NHMAX+1
	1275	do 101 n=1,NL
	1276	H(n)=0.
	1277	101 continue
	1278
	1279	NFILT = 2*NH+1
	1280	print *,' start optfil, NFILT=', nfilt
	1281
	1282	!
	1283	! 930325 PL & XYH : the upper limit is changed from 64 to 128.
	1284	IF(NFILT.LE.0 .OR. NFILT.GT.128 ) THEN
	1285	WRITE(6,*) 'NH=',NH
	1286	CALL wrf_error_fatal (' Sorry, error 1 in call to OPTFIL ')
	1287	ENDIF
	1288	!
	1289	! The following four should always be the same.
	1290	JTYPE = 1
	1291	NBANDS = 2
	1292	!CC JPRINT = 0
	1293	LGRID = 16
	1294	!
	1295	! calculate transition frequencies.
	1296	DT = ABS(DELTAT)
	1297	FS = DT/(TAUS*3600.)
	1298	FP = DT/(TAUP*3600.)
	1299	IF(FS.GT.0.5) then
	1300	! print *,' FS too large in OPTFIL '
	1301	CALL wrf_error_fatal (' FS too large in OPTFIL ')
	1302	! return
	1303	end if
	1304	IF(FP.LT.0.0) then
	1305	! print *, ' FP too small in OPTFIL '
	1306	CALL wrf_error_fatal (' FP too small in OPTFIL ')
	1307	! return
	1308	end if
	1309	!
	1310	! Relative Weights in pass- and stop-bands.
	1311	WTP = 1.0
	1312	WTS = 1.0
	1313	!
	1314	!CC NOTE: (FP,FC,FS) is an arithmetic progression, so
	1315	!CC (1/FS,1/FC,1/FP) is a harmonic one.
	1316	!CC TAUP = 1/( (1/TAUC)-(1/DTAU) )
	1317	!CC TAUS = 1/( (1/TAUC)+(1/DTAU) )
	1318	!CC TAUC : Cut-off Period (hours).
	1319	!CC DTAU : Transition half-width (hours).
	1320	!CC FC = 1/TAUC ; DF = 1/DTAU
	1321	!CC FP = FC - DF : FS = FC + DF
	1322	!
	1323	IF ( LPRINT ) THEN
	1324	TAUC = 2./((1/TAUS)+(1/TAUP))
	1325	DTAU = 2./((1/TAUS)-(1/TAUP))
	1326	FC = DT/(TAUC*3600.)
	1327	DF = DT/(DTAU*3600.)
	1328	WRITE(6,*) ' DT ' , dt
	1329	WRITE(6,*) ' TAUS, TAUP ' , TAUS,TAUP
	1330	WRITE(6,*) ' TAUC, DTAU ' , TAUC,DTAU
	1331	WRITE(6,*) ' FP, FS ' , FP, FS
	1332	WRITE(6,*) ' FC, DF ' , FC, DF
	1333	WRITE(6,*) ' WTS, WTP ' , WTS, WTP
	1334	ENDIF
	1335	!
	1336	! Fill the control vectors for MCCPAR
	1337	EDGE(1) = 0.0
	1338	EDGE(2) = FP
	1339	EDGE(3) = FS
	1340	EDGE(4) = 0.5
	1341	FX(1) = 1.0
	1342	FX(2) = 0.0
	1343	WTX(1) = WTP
	1344	WTX(2) = WTS
	1345
	1346	CALL MCCPAR(NFILT,JTYPE,NBANDS,LPRINT,LGRID, &
	1347	EDGE,FX,WTX,DEVIAT, h )
	1348	!
	1349	! Save the deviations in the pass- and stop-bands.
	1350	DP = DEVIAT(1)
	1351	DS = DEVIAT(2)
	1352	!
	1353	! Fill out the array H (only first half filled in MCCPAR).
	1354	IF(MOD(NFILT,2).EQ.0) THEN
	1355	NHALF = ( NFILT )/2
	1356	ELSE
	1357	NHALF = (NFILT+1)/2
	1358	ENDIF
	1359	DO 100 nn=1,NHALF
	1360	H(NFILT+1-nn) = h(nn)
	1361	100 CONTINUE
	1362
	1363	! normalize the sums to be unity
	1364	sumh = 0
	1365	do 150 n=1,NFILT
	1366	sumh = sumh + H(n)
	1367	150 continue
	1368	print *,'SUMH =', sumh
	1369
	1370	do 200 n=1,NFILT
	1371	H(n) = H(n)/sumh
	1372	200 continue
	1373	do 201 n=1,NFILT
	1374	grid%hcoeff(n)=H(n)
	1375	201 continue
	1376	! print *,'HCOEFF(n) ', grid%hcoeff
	1377	!
	1378	END SUBROUTINE optfil
	1379
	1380
	1381	SUBROUTINE MCCPAR (NFILT,JTYPE,NBANDS,LPRINT,LGRID, &
	1382	EDGE,FX,WTX,DEVIAT,h )
	1383
	1384	! PROGRAM FOR THE DESIGN OF LINEAR PHASE FINITE IMPULSE
	1385	! REPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM
	1386	!
	1387	!************************************************************
	1388	!* Reference: McClellan, J.H., T.W. Parks and L.R.Rabiner, *
	1389	!* 1973: A computer program for designing *
	1390	!* optimum FIR linear phase digital filters. *
	1391	!* IEEE Trans. on Audio and Electroacoustics, *
	1392	!* Vol AU-21, No. 6, 506-526. *
	1393	!************************************************************
	1394	!
	1395	! THREE TYPES OF FILTERS ARE INCLUDED -- BANDPASS FILTERS
	1396	! DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS
	1397	!
	1398	!---------------------------------------------------------------
	1399	!
	1400	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
	1401	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
	1402	DIMENSION H(66)
	1403	DIMENSION DES(1045),GRID(1045),WT(1045)
	1404	DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
	1405	DOUBLE PRECISION PI2,PI
	1406	DOUBLE PRECISION AD,DEV,X,Y
	1407	LOGICAL LPRINT
	1408
	1409	PI = 3.141592653589793
	1410	PI2 = 6.283185307179586
	1411
	1412	! ......
	1413
	1414	NFMAX = 128
	1415	100 CONTINUE
	1416
	1417	! PROGRAM INPUT SECTION
	1418
	1419	!CC READ(5,*) NFILT,JTYPE,NBANDS,JPRINT,LGRID
	1420
	1421	IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) THEN
	1422	CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
	1423	END IF
	1424	IF(NBANDS.LE.0) NBANDS = 1
	1425
	1426	! ....
	1427
	1428	IF(LGRID.LE.0) LGRID = 16
	1429	JB = 2*NBANDS
	1430	!cc READ(5,*) (EDGE(J),J=1,JB)
	1431	!cc READ(5,*) (FX(J),J=1,NBANDS)
	1432	!cc READ(5,*) (WTX(J),J=1,NBANDS)
	1433	IF(JTYPE.EQ.0) THEN
	1434	CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
	1435	END IF
	1436	NEG = 1
	1437	IF(JTYPE.EQ.1) NEG = 0
	1438	NODD = NFILT/2
	1439	NODD = NFILT-2*NODD
	1440	NFCNS = NFILT/2
	1441	IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS = NFCNS+1
	1442
	1443	! ...
	1444
	1445	GRID(1) = EDGE(1)
	1446	DELF = LGRID*NFCNS
	1447	DELF = 0.5/DELF
	1448	IF(NEG.EQ.0) GOTO 135
	1449	IF(EDGE(1).LT.DELF) GRID(1) = DELF
	1450	135 CONTINUE
	1451	J = 1
	1452	L = 1
	1453	LBAND = 1
	1454	140 FUP = EDGE(L+1)
	1455	145 TEMP = GRID(J)
	1456
	1457	! ....
	1458
	1459	DES(J) = EFF(TEMP,FX,WTX,LBAND,JTYPE)
	1460	WT(J) = WATE(TEMP,FX,WTX,LBAND,JTYPE)
	1461	J = J+1
	1462	GRID(J) = TEMP+DELF
	1463	IF(GRID(J).GT.FUP) GOTO 150
	1464	GOTO 145
	1465	150 GRID(J-1) = FUP
	1466	DES(J-1) = EFF(FUP,FX,WTX,LBAND,JTYPE)
	1467	WT(J-1) = WATE(FUP,FX,WTX,LBAND,JTYPE)
	1468	LBAND = LBAND+1
	1469	L = L+2
	1470	IF(LBAND.GT.NBANDS) GOTO 160
	1471	GRID(J) = EDGE(L)
	1472	GOTO 140
	1473	160 NGRID = J-1
	1474	IF(NEG.NE.NODD) GOTO 165
	1475	IF(GRID(NGRID).GT.(0.5-DELF)) NGRID = NGRID-1
	1476	165 CONTINUE
	1477
	1478	! ......
	1479
	1480	IF(NEG) 170,170,180
	1481	170 IF(NODD.EQ.1) GOTO 200
	1482	DO 175 J=1,NGRID
	1483	CHANGE = DCOS(PI*GRID(J))
	1484	DES(J) = DES(J)/CHANGE
	1485	WT(J) = WT(J)*CHANGE
	1486	175 CONTINUE
	1487	GOTO 200
	1488	180 IF(NODD.EQ.1) GOTO 190
	1489	DO 185 J = 1,NGRID
	1490	CHANGE = DSIN(PI*GRID(J))
	1491	DES(J) = DES(J)/CHANGE
	1492	WT(J) = WT(J)*CHANGE
	1493	185 CONTINUE
	1494	GOTO 200
	1495	190 DO 195 J =1,NGRID
	1496	CHANGE = DSIN(PI2*GRID(J))
	1497	DES(J) = DES(J)/CHANGE
	1498	WT(J) = WT(J)*CHANGE
	1499	195 CONTINUE
	1500
	1501	! ......
	1502
	1503	200 TEMP = FLOAT(NGRID-1)/FLOAT(NFCNS)
	1504	DO 210 J = 1,NFCNS
	1505	IEXT(J) = (J-1)*TEMP+1
	1506	210 CONTINUE
	1507	IEXT(NFCNS+1) = NGRID
	1508	NM1 = NFCNS-1
	1509	NZ = NFCNS+1
	1510
	1511	! CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM
	1512
	1513	CALL REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
	1514
	1515	! CALCULATE THE IMPULSE RESPONSE
	1516
	1517	IF(NEG) 300,300,320
	1518	300 IF(NODD.EQ.0) GOTO 310
	1519	DO 305 J=1,NM1
	1520	H(J) = 0.5*ALPHA(NZ-J)
	1521	305 CONTINUE
	1522	H(NFCNS)=ALPHA(1)
	1523	GOTO 350
	1524	310 H(1) = 0.25*ALPHA(NFCNS)
	1525	DO 315 J = 2,NM1
	1526	H(J) = 0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
	1527	315 CONTINUE
	1528	H(NFCNS) = 0.5ALPHA(1)+0.25ALPHA(2)
	1529	GOTO 350
	1530	320 IF(NODD.EQ.0) GOTO 330
	1531	H(1) = 0.25*ALPHA(NFCNS)
	1532	H(2) = 0.25*ALPHA(NM1)
	1533	DO 325 J = 3,NM1
	1534	H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
	1535	325 CONTINUE
	1536	H(NFCNS) = 0.5ALPHA(1)-0.25ALPHA(3)
	1537	H(NZ) = 0.0
	1538	GOTO 350
	1539	330 H(1) = 0.25*ALPHA(NFCNS)
	1540	DO 335 J =2,NM1
	1541	H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
	1542	335 CONTINUE
	1543	H(NFCNS) = 0.5ALPHA(1)-0.25ALPHA(2)
	1544
	1545	! PROGRAM OUTPUT SECTION
	1546
	1547	350 CONTINUE
	1548	!
	1549	IF(LPRINT) THEN
	1550
	1551	print , '***************************************************'
	1552	print *, 'FINITE IMPULSE RESPONSE (FIR)'
	1553	print *, 'LINEAR PHASE DIGITAL FILTER DESIGN'
	1554	print *, 'REMEZ EXCHANGE ALGORITHM'
	1555
	1556	IF(JTYPE.EQ.1) WRITE(6,365)
	1557	365 FORMAT(25X,'BANDPASS FILTER'/)
	1558
	1559	IF(JTYPE.EQ.2) WRITE(6,370)
	1560	370 FORMAT(25X,'DIFFERENTIATOR '/)
	1561
	1562	IF(JTYPE.EQ.3) WRITE(6,375)
	1563	375 FORMAT(25X,'HILBERT TRANSFORMER '/)
	1564
	1565	WRITE(6,378) NFILT
	1566	378 FORMAT(15X,'FILTER LENGTH =',I3/)
	1567
	1568	WRITE(6,380)
	1569	380 FORMAT(15X,'*** IMPULSE RESPONSE ***')
	1570
	1571	DO 381 J = 1,NFCNS
	1572	K = NFILT+1-J
	1573	IF(NEG.EQ.0) WRITE(6,382) J,H(J),K
	1574	IF(NEG.EQ.1) WRITE(6,383) J,H(J),K
	1575	381 CONTINUE
	1576	382 FORMAT(20X,'H(',I3,') = ',E15.8,' = H(',I4,')')
	1577	383 FORMAT(20X,'H(',I3,') = ',E15.8,' = -H(',I4,')')
	1578
	1579	IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ
	1580	384 FORMAT(20X,'H(',I3,') = 0.0')
	1581
	1582	DO 450 K=1,NBANDS,4
	1583	KUP = K+3
	1584	IF(KUP.GT.NBANDS) KUP = NBANDS
	1585	print *
	1586	WRITE(6,385) (J,J=K,KUP)
	1587	385 FORMAT(24X,4('BAND',I3,8X))
	1588	WRITE(6,390) (EDGE(2*J-1),J=K,KUP)
	1589	390 FORMAT(2X,'LOWER BAND EDGE',5F15.8)
	1590	WRITE(6,395) (EDGE(2*J),J=K,KUP)
	1591	395 FORMAT(2X,'UPPER BAND EDGE',5F15.8)
	1592	IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
	1593	400 FORMAT(2X,'DESIRED VALUE',2X,5F15.8)
	1594	IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
	1595	405 FORMAT(2X,'DESIRED SLOPE',2X,5F15.8)
	1596	WRITE(6,410) (WTX(J),J=K,KUP)
	1597	410 FORMAT(2X,'WEIGHTING',6X,5F15.8)
	1598	DO 420 J = K,KUP
	1599	DEVIAT(J) = DEV/WTX(J)
	1600	420 CONTINUE
	1601	WRITE(6,425) (DEVIAT(J),J=K,KUP)
	1602	425 FORMAT(2X,'DEVIATION',6X,5F15.8)
	1603	IF(JTYPE.NE.1) GOTO 450
	1604	DO 430 J = K,KUP
	1605	DEVIAT(J) = 20.0*ALOG10(DEVIAT(J))
	1606	430 CONTINUE
	1607	WRITE(6,435) (DEVIAT(J),J=K,KUP)
	1608	435 FORMAT(2X,'DEVIATION IN DB',5F15.8)
	1609	450 CONTINUE
	1610	print *, 'EXTREMAL FREQUENCIES'
	1611	WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
	1612	455 FORMAT((2X,5F15.7))
	1613	WRITE(6,460)
	1614	460 FORMAT(1X,70(1H*))
	1615	!
	1616	ENDIF
	1617	!
	1618	!CC IF(NFILT.NE.0) GOTO 100 ! removal of re-run loop.
	1619	!
	1620	END SUBROUTINE mccpar
	1621
	1622
	1623	FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE)
	1624	DIMENSION FX(5),WTX(5)
	1625	IF(JTYPE.EQ.2) GOTO 1
	1626	EFF = FX(LBAND)
	1627	RETURN
	1628	1 EFF = FX(LBAND)*TEMP
	1629	END FUNCTION eff
	1630
	1631
	1632	FUNCTION WATE(TEMP,FX,WTX,LBAND,JTYPE)
	1633	DIMENSION FX(5),WTX(5)
	1634	IF(JTYPE.EQ.2) GOTO 1
	1635	WATE = WTX(LBAND)
	1636	RETURN
	1637	1 IF(FX(LBAND).LT.0.0001) GOTO 2
	1638	WATE = WTX(LBAND)/TEMP
	1639	RETURN
	1640	2 WATE = WTX(LBAND)
	1641	END FUNCTION wate
	1642
	1643
	1644	! SUBROUTINE ERROR
	1645	!! WRITE(6,)' * ERROR IN INPUT DATA **'
	1646	! CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
	1647	! END SUBROUTINE error
	1648
	1649
	1650	SUBROUTINE REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
	1651	! THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
	1652	! FOR THE WEIGHTED CHEBCHEV APPROXIMATION OF A CONTINUOUS
	1653	! FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE
	1654	! ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
	1655	! DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
	1656	! GRID, THE NUMBER OF COSINES, AND THE INITIAL GUESS OF THE
	1657	! EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV
	1658	! ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
	1659	! FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
	1660	! THE COEFFICIENTS OF THE BEST APPROXIMATION.
	1661	!
	1662	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
	1663	DIMENSION EDGE(20)
	1664	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
	1665	DIMENSION DES(1045),GRID(1045),WT(1045)
	1666	DIMENSION A(66),P(65),Q(65)
	1667	DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
	1668	DOUBLE PRECISION AD,DEV,X,Y
	1669	DOUBLE PRECISION, EXTERNAL :: D, GEE
	1670	!
	1671	! THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
	1672	!
	1673	ITRMAX=25
	1674	DEVL=-1.0
	1675	NZ=NFCNS+1
	1676	NZZ=NFCNS+2
	1677	NITER=0
	1678	100 CONTINUE
	1679	IEXT(NZZ)=NGRID+1
	1680	NITER=NITER+1
	1681	IF(NITER.GT.ITRMAX) GO TO 400
	1682	DO 110 J=1,NZ
	1683	DTEMP=GRID(IEXT(J))
	1684	DTEMP=DCOS(DTEMP*PI2)
	1685	110 X(J)=DTEMP
	1686	JET=(NFCNS-1)/15+1
	1687	DO 120 J=1,NZ
	1688	120 AD(J)=D(J,NZ,JET,X)
	1689	DNUM=0.0
	1690	DDEN=0.0
	1691	K=1
	1692	DO 130 J=1,NZ
	1693	L=IEXT(J)
	1694	DTEMP=AD(J)*DES(L)
	1695	DNUM=DNUM+DTEMP
	1696	DTEMP=K*AD(J)/WT(L)
	1697	DDEN=DDEN+DTEMP
	1698	130 K=-K
	1699	DEV=DNUM/DDEN
	1700	NU=1
	1701	IF(DEV.GT.0.0) NU=-1
	1702	DEV=-NU*DEV
	1703	K=NU
	1704	DO 140 J=1,NZ
	1705	L=IEXT(J)
	1706	DTEMP=K*DEV/WT(L)
	1707	Y(J)=DES(L)+DTEMP
	1708	140 K=-K
	1709	IF(DEV.GE.DEVL) GO TO 150
	1710	WRITE(6,) ' ***** FAILURE TO CONVERGE ********* '
	1711	WRITE(6,*) ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR '
	1712	WRITE(6,*) ' THE IMPULSE RESPONSE MAY BE CORRECT '
	1713	WRITE(6,*) ' CHECK WITH A FREQUENCY RESPONSE '
	1714	WRITE(6,) ' *************************************** '
	1715	GO TO 400
	1716	150 DEVL=DEV
	1717	JCHNGE=0
	1718	K1=IEXT(1)
	1719	KNZ=IEXT(NZ)
	1720	KLOW=0
	1721	NUT=-NU
	1722	J=1
	1723	!
	1724	! SEARCH FOR THE EXTERMAL FREQUENCIES OF THE BEST
	1725	! APPROXIMATION.
	1726
	1727	200 IF(J.EQ.NZZ) YNZ=COMP
	1728	IF(J.GE.NZZ) GO TO 300
	1729	KUP=IEXT(J+1)
	1730	L=IEXT(J)+1
	1731	NUT=-NUT
	1732	IF(J.EQ.2) Y1=COMP
	1733	COMP=DEV
	1734	IF(L.GE.KUP) GO TO 220
	1735	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1736	ERR=(ERR-DES(L))*WT(L)
	1737	DTEMP=NUT*ERR-COMP
	1738	IF(DTEMP.LE.0.0) GO TO 220
	1739	COMP=NUT*ERR
	1740	210 L=L+1
	1741	IF(L.GE.KUP) GO TO 215
	1742	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1743	ERR=(ERR-DES(L))*WT(L)
	1744	DTEMP=NUT*ERR-COMP
	1745	IF(DTEMP.LE.0.0) GO TO 215
	1746	COMP=NUT*ERR
	1747	GO TO 210
	1748	215 IEXT(J)=L-1
	1749	J=J+1
	1750	KLOW=L-1
	1751	JCHNGE=JCHNGE+1
	1752	GO TO 200
	1753	220 L=L-1
	1754	225 L=L-1
	1755	IF(L.LE.KLOW) GO TO 250
	1756	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1757	ERR=(ERR-DES(L))*WT(L)
	1758	DTEMP=NUT*ERR-COMP
	1759	IF(DTEMP.GT.0.0) GO TO 230
	1760	IF(JCHNGE.LE.0) GO TO 225
	1761	GO TO 260
	1762	230 COMP=NUT*ERR
	1763	235 L=L-1
	1764	IF(L.LE.KLOW) GO TO 240
	1765	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1766	ERR=(ERR-DES(L))*WT(L)
	1767	DTEMP=NUT*ERR-COMP
	1768	IF(DTEMP.LE.0.0) GO TO 240
	1769	COMP=NUT*ERR
	1770	GO TO 235
	1771	240 KLOW=IEXT(J)
	1772	IEXT(J)=L+1
	1773	J=J+1
	1774	JCHNGE=JCHNGE+1
	1775	GO TO 200
	1776	250 L=IEXT(J)+1
	1777	IF(JCHNGE.GT.0) GO TO 215
	1778	255 L=L+1
	1779	IF(L.GE.KUP) GO TO 260
	1780	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1781	ERR=(ERR-DES(L))*WT(L)
	1782	DTEMP=NUT*ERR-COMP
	1783	IF(DTEMP.LE.0.0) GO TO 255
	1784	COMP=NUT*ERR
	1785	GO TO 210
	1786	260 KLOW=IEXT(J)
	1787	J=J+1
	1788	GO TO 200
	1789	300 IF(J.GT.NZZ) GO TO 320
	1790	IF(K1.GT.IEXT(1)) K1=IEXT(1)
	1791	IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
	1792	NUT1=NUT
	1793	NUT=-NU
	1794	L=0
	1795	KUP=K1
	1796	COMP=YNZ*(1.00001)
	1797	LUCK=1
	1798	310 L=L+1
	1799	IF(L.GE.KUP) GO TO 315
	1800	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1801	ERR=(ERR-DES(L))*WT(L)
	1802	DTEMP=NUT*ERR-COMP
	1803	IF(DTEMP.LE.0.0) GO TO 310
	1804	COMP=NUT*ERR
	1805	J=NZZ
	1806	GO TO 210
	1807	315 LUCK=6
	1808	GO TO 325
	1809	320 IF(LUCK.GT.9) GO TO 350
	1810	IF(COMP.GT.Y1) Y1=COMP
	1811	K1=IEXT(NZZ)
	1812	325 L=NGRID+1
	1813	KLOW=KNZ
	1814	NUT=-NUT1
	1815	COMP=Y1*(1.00001)
	1816	330 L=L-1
	1817	IF(L.LE.KLOW) GO TO 340
	1818	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
	1819	ERR=(ERR-DES(L))*WT(L)
	1820	DTEMP=NUT*ERR-COMP
	1821	IF(DTEMP.LE.0.0) GO TO 330
	1822	J=NZZ
	1823	COMP=NUT*ERR
	1824	LUCK=LUCK+10
	1825	GO TO 235
	1826	340 IF(LUCK.EQ.6) GO TO 370
	1827	DO 345 J=1,NFCNS
	1828	345 IEXT(NZZ-J)=IEXT(NZ-J)
	1829	IEXT(1)=K1
	1830	GO TO 100
	1831	350 KN=IEXT(NZZ)
	1832	DO 360 J=1,NFCNS
	1833	360 IEXT(J)=IEXT(J+1)
	1834	IEXT(NZ)=KN
	1835	GO TO 100
	1836	370 IF(JCHNGE.GT.0) GO TO 100
	1837	!
	1838	! CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
	1839	! USING THE INVERSE DISCRETE FOURIER TRANSFORM.
	1840	!
	1841	400 CONTINUE
	1842	NM1=NFCNS-1
	1843	FSH=1.0E-06
	1844	GTEMP=GRID(1)
	1845	X(NZZ)=-2.0
	1846	CN=2*NFCNS-1
	1847	DELF=1.0/CN
	1848	L=1
	1849	KKK=0
	1850	IF(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1
	1851	IF(NFCNS.LE.3) KKK=1
	1852	IF(KKK.EQ.1) GO TO 405
	1853	DTEMP=DCOS(PI2*GRID(1))
	1854	DNUM=DCOS(PI2*GRID(NGRID))
	1855	AA=2.0/(DTEMP-DNUM)
	1856	BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
	1857	405 CONTINUE
	1858	DO 430 J=1,NFCNS
	1859	FT=(J-1)*DELF
	1860	XT=DCOS(PI2*FT)
	1861	IF(KKK.EQ.1) GO TO 410
	1862	XT=(XT-BB)/AA
	1863	! original : FT=ARCOS(XT)/PI2
	1864	FT=ACOS(XT)/PI2
	1865	410 XE=X(L)
	1866	IF(XT.GT.XE) GO TO 420
	1867	IF((XE-XT).LT.FSH) GO TO 415
	1868	L=L+1
	1869	GO TO 410
	1870	415 A(J)=Y(L)
	1871	GO TO 425
	1872	420 IF((XT-XE).LT.FSH) GO TO 415
	1873	GRID(1)=FT
	1874	A(J)=GEE(1,NZ,GRID,PI2,X,Y,AD)
	1875	425 CONTINUE
	1876	IF(L.GT.1) L=L-1
	1877	430 CONTINUE
	1878	GRID(1)=GTEMP
	1879	DDEN=PI2/CN
	1880	DO 510 J=1,NFCNS
	1881	DTEMP=0.0
	1882	DNUM=(J-1)*DDEN
	1883	IF(NM1.LT.1) GO TO 505
	1884	DO 500 K=1,NM1
	1885	500 DTEMP=DTEMP+A(K+1)DCOS(DNUMK)
	1886	505 DTEMP=2.0*DTEMP+A(1)
	1887	510 ALPHA(J)=DTEMP
	1888	DO 550 J=2,NFCNS
	1889	550 ALPHA(J)=2*ALPHA(J)/CN
	1890	ALPHA(1)=ALPHA(1)/CN
	1891	IF(KKK.EQ.1) GO TO 545
	1892	P(1)=2.0ALPHA(NFCNS)BB+ALPHA(NM1)
	1893	P(2)=2.0AAALPHA(NFCNS)
	1894	Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
	1895	DO 540 J=2,NM1
	1896	IF(J.LT.NM1) GO TO 515
	1897	AA=0.5*AA
	1898	BB=0.5*BB
	1899	515 CONTINUE
	1900	P(J+1)=0.0
	1901	DO 520 K=1,J
	1902	A(K)=P(K)
	1903	520 P(K)=2.0BBA(K)
	1904	P(2)=P(2)+A(1)2.0AA
	1905	JM1=J-1
	1906	DO 525 K=1,JM1
	1907	525 P(K)=P(K)+Q(K)+AA*A(K+1)
	1908	JP1=J+1
	1909	DO 530 K=3,JP1
	1910	530 P(K)=P(K)+AA*A(K-1)
	1911	IF(J.EQ.NM1) GO TO 540
	1912	DO 535 K=1,J
	1913	535 Q(K)=-A(K)
	1914	Q(1)=Q(1)+ALPHA(NFCNS-1-J)
	1915	540 CONTINUE
	1916	DO 543 J=1,NFCNS
	1917	543 ALPHA(J)=P(J)
	1918	545 CONTINUE
	1919	IF(NFCNS.GT.3) RETURN
	1920	ALPHA(NFCNS+1)=0.0
	1921	ALPHA(NFCNS+2)=0.0
	1922	END SUBROUTINE remez
	1923
	1924	DOUBLE PRECISION FUNCTION D(K,N,M,X)
	1925	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
	1926	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
	1927	DIMENSION DES(1045),GRID(1045),WT(1045)
	1928	DOUBLE PRECISION AD,DEV,X,Y
	1929	DOUBLE PRECISION Q
	1930	DOUBLE PRECISION PI2
	1931	D = 1.0
	1932	Q = X(K)
	1933	DO 3 L = 1,M
	1934	DO 2 J = L,N,M
	1935	IF(J-K) 1,2,1
	1936	1 D = 2.0D(Q-X(J))
	1937	2 CONTINUE
	1938	3 CONTINUE
	1939	D = 1.0/D
	1940	END FUNCTION D
	1941
	1942
	1943	DOUBLE PRECISION FUNCTION GEE(K,N,GRID,PI2,X,Y,AD)
	1944	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
	1945	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
	1946	DIMENSION DES(1045),GRID(1045),WT(1045)
	1947	DOUBLE PRECISION AD,DEV,X,Y
	1948	DOUBLE PRECISION P,C,D,XF
	1949	DOUBLE PRECISION PI2
	1950	P = 0.0
	1951	XF = GRID(K)
	1952	XF = DCOS(PI2*XF)
	1953	D = 0.0
	1954	DO 1 J =1,N
	1955	C = XF-X(J)
	1956	C = AD(J)/C
	1957	D = D+C
	1958	P = P+C*Y(J)
	1959	1 CONTINUE
	1960	GEE = P/D
	1961	END FUNCTION GEE
	1962
	1963	#endif

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

Download in other formats:

Original Format