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(N*OMEGAC*DELTAT)/(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(N*PI/(NSTEPS+1)) / ( N*PI/(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(N*PI/(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.50*COS( N*PI/(NSTEPS)) & |
---|
716 | + 0.08*COS(2*N*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 = ALPHA*SQRT(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+Y*P7))))) |
---|
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(X0**2-1))**(FLOAT(N-1)) |
---|
874 | TERM2 = (X0 - SQRT(X0**2-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 = X0*COS(I*PI/N) |
---|
886 | CALL CHEBY(T,NM1,ARG) |
---|
887 | TERM1 = T(NM1) |
---|
888 | TERM2 = COS(2*NT*PI*I/N) |
---|
889 | SUM = SUM + 2*TERM1*TERM2 |
---|
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 = 2*PI*ABS(DELTAT/TAUS) |
---|
943 | X0 = 1/COS(THETAS/2) |
---|
944 | TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1)) |
---|
945 | TERM2 = (X0 - SQRT(X0**2-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 = X0*COS(I*PI/N) |
---|
959 | CALL CHEBY(T,NM1,ARG) |
---|
960 | TERM1 = T(NM1) |
---|
961 | TERM2 = COS(2*NT*PI*I/N) |
---|
962 | SUM = SUM + R*2*TERM1*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) = 2*X*T(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 = 2*PI*DTT/(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 = (MUC*SIGMA/(1+MUC*SIGMA))**NORD |
---|
1184 | ZN = (1-MUC*SIGMA)/(1+MUC*SIGMA) |
---|
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.5*ALPHA(1)+0.25*ALPHA(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.5*ALPHA(1)-0.25*ALPHA(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.5*ALPHA(1)-0.25*ALPHA(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(DNUM*K) |
---|
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.0*ALPHA(NFCNS)*BB+ALPHA(NM1) |
---|
1893 | P(2)=2.0*AA*ALPHA(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.0*BB*A(K) |
---|
1904 | P(2)=P(2)+A(1)*2.0*AA |
---|
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.0*D*(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 |
---|