Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dfi.F @ 990

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 81.5 KB

Line
1	SUBROUTINE dfi_accumulate( grid )
2
3	USE module_domain, ONLY : 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	CHARACTER*80 mess
15
16	! Input data.
17
18	TYPE(domain) , POINTER :: grid
19
20	IF ( grid%dfi_opt .EQ. DFI_NODFI .OR. grid%dfi_stage .EQ. DFI_FST ) RETURN
21
22	#if (EM_CORE == 1)
23
24
25	hn = grid%hcoeff(grid%itimestep+1)
26
27	! accumulate dynamic variables
28	grid%dfi_mu(:,:) = grid%dfi_mu(:,:) + grid%mu_2(:,:) * hn
29	grid%dfi_u(:,:,:) = grid%dfi_u(:,:,:) + grid%u_2(:,:,:) * hn
30	grid%dfi_v(:,:,:) = grid%dfi_v(:,:,:) + grid%v_2(:,:,:) * hn
31	grid%dfi_w(:,:,:) = grid%dfi_w(:,:,:) + grid%w_2(:,:,:) * hn
32	grid%dfi_ww(:,:,:) = grid%dfi_ww(:,:,:) + grid%ww(:,:,:) * hn
33	grid%dfi_t(:,:,:) = grid%dfi_t(:,:,:) + grid%t_2(:,:,:) * hn
34	grid%dfi_phb(:,:,:) = grid%dfi_phb(:,:,:) + grid%phb(:,:,:) * hn
35	grid%dfi_ph0(:,:,:) = grid%dfi_ph0(:,:,:) + grid%ph0(:,:,:) * hn
36	grid%dfi_php(:,:,:) = grid%dfi_php(:,:,:) + grid%php(:,:,:) * hn
37	grid%dfi_p(:,:,:) = grid%dfi_p(:,:,:) + grid%p(:,:,:) * hn
38	grid%dfi_ph(:,:,:) = grid%dfi_ph(:,:,:) + grid%ph_2(:,:,:) * hn
39	grid%dfi_tke(:,:,:) = grid%dfi_tke(:,:,:) + grid%tke_2(:,:,:) * hn
40	grid%dfi_al(:,:,:) = grid%dfi_al(:,:,:) + grid%al(:,:,:) * hn
41	grid%dfi_alt(:,:,:) = grid%dfi_alt(:,:,:) + grid%alt(:,:,:) * hn
42	grid%dfi_pb(:,:,:) = grid%dfi_pb(:,:,:) + grid%pb(:,:,:) * hn
43	! neg. check on hydrometeor and scalar variables
44	grid%moist(:,:,:,:) = max(0.,grid%moist(:,:,:,:))
45	grid%dfi_scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:))
46	#if ( WRF_DFI_RADAR == 1 )
47	IF ( grid%dfi_radar .EQ. 0 ) then
48	grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
49	ELSE
50	grid%dfi_moist(:,:,:,P_QV) = grid%dfi_moist(:,:,:,P_QV) + grid%moist(:,:,:,P_QV) * hn
51	ENDIF
52	#else
53	grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
54	#endif
55	grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
56
57	! accumulate DFI coefficient
58	grid%hcoeff_tot = grid%hcoeff_tot + hn
59	#endif
60
61	#if (NMM_CORE == 1)
62	hn = grid%hcoeff(grid%ntsd+1)
63	grid%dfi_pd(:,:) = grid%dfi_pd(:,:) + grid%pd(:,:) * hn
64	grid%dfi_pint(:,:,:) = grid%dfi_pint(:,:,:) + grid%pint(:,:,:) * hn
65	grid%dfi_dwdt(:,:,:) = grid%dfi_dwdt(:,:,:) + grid%dwdt(:,:,:) * hn
66	grid%dfi_t(:,:,:) = grid%dfi_t(:,:,:) + grid%t(:,:,:) * hn
67	grid%dfi_q(:,:,:) = grid%dfi_q(:,:,:) + grid%q(:,:,:) * hn
68	grid%dfi_q2(:,:,:) = grid%dfi_q2(:,:,:) + grid%q2(:,:,:) * hn
69	grid%dfi_cwm(:,:,:) = grid%dfi_cwm(:,:,:) + grid%cwm(:,:,:) * hn
70	grid%dfi_u(:,:,:) = grid%dfi_u(:,:,:) + grid%u(:,:,:) * hn
71	grid%dfi_v(:,:,:) = grid%dfi_v(:,:,:) + grid%v(:,:,:) * hn
72	grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
73	grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
74	! accumulate DFI coefficient
75	grid%hcoeff_tot = grid%hcoeff_tot + hn
76	write(mess,*) 'grid%hcoeff_tot: ', grid%hcoeff_tot
77	call wrf_message(mess)
78	#endif
79
80	END SUBROUTINE dfi_accumulate
81
82	SUBROUTINE dfi_bck_init ( grid )
83
84	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
85	USE module_utility
86	USE module_state_description
87
88	IMPLICIT NONE
89
90	TYPE (domain) , POINTER :: grid
91	INTEGER rc
92	CHARACTER*80 mess
93
94	INTERFACE
95	SUBROUTINE Setup_Timekeeping(grid)
96	USE module_domain, ONLY : domain
97	TYPE (domain), POINTER :: grid
98	END SUBROUTINE Setup_Timekeeping
99
100	SUBROUTINE dfi_save_arrays(grid)
101	USE module_domain, ONLY : domain
102	TYPE (domain), POINTER :: grid
103	END SUBROUTINE dfi_save_arrays
104
105	SUBROUTINE dfi_clear_accumulation(grid)
106	USE module_domain, ONLY : domain
107	TYPE (domain), POINTER :: grid
108	END SUBROUTINE dfi_clear_accumulation
109
110	SUBROUTINE optfil_driver(grid)
111	USE module_domain, ONLY : domain
112	TYPE (domain), POINTER :: grid
113	END SUBROUTINE optfil_driver
114
115	SUBROUTINE start_domain(grid, allowed_to_read)
116	USE module_domain, ONLY : domain
117	TYPE (domain) :: grid
118	LOGICAL, INTENT(IN) :: allowed_to_read
119	END SUBROUTINE start_domain
120	END INTERFACE
121
122	grid%dfi_stage = DFI_BCK
123
124	! Negate time step
125	IF ( grid%time_step_dfi .gt. 0 ) THEN
126	CALL nl_set_time_step ( 1, -grid%time_step_dfi )
127	ELSE
128	CALL nl_set_time_step ( 1, -grid%time_step )
129	CALL nl_set_time_step_fract_num ( 1, -grid%time_step_fract_num )
130	ENDIF
131
132	CALL Setup_Timekeeping (grid)
133
134	! set physics options to zero
135	CALL nl_set_mp_physics( grid%id, 0 )
136	CALL nl_set_ra_lw_physics( grid%id, 0 )
137	CALL nl_set_ra_sw_physics( grid%id, 0 )
138	CALL nl_set_sf_surface_physics( grid%id, 0 )
139	CALL nl_set_sf_sfclay_physics( grid%id, 0 )
140	CALL nl_set_sf_urban_physics( grid%id, 0 )
141	CALL nl_set_bl_pbl_physics( grid%id, 0 )
142	CALL nl_set_cu_physics( grid%id, 0 )
143	CALL nl_set_damp_opt( grid%id, 0 )
144	CALL nl_set_sst_update( grid%id, 0 )
145	CALL nl_set_fractional_seaice( grid%id, 0 )
146	CALL nl_set_gwd_opt( grid%id, 0 )
147	CALL nl_set_feedback( grid%id, 0 )
148	! set bc
149	#if (EM_CORE == 1)
150	CALL nl_set_diff_6th_opt( grid%id, 0 )
151	CALL nl_set_constant_bc(1, grid%constant_bc)
152	CALL nl_set_use_adaptive_time_step( grid%id, .false. )
153	#endif
154
155	#ifdef WRF_CHEM
156	! set chemistry option to zero
157	CALL nl_set_chem_opt (grid%id, 0)
158	CALL nl_set_aer_ra_feedback (grid%id, 0)
159	CALL nl_set_io_form_auxinput5 (grid%id, 0)
160	CALL nl_set_io_form_auxinput7 (grid%id, 0)
161	CALL nl_set_io_form_auxinput8 (grid%id, 0)
162	#endif
163
164	! set diffusion to zero for backward integration
165
166	#if (EM_CORE == 1)
167	CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
168	CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
169	IF ( grid%moist_adv_opt == 2 ) THEN
170	CALL nl_set_moist_adv_opt( grid%id, 0)
171	ENDIF
172	#endif
173
174	#if (NMM_CORE == 1)
175
176	!
177	! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
178	!
179	CALL nl_get_time_step( grid%id, grid%time_step )
180	if (grid%time_step .lt. 0) then
181	! DT =-DT
182
183	write(mess,*) 'changing signs for backward integration'
184	call wrf_message(mess)
185	grid%CPGFV = -grid%CPGFV
186	grid%EN = -grid%EN
187	grid%ENT = -grid%ENT
188	grid%F4D = -grid%F4D
189	grid%F4Q = -grid%F4Q
190	grid%EF4T = -grid%EF4T
191
192	grid%EM(:) = -grid%EM(:)
193	grid%EMT(:) = -grid%EMT(:)
194	grid%F4Q2(:) = -grid%F4Q2(:)
195
196	grid%WPDAR(:,:)= -grid%WPDAR(:,:)
197	grid%CPGFU(:,:)= -grid%CPGFU(:,:)
198	grid%CURV(:,:)= -grid%CURV(:,:)
199	grid%FCP(:,:)= -grid%FCP(:,:)
200	grid%FAD(:,:)= -grid%FAD(:,:)
201	grid%F(:,:)= -grid%F(:,:)
202	endif
203	#endif
204
205	grid%start_subtime = domain_get_start_time ( grid )
206	grid%stop_subtime = domain_get_stop_time ( grid )
207
208	CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
209
210	CALL dfi_save_arrays ( grid )
211	CALL dfi_clear_accumulation( grid )
212	CALL optfil_driver(grid)
213
214	!tgs need to call start_domain here to reset bc initialization for negative dt
215	CALL start_domain ( grid , .TRUE. )
216
217	END SUBROUTINE dfi_bck_init
218
219
220	SUBROUTINE dfi_fwd_init ( grid )
221
222	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
223	USE module_utility
224	USE module_state_description
225
226	IMPLICIT NONE
227
228	TYPE (domain) , POINTER :: grid
229	INTEGER rc
230	CHARACTER*80 mess
231
232	INTERFACE
233	SUBROUTINE Setup_Timekeeping(grid)
234	USE module_domain, ONLY : domain
235	TYPE (domain), POINTER :: grid
236	END SUBROUTINE Setup_Timekeeping
237
238	SUBROUTINE dfi_save_arrays(grid)
239	USE module_domain, ONLY : domain
240	TYPE (domain), POINTER :: grid
241	END SUBROUTINE dfi_save_arrays
242
243	SUBROUTINE dfi_clear_accumulation(grid)
244	USE module_domain, ONLY : domain
245	TYPE (domain), POINTER :: grid
246	END SUBROUTINE dfi_clear_accumulation
247
248	SUBROUTINE optfil_driver(grid)
249	USE module_domain, ONLY : domain
250	TYPE (domain), POINTER :: grid
251	END SUBROUTINE optfil_driver
252
253	SUBROUTINE start_domain(grid, allowed_to_read)
254	USE module_domain, ONLY : domain
255	TYPE (domain) :: grid
256	LOGICAL, INTENT(IN) :: allowed_to_read
257	END SUBROUTINE start_domain
258	END INTERFACE
259
260	grid%dfi_stage = DFI_FWD
261
262	! for Setup_Timekeeping to use when setting the clock
263
264	IF ( grid%time_step_dfi .gt. 0 ) THEN
265	CALL nl_set_time_step ( grid%id, grid%time_step_dfi )
266	ELSE
267
268	CALL nl_get_time_step( grid%id, grid%time_step )
269	CALL nl_get_time_step_fract_num( grid%id, grid%time_step_fract_num )
270
271	grid%time_step = abs(grid%time_step)
272	CALL nl_set_time_step( grid%id, grid%time_step )
273
274	grid%time_step_fract_num = abs(grid%time_step_fract_num)
275	CALL nl_set_time_step_fract_num( grid%id, grid%time_step_fract_num )
276
277	ENDIF
278
279	grid%itimestep=0
280	grid%xtime=0.
281
282	! reset physics options to normal
283	CALL nl_set_mp_physics( grid%id, grid%mp_physics)
284	CALL nl_set_ra_lw_physics( grid%id, grid%ra_lw_physics)
285	CALL nl_set_ra_sw_physics( grid%id, grid%ra_sw_physics)
286	CALL nl_set_sf_surface_physics( grid%id, grid%sf_surface_physics)
287	CALL nl_set_sf_sfclay_physics( grid%id, grid%sf_sfclay_physics)
288	CALL nl_set_sf_urban_physics( grid%id, grid%sf_urban_physics)
289	CALL nl_set_bl_pbl_physics( grid%id, grid%bl_pbl_physics)
290	CALL nl_set_cu_physics( grid%id, grid%cu_physics)
291	CALL nl_set_damp_opt( grid%id, grid%damp_opt)
292	CALL nl_set_sst_update( grid%id, 0)
293	CALL nl_set_fractional_seaice( grid%id, grid%fractional_seaice)
294	CALL nl_set_gwd_opt( grid%id, grid%gwd_opt)
295	CALL nl_set_feedback( grid%id, 0 )
296	#if (EM_CORE == 1)
297	CALL nl_set_diff_6th_opt( grid%id, grid%diff_6th_opt)
298	CALL nl_set_use_adaptive_time_step( grid%id, .false. )
299	! set bc
300	CALL nl_set_constant_bc(grid%id, grid%constant_bc)
301	#endif
302
303	#if (NMM_CORE == 1)
304	!
305	! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
306	!
307	!mptest CALL nl_get_time_step( grid%id, grid%time_step )
308
309
310	!!! here need to key off something else being the "wrong" sign
311	if (grid%cpgfv .gt. 0) then
312
313	write(mess,*) 'changing signs for forward integration'
314	call wrf_message(mess)
315	grid%CPGFV = -grid%CPGFV
316	grid%EN = -grid%EN
317	grid%ENT = -grid%ENT
318	grid%F4D = -grid%F4D
319	grid%F4Q = -grid%F4Q
320	grid%EF4T = -grid%EF4T
321
322	grid%EM(:) = -grid%EM(:)
323	grid%EMT(:) = -grid%EMT(:)
324	grid%F4Q2(:) = -grid%F4Q2(:)
325
326	grid%WPDAR(:,:)= -grid%WPDAR(:,:)
327	grid%CPGFU(:,:)= -grid%CPGFU(:,:)
328	grid%CURV(:,:)= -grid%CURV(:,:)
329	grid%FCP(:,:)= -grid%FCP(:,:)
330	grid%FAD(:,:)= -grid%FAD(:,:)
331	grid%F(:,:)= -grid%F(:,:)
332	endif
333	#endif
334
335
336	!#ifdef WRF_CHEM
337	! ! reset chem option to normal
338	! CALL nl_set_chem_opt( grid%id, grid%chem_opt)
339	! CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
340	!#endif
341
342	#if (EM_CORE == 1)
343	! reset km_opt to norlmal
344	CALL nl_set_km_opt( grid%id, grid%km_opt)
345	CALL nl_set_moist_adv_opt( grid%id, grid%moist_adv_opt)
346	#endif
347
348
349	CALL Setup_Timekeeping (grid)
350	grid%start_subtime = domain_get_start_time ( head_grid )
351	grid%stop_subtime = domain_get_stop_time ( head_grid )
352
353	CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
354
355	IF ( grid%dfi_opt .EQ. DFI_DFL ) THEN
356	CALL dfi_save_arrays ( grid )
357	END IF
358	CALL dfi_clear_accumulation( grid )
359	CALL optfil_driver(grid)
360
361	!tgs need to call it here to reset bc initialization for positive time_step
362	CALL start_domain ( grid , .TRUE. )
363
364	END SUBROUTINE dfi_fwd_init
365
366	SUBROUTINE dfi_fst_init ( grid )
367
368	USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time
369	USE module_state_description
370
371	IMPLICIT NONE
372
373	TYPE (domain) , POINTER :: grid
374	CHARACTER (LEN=80) :: wrf_error_message
375
376	INTERFACE
377	SUBROUTINE Setup_Timekeeping(grid)
378	USE module_domain, ONLY : domain
379	TYPE (domain), POINTER :: grid
380	END SUBROUTINE Setup_Timekeeping
381
382	SUBROUTINE dfi_save_arrays(grid)
383	USE module_domain, ONLY : domain
384	TYPE (domain), POINTER :: grid
385	END SUBROUTINE dfi_save_arrays
386
387	SUBROUTINE dfi_clear_accumulation(grid)
388	USE module_domain, ONLY : domain
389	TYPE (domain), POINTER :: grid
390	END SUBROUTINE dfi_clear_accumulation
391
392	SUBROUTINE optfil_driver(grid)
393	USE module_domain, ONLY : domain
394	TYPE (domain), POINTER :: grid
395	END SUBROUTINE optfil_driver
396
397	SUBROUTINE start_domain(grid, allowed_to_read)
398	USE module_domain, ONLY : domain
399	TYPE (domain) :: grid
400	LOGICAL, INTENT(IN) :: allowed_to_read
401	END SUBROUTINE start_domain
402	END INTERFACE
403
404	grid%dfi_stage = DFI_FST
405
406	! reset time_step to normal and adaptive_time_step
407	CALL nl_set_time_step( grid%id, grid%time_step )
408
409	grid%itimestep=0
410	grid%xtime=0. ! BUG: This will probably not work for all DFI options
411	! only use adaptive time stepping for forecast
412	#if (EM_CORE == 1)
413	if (grid%id == 1) then
414	CALL nl_set_use_adaptive_time_step( 1, grid%use_adaptive_time_step )
415	endif
416	CALL nl_set_sst_update( grid%id, grid%sst_update)
417	! reset to normal bc
418	CALL nl_set_constant_bc(grid%id, .false.)
419	#endif
420	CALL nl_set_feedback( grid%id, grid%feedback )
421
422	#ifdef WRF_CHEM
423	! reset chem option to normal
424	CALL nl_set_chem_opt( grid%id, grid%chem_opt)
425	CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
426	CALL nl_set_io_form_auxinput5 (grid%id, grid%io_form_auxinput5)
427	CALL nl_set_io_form_auxinput7 (grid%id, grid%io_form_auxinput7)
428	CALL nl_set_io_form_auxinput8 (grid%id, grid%io_form_auxinput8)
429	#endif
430
431
432	CALL Setup_Timekeeping (grid)
433	grid%start_subtime = domain_get_start_time ( grid )
434	grid%stop_subtime = domain_get_stop_time ( grid )
435
436	CALL start_domain ( grid , .TRUE. )
437
438	END SUBROUTINE dfi_fst_init
439
440
441	SUBROUTINE dfi_write_initialized_state( grid )
442
443	! Driver layer
444	USE module_domain, ONLY : domain, head_grid
445	USE module_io_domain
446	USE module_configure, ONLY : grid_config_rec_type, model_config_rec
447
448	IMPLICIT NONE
449
450	TYPE (domain) , POINTER :: grid
451	INTEGER :: fid, ierr
452	CHARACTER (LEN=80) :: wrf_error_message
453	CHARACTER (LEN=80) :: rstname
454	CHARACTER (LEN=132) :: message
455
456	TYPE (grid_config_rec_type) :: config_flags
457
458	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
459
460	WRITE (wrf_err_message,'(A,I4)') 'Writing out initialized model state'
461	CALL wrf_message(TRIM(wrf_err_message))
462
463	WRITE (rstname,'(A,I2.2)')'wrfinput_initialized_d',grid%id
464	CALL open_w_dataset ( fid, TRIM(rstname), grid, config_flags, output_input, "DATASET=INPUT", ierr )
465	IF ( ierr .NE. 0 ) THEN
466	WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
467	CALL WRF_ERROR_FATAL ( message )
468	END IF
469	CALL output_input ( fid, grid, config_flags, ierr )
470	CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
471
472	END SUBROUTINE dfi_write_initialized_state
473
474	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
475	! DFI array reset group of functions
476	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
477
478	SUBROUTINE wrf_dfi_array_reset ( )
479
480	USE module_domain, ONLY : domain, head_grid, set_current_grid_ptr
481
482	IMPLICIT NONE
483
484	INTERFACE
485	RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid)
486	USE module_domain, ONLY : domain
487	TYPE (domain), POINTER :: grid
488	END SUBROUTINE dfi_array_reset_recurse
489	END INTERFACE
490
491	! Copy filtered arrays back into state arrays in grid structure, and
492	! restore original land surface fields
493
494	CALL dfi_array_reset_recurse(head_grid)
495
496	CALL set_current_grid_ptr( head_grid )
497
498	END SUBROUTINE wrf_dfi_array_reset
499
500	SUBROUTINE dfi_array_reset( grid )
501
502	USE module_domain, ONLY : domain
503	! USE module_configure
504	! USE module_driver_constants
505	! USE module_machine
506	! USE module_dm
507	USE module_model_constants
508	USE module_state_description
509
510	IMPLICIT NONE
511
512	INTEGER :: its, ite, jts, jte, kts, kte, &
513	i, j, k
514
515	! Input data.
516	TYPE(domain) , POINTER :: grid
517
518	! local
519	! real p1000mb,eps,svp1,svp2,svp3,svpt0
520	real eps
521	! parameter (p1000mb = 1.e+05, eps=0.622,svp1=0.6112,svp3=29.65,svpt0=273.15)
522	parameter (eps=0.622)
523	REAL es,qs,pol,tx,temp,pres,rslf
524	CHARACTER*80 mess
525
526	IF ( grid%dfi_opt .EQ. DFI_NODFI ) RETURN
527
528
529	! Set dynamic variables
530	! divide by total DFI coefficient
531
532	#if (EM_CORE == 1)
533	grid%mu_2(:,:) = grid%dfi_mu(:,:) / grid%hcoeff_tot
534	grid%u_2(:,:,:) = grid%dfi_u(:,:,:) / grid%hcoeff_tot
535	grid%v_2(:,:,:) = grid%dfi_v(:,:,:) / grid%hcoeff_tot
536	grid%w_2(:,:,:) = grid%dfi_w(:,:,:) / grid%hcoeff_tot
537	grid%ww(:,:,:) = grid%dfi_ww(:,:,:) / grid%hcoeff_tot
538	grid%t_2(:,:,:) = grid%dfi_t(:,:,:) / grid%hcoeff_tot
539	grid%phb(:,:,:) = grid%dfi_phb(:,:,:) / grid%hcoeff_tot
540	grid%ph0(:,:,:) = grid%dfi_ph0(:,:,:) / grid%hcoeff_tot
541	grid%php(:,:,:) = grid%dfi_php(:,:,:) / grid%hcoeff_tot
542	grid%p(:,:,:) = grid%dfi_p(:,:,:) / grid%hcoeff_tot
543	grid%ph_2(:,:,:) = grid%dfi_ph(:,:,:) / grid%hcoeff_tot
544	grid%tke_2(:,:,:) = grid%dfi_tke(:,:,:) / grid%hcoeff_tot
545	grid%al(:,:,:) = grid%dfi_al(:,:,:) / grid%hcoeff_tot
546	grid%alt(:,:,:) = grid%dfi_alt(:,:,:) / grid%hcoeff_tot
547	grid%pb(:,:,:) = grid%dfi_pb(:,:,:) / grid%hcoeff_tot
548	#if ( WRF_DFI_RADAR == 1 )
549	IF ( grid%dfi_radar .EQ. 0 ) then ! tgs no radar assimilation
550	grid%moist(:,:,:,:) = max(0.,grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot)
551	ELSE
552	grid%moist(:,:,:,P_QV) = max(0.,grid%dfi_moist(:,:,:,P_QV) / grid%hcoeff_tot)
553	ENDIF
554	#else
555	grid%moist(:,:,:,:) = max(0.,grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot)
556	#endif
557	grid%scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot)
558
559	! restore initial fields
560	grid%SNOW (:,:) = grid%dfi_SNOW (:,:)
561	grid%SNOWH (:,:) = grid%dfi_SNOWH (:,:)
562	grid%SNOWC (:,:) = grid%dfi_SNOWC (:,:)
563	grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
564	grid%TSK (:,:) = grid%dfi_TSK (:,:)
565
566	grid%TSLB (:,:,:) = grid%dfi_TSLB (:,:,:)
567	grid%SMOIS (:,:,:) = grid%dfi_SMOIS (:,:,:)
568	IF ( grid%sf_surface_physics .EQ. 3 ) then
569	grid%QVG (:,:) = grid%dfi_QVG (:,:)
570	grid%TSNAV (:,:) = grid%dfi_TSNAV (:,:)
571	grid%SOILT1(:,:) = grid%dfi_SOILT1(:,:)
572	grid%SMFR3D(:,:,:) = grid%dfi_SMFR3D (:,:,:)
573	grid%KEEPFR3DFLAG(:,:,:) = grid%dfi_KEEPFR3DFLAG(:,:,:)
574	ENDIF
575
576	! restore analized hydrometeor fileds
577	#if ( WRF_DFI_RADAR == 1 )
578	IF ( grid%dfi_radar .EQ. 1 ) then
579	! grid%moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) !tgs
580	grid%moist(:,:,:,P_QC) = grid%dfi_moist(:,:,:,P_QC) !tgs
581	grid%moist(:,:,:,P_QR) = grid%dfi_moist(:,:,:,P_QR) !tgs
582	grid%moist(:,:,:,P_QI) = grid%dfi_moist(:,:,:,P_QI) !tgs
583	grid%moist(:,:,:,P_QS) = grid%dfi_moist(:,:,:,P_QS) !tgs
584	grid%moist(:,:,:,P_QG) = grid%dfi_moist(:,:,:,P_QG) !tgs
585
586	if(grid%dfi_stage .EQ. DFI_FWD) then
587	!tgs change QV to restore initial RH field after the diabatic DFI
588	its = grid%sp31 ; ite = grid%ep31 ;
589	kts = grid%sp32 ; kte = grid%ep32 ;
590	jts = grid%sp33 ; jte = grid%ep33 ;
591	DO j=jts,jte
592	DO i=its,ite
593	do k = kts , kte
594	temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
595	** (r_d / Cp)
596	pres = grid%p(i,k,j)+grid%pb(i,k,j)
597	!tgs rslf - function to compute qs from Thompson microphysics
598	qs = rslf(pres, temp)
599
600	! if(i.eq. 178 .and. j.eq. 148 .and. k.eq.11) then
601	! print *,'temp,pres,qs-thomp',temp,pres,qs
602	! endif
603
604	IF(grid%moist(i,k,j,P_QC) .GT. 1.e-6 .or. &
605	grid%moist(i,k,j,P_QI) .GT. 1.e-6) THEN
606	grid%moist (i,k,j,P_QV) = MAX(0.,grid%dfi_rh(i,k,j)*QS)
607	ENDIF
608
609	! if(i.eq. 178 .and. j.eq. 148 .and. k.eq.11) then
610	! print *,'temp,pres,qs,grid%moist (i,k,j,P_QV)',temp,pres,qs, &
611	! grid%moist(i,k,j,P_QV)
612	! endif
613	enddo
614	ENDDO
615	ENDDO
616	endif
617
618	ENDIF
619	#endif
620	#endif
621
622	#if (NMM_CORE == 1)
623	write(mess,*) ' divide by grid%hcoeff_tot: ', grid%hcoeff_tot
624	call wrf_message(mess)
625	if (grid%hcoeff_tot .lt. 0.0001) then
626	call wrf_error_fatal("bad grid%hcoeff_tot")
627	endif
628	grid%pd(:,:) = grid%dfi_pd(:,:) / grid%hcoeff_tot
629	grid%pint(:,:,:) = grid%dfi_pint(:,:,:) / grid%hcoeff_tot
630	! grid%dwdt(:,:,:) = grid%dfi_dwdt(:,:,:) / grid%hcoeff_tot
631	grid%t(:,:,:) = grid%dfi_t(:,:,:) / grid%hcoeff_tot
632	grid%q(:,:,:) = grid%dfi_q(:,:,:) / grid%hcoeff_tot
633	grid%q2(:,:,:) = grid%dfi_q2(:,:,:) / grid%hcoeff_tot
634	grid%cwm(:,:,:) = grid%dfi_cwm(:,:,:) / grid%hcoeff_tot
635	grid%u(:,:,:) = grid%dfi_u(:,:,:) / grid%hcoeff_tot
636	grid%v(:,:,:) = grid%dfi_v(:,:,:) / grid%hcoeff_tot
637	grid%moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot
638	grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot
639
640	! restore initial fields
641	grid%SNOW(:,:) = grid%dfi_SNOW(:,:)
642	grid%SNOWH(:,:) = grid%dfi_SNOWH(:,:)
643	! grid%SNOWC(:,:) = grid%dfi_SNOWC(:,:)
644	grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
645	grid%NMM_TSK(:,:) = grid%dfi_NMM_TSK(:,:)
646	! save soil fields
647	grid%STC(:,:,:) = grid%dfi_STC(:,:,:)
648	grid%SMC(:,:,:) = grid%dfi_SMC(:,:,:)
649	grid%SH2O(:,:,:) = grid%dfi_SH2O(:,:,:)
650	#endif
651
652
653	END SUBROUTINE dfi_array_reset
654
655	SUBROUTINE optfil_driver( grid )
656
657	USE module_domain, ONLY : domain
658	USE module_utility
659	! USE module_wrf_error
660	! USE module_timing
661	! USE module_date_time
662	! USE module_configure
663	USE module_state_description
664	USE module_model_constants
665
666	IMPLICIT NONE
667
668	TYPE (domain) , POINTER :: grid
669
670	! Local variables
671	integer :: nstep2, nstepmax, rundfi, i, rc,hr,min,sec
672	integer :: yr,jday
673	real :: timestep, tauc
674	TYPE(WRFU_TimeInterval) :: run_interval
675	CHARACTER*80 mess
676
677	timestep=abs(grid%dt)
678	run_interval = grid%stop_subtime - grid%start_subtime
679	CALL WRFU_TimeGet(grid%start_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
680	CALL WRFU_TimeGet(grid%stop_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
681
682	CALL WRFU_TimeIntervalGet( run_interval, S=rundfi, rc=rc )
683	rundfi = abs(rundfi)
684
685	nstep2= ceiling((1.0 + real(rundfi)/timestep) / 2.0)
686
687	! nstep2 is equal to a half of timesteps per initialization period,
688	! should not exceed nstepmax
689
690	tauc = real(grid%dfi_cutoff_seconds)
691
692	! Get DFI coefficient
693	grid%hcoeff(:) = 0.0
694	IF ( grid%dfi_nfilter < 0 .OR. grid%dfi_nfilter > 8 ) THEN
695	write(mess,*) 'Invalid filter specified in namelist.'
696	call wrf_message(mess)
697	ELSE
698	call dfcoef(nstep2-1, grid%dt, tauc, grid%dfi_nfilter, grid%hcoeff)
699	END IF
700
701	IF ( MOD(int(1.0 + real(rundfi)/timestep),2) /= 0 ) THEN
702	DO i=1,nstep2-1
703	grid%hcoeff(2*nstep2-i) = grid%hcoeff(i)
704	END DO
705	ELSE
706	DO i=1,nstep2
707	grid%hcoeff(2*nstep2-i+1) = grid%hcoeff(i)
708	END DO
709	END IF
710
711	END SUBROUTINE optfil_driver
712
713
714	SUBROUTINE dfi_clear_accumulation( grid )
715
716	USE module_domain, ONLY : domain
717	! USE module_configure
718	! USE module_driver_constants
719	! USE module_machine
720	! USE module_dm
721	! USE module_model_constants
722	USE module_state_description
723
724	IMPLICIT NONE
725
726	! Input data.
727	TYPE(domain) , POINTER :: grid
728
729	#if (EM_CORE == 1)
730	grid%dfi_mu(:,:) = 0.
731	grid%dfi_u(:,:,:) = 0.
732	grid%dfi_v(:,:,:) = 0.
733	grid%dfi_w(:,:,:) = 0.
734	grid%dfi_ww(:,:,:) = 0.
735	grid%dfi_t(:,:,:) = 0.
736	grid%dfi_phb(:,:,:) = 0.
737	grid%dfi_ph0(:,:,:) = 0.
738	grid%dfi_php(:,:,:) = 0.
739	grid%dfi_p(:,:,:) = 0.
740	grid%dfi_ph(:,:,:) = 0.
741	grid%dfi_tke(:,:,:) = 0.
742	grid%dfi_al(:,:,:) = 0.
743	grid%dfi_alt(:,:,:) = 0.
744	grid%dfi_pb(:,:,:) = 0.
745	#if ( WRF_DFI_RADAR == 1 )
746	IF ( grid%dfi_radar .EQ. 0 ) then
747	grid%dfi_moist(:,:,:,:) = 0.
748	ELSE
749	grid%dfi_moist(:,:,:,P_QV) = 0.
750	ENDIF
751	#else
752	grid%dfi_moist(:,:,:,:) = 0.
753	#endif
754	grid%dfi_scalar(:,:,:,:) = 0.
755	#endif
756
757	#if (NMM_CORE == 1)
758	grid%dfi_pd(:,:) = 0.
759	grid%dfi_pint(:,:,:) = 0.
760	grid%dfi_dwdt(:,:,:) = 0.
761	grid%dfi_t(:,:,:) = 0.
762	grid%dfi_q(:,:,:) = 0.
763	grid%dfi_q2(:,:,:) = 0.
764	grid%dfi_cwm(:,:,:) = 0.
765	grid%dfi_u(:,:,:) = 0.
766	grid%dfi_v(:,:,:) = 0.
767	grid%dfi_moist(:,:,:,:) = 0.
768	grid%dfi_scalar(:,:,:,:) = 0.
769	#endif
770
771	grid%hcoeff_tot = 0.0
772
773	END SUBROUTINE dfi_clear_accumulation
774
775
776	SUBROUTINE dfi_save_arrays( grid )
777
778	USE module_domain, ONLY : domain
779	! USE module_configure
780	! USE module_driver_constants
781	! USE module_machine
782	! USE module_dm
783	USE module_model_constants
784	USE module_state_description
785
786	IMPLICIT NONE
787
788	INTEGER :: its, ite, jts, jte, kts, kte, &
789	i, j, k
790
791	! Input data.
792	TYPE(domain) , POINTER :: grid
793	! local
794
795	REAL es,qs,pol,tx,temp,pres,rslf
796
797	#if (EM_CORE == 1)
798	! save surface 2-D fields
799	grid%dfi_SNOW(:,:) = grid%SNOW(:,:)
800	grid%dfi_SNOWH(:,:) = grid%SNOWH(:,:)
801	grid%dfi_SNOWC(:,:) = grid%SNOWC(:,:)
802	grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
803	grid%dfi_TSK(:,:) = grid%TSK(:,:)
804
805	! save soil fields
806	grid%dfi_TSLB(:,:,:) = grid%TSLB(:,:,:)
807	grid%dfi_SMOIS(:,:,:) = grid%SMOIS(:,:,:)
808	! RUC LSM only, need conditional
809	IF ( grid%sf_surface_physics .EQ. 3 ) then
810	grid%dfi_QVG(:,:) = grid%QVG(:,:)
811	grid%dfi_SOILT1(:,:) = grid%SOILT1(:,:)
812	grid%dfi_TSNAV(:,:) = grid%TSNAV(:,:)
813	grid%dfi_SMFR3D(:,:,:) = grid%SMFR3D(:,:,:)
814	grid%dfi_KEEPFR3DFLAG(:,:,:) = grid%KEEPFR3DFLAG(:,:,:)
815	ENDIF
816	#endif
817
818	#if (NMM_CORE == 1)
819	! save surface 2-D fields
820	grid%dfi_SNOW(:,:) = grid%SNOW(:,:)
821	grid%dfi_SNOWH(:,:) = grid%SNOWH(:,:)
822	! grid%dfi_SNOWC(:,:) = grid%SNOWC(:,:)
823	grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
824	grid%dfi_NMM_TSK(:,:) = grid%NMM_TSK(:,:)
825	! save soil fields
826	grid%dfi_STC(:,:,:) = grid%STC(:,:,:)
827	grid%dfi_SMC(:,:,:) = grid%SMC(:,:,:)
828	grid%dfi_SH2O(:,:,:) = grid%SH2O(:,:,:)
829	#endif
830
831
832	! save hydrometeor fields
833	#if (EM_CORE == 1)
834	#if ( WRF_DFI_RADAR == 1 )
835	IF ( grid%dfi_radar .EQ. 1 ) then !tgs
836	! grid%dfi_moist(:,:,:,:) = grid%moist(:,:,:,:)
837	grid%dfi_moist(:,:,:,P_QC) = grid%moist(:,:,:,P_QC)
838	grid%dfi_moist(:,:,:,P_QR) = grid%moist(:,:,:,P_QR)
839	grid%dfi_moist(:,:,:,P_QI) = grid%moist(:,:,:,P_QI)
840	grid%dfi_moist(:,:,:,P_QS) = grid%moist(:,:,:,P_QS)
841	grid%dfi_moist(:,:,:,P_QG) = grid%moist(:,:,:,P_QG)
842
843	if(grid%dfi_stage .EQ. DFI_BCK) then
844	! compute initial RH field to be reintroduced after the diabatic DFI
845	its = grid%sp31 ; ite = grid%ep31 ;
846	kts = grid%sp32 ; kte = grid%ep32 ;
847	jts = grid%sp33 ; jte = grid%ep33 ;
848	DO j=jts,jte
849	DO i=its,ite
850	do k = kts , kte
851	temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
852	** (r_d / Cp)
853	pres = grid%p(i,k,j)+grid%pb(i,k,j)
854	!tgs rslf - function to compute qs from Thompson microphysics
855	qs = rslf(pres, temp)
856	grid%dfi_rh (i,k,j) = MIN(1.,MAX(0.,grid%moist(i,k,j,P_QV)/qs))
857
858	!tgs saturation check for points with water or ice clouds
859	IF ((grid%moist (i,k,j,P_QC) .GT. 1.e-6 .or. &
860	grid%moist (i,k,j,P_QI) .GT. 1.e-6) .and. &
861	grid%dfi_rh (i,k,j) .lt. 1.) THEN
862	grid%dfi_rh (i,k,j)=1.
863	ENDIF
864
865	end do
866	END DO
867	ENDDO
868	endif
869
870	ENDIF
871	#endif
872	#endif
873
874	END SUBROUTINE dfi_save_arrays
875
876
877	SUBROUTINE dfcoef (NSTEPS,DT,TAUC,NFILT,H)
878	!
879	! calculate filter weights with selected window.
880	!
881	! peter lynch and xiang-yu huang
882	!
883	! ref: see hamming, r.w., 1989: digital filters,
884	! prentice-hall international. 3rd edition.
885	!
886	! input: nsteps - number of timesteps
887	! forward or backward.
888	! dt - time step in seconds.
889	! tauc - cut-off period in seconds.
890	! nfilt - indicator for selected filter.
891	!
892	! output: h - array(0:nsteps) with the
893	! required filter weights
894	!
895	!------------------------------------------------------------
896
897	implicit none
898
899	integer, intent(in) :: nsteps, nfilt
900	real , intent(in) :: dt, tauc
901	real, intent(out) :: h(1:nsteps+1)
902
903	! Local data
904
905	integer :: n
906	real :: pi, omegac, x, window, deltat
907	real :: hh(0:nsteps)
908
909	pi=4*ATAN(1.)
910	deltat=ABS(dt)
911
912	! windows are defined by a call and held in hh.
913
914	if ( nfilt .eq. -1) call debug (nsteps,h)
915
916	IF ( NFILT .EQ. 0 ) CALL UNIFORM (NSTEPS,HH)
917	IF ( NFILT .EQ. 1 ) CALL LANCZOS (NSTEPS,HH)
918	IF ( NFILT .EQ. 2 ) CALL HAMMING (NSTEPS,HH)
919	IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
920	IF ( NFILT .EQ. 4 ) CALL KAISER (NSTEPS,HH)
921	IF ( NFILT .EQ. 5 ) CALL POTTER2 (NSTEPS,HH)
922	IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
923
924	IF ( NFILT .LE. 6 ) THEN ! sinc-windowed filters
925
926	! calculate the cutoff frequency
927	OMEGAC = 2.*PI/TAUC
928
929	DO N=0,NSTEPS
930	WINDOW = HH(N)
931	IF ( N .EQ. 0 ) THEN
932	X = (OMEGAC*DELTAT/PI)
933	ELSE
934	X = SIN(NOMEGACDELTAT)/(N*PI)
935	END IF
936	HH(N) = X*WINDOW
937	END DO
938
939	! normalize the sums to be unity
940	CALL NORMLZ(HH,NSTEPS)
941
942	DO N=0,NSTEPS
943	H(N+1) = HH(NSTEPS-N)
944	END DO
945
946	ELSE IF ( NFILT .EQ. 7 ) THEN ! dolph filter
947
948	CALL DOLPH(DT,TAUC,NSTEPS,H)
949
950	ELSE IF ( NFILT .EQ. 8 ) THEN ! 2nd order, 2nd type quick start filter (RHO)
951
952	CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS)
953
954	END IF
955
956	RETURN
957
958	END SUBROUTINE dfcoef
959
960
961	SUBROUTINE NORMLZ(HH,NMAX)
962
963	! normalize the sum of hh to be unity
964
965	implicit none
966
967	integer, intent(in) :: nmax
968	real , dimension(0:nmax), intent(out) :: hh
969
970	! local data
971	real :: sumhh
972	integer :: n
973
974	SUMHH = HH(0)
975	DO N=1,NMAX
976	SUMHH = SUMHH + 2*HH(N)
977	ENDDO
978	DO N=0,NMAX
979	HH(N) = HH(N)/SUMHH
980	ENDDO
981
982	RETURN
983
984	END subroutine normlz
985
986
987	subroutine debug(nsteps, ww)
988
989	implicit none
990
991	integer, intent(in) :: nsteps
992	real , dimension(0:nsteps), intent(out) :: ww
993	integer :: n
994
995	do n=0,nsteps
996	ww(n)=0
997	end do
998
999	ww(int(nsteps/2))=1
1000
1001	return
1002
1003	end subroutine debug
1004
1005
1006	SUBROUTINE UNIFORM(NSTEPS,WW)
1007
1008	! define uniform or rectangular window function.
1009
1010	implicit none
1011
1012	integer, intent(in) :: nsteps
1013	real , dimension(0:nsteps), intent(out) :: ww
1014
1015	integer :: n
1016
1017	DO N=0,NSTEPS
1018	WW(N) = 1.
1019	ENDDO
1020
1021	RETURN
1022
1023	END subroutine uniform
1024
1025
1026	SUBROUTINE LANCZOS(NSTEPS,WW)
1027
1028	! define (genaralised) lanczos window function.
1029
1030	implicit none
1031
1032	integer, parameter :: nmax = 1000
1033	integer, intent(in) :: nsteps
1034	real , dimension(0:nmax), intent(out) :: ww
1035	integer :: n
1036	real :: power, pi, w
1037
1038	! (for the usual lanczos window, power = 1 )
1039	POWER = 1
1040
1041	PI=4*ATAN(1.)
1042	DO N=0,NSTEPS
1043	IF ( N .EQ. 0 ) THEN
1044	W = 1.0
1045	ELSE
1046	W = SIN(NPI/(NSTEPS+1)) / ( NPI/(NSTEPS+1))
1047	ENDIF
1048	WW(N) = W**POWER
1049	ENDDO
1050
1051	RETURN
1052
1053	END SUBROUTINE lanczos
1054
1055
1056	SUBROUTINE HAMMING(NSTEPS,WW)
1057
1058	! define (genaralised) hamming window function.
1059
1060	implicit none
1061
1062	integer, intent(in) :: nsteps
1063	real, dimension(0:nsteps) :: ww
1064	integer :: n
1065	real :: alpha, pi, w
1066
1067	! (for the usual hamming window, alpha=0.54,
1068	! for the hann window, alpha=0.50).
1069	ALPHA=0.54
1070
1071	PI=4*ATAN(1.)
1072	DO N=0,NSTEPS
1073	IF ( N .EQ. 0 ) THEN
1074	W = 1.0
1075	ELSE
1076	W = ALPHA + (1-ALPHA)COS(NPI/(NSTEPS))
1077	ENDIF
1078	WW(N) = W
1079	ENDDO
1080
1081	RETURN
1082
1083	END SUBROUTINE hamming
1084
1085
1086	SUBROUTINE BLACKMAN(NSTEPS,WW)
1087
1088	! define blackman window function.
1089
1090	implicit none
1091
1092	integer, intent(in) :: nsteps
1093	real, dimension(0:nsteps) :: ww
1094	integer :: n
1095
1096	real :: pi, w
1097
1098	PI=4*ATAN(1.)
1099	DO N=0,NSTEPS
1100	IF ( N .EQ. 0 ) THEN
1101	W = 1.0
1102	ELSE
1103	W = 0.42 + 0.50COS( NPI/(NSTEPS)) &
1104	+ 0.08COS(2N*PI/(NSTEPS))
1105	ENDIF
1106	WW(N) = W
1107	ENDDO
1108
1109	RETURN
1110
1111	END SUBROUTINE blackman
1112
1113
1114	SUBROUTINE KAISER(NSTEPS,WW)
1115
1116	! define kaiser window function.
1117
1118	implicit none
1119
1120	real, external :: bessi0
1121
1122	integer, intent(in) :: nsteps
1123	real, dimension(0:nsteps) :: ww
1124	integer :: n
1125	real :: alpha, xi0a, xn, as
1126
1127	ALPHA = 1
1128
1129	XI0A = BESSI0(ALPHA)
1130	DO N=0,NSTEPS
1131	XN = N
1132	AS = ALPHASQRT(1.-(XN/NSTEPS)*2)
1133	WW(N) = BESSI0(AS) / XI0A
1134	ENDDO
1135
1136	RETURN
1137
1138	END SUBROUTINE kaiser
1139
1140
1141	REAL FUNCTION BESSI0(X)
1142
1143	! from numerical recipes (press, et al.)
1144
1145	implicit none
1146
1147	real(8) :: Y
1148	real(8) :: P1 = 1.0d0
1149	real(8) :: P2 = 3.5156230D0
1150	real(8) :: P3 = 3.0899424D0
1151	real(8) :: P4 = 1.2067492D0
1152	real(8) :: P5 = 0.2659732D0
1153	real(8) :: P6 = 0.360768D-1
1154	real(8) :: P7 = 0.45813D-2
1155
1156	real*8 :: Q1 = 0.39894228D0
1157	real*8 :: Q2 = 0.1328592D-1
1158	real*8 :: Q3 = 0.225319D-2
1159	real*8 :: Q4 = -0.157565D-2
1160	real*8 :: Q5 = 0.916281D-2
1161	real*8 :: Q6 = -0.2057706D-1
1162	real*8 :: Q7 = 0.2635537D-1
1163	real*8 :: Q8 = -0.1647633D-1
1164	real*8 :: Q9 = 0.392377D-2
1165
1166	real :: x, ax
1167
1168
1169	IF (ABS(X).LT.3.75) THEN
1170	Y=(X/3.75)**2
1171	BESSI0=P1+Y(P2+Y(P3+Y(P4+Y(P5+Y(P6+YP7)))))
1172	ELSE
1173	AX=ABS(X)
1174	Y=3.75/AX
1175	BESSI0=(EXP(AX)/SQRT(AX))(Q1+Y(Q2+Y(Q3+Y(Q4 &
1176	+Y(Q5+Y(Q6+Y(Q7+Y(Q8+Y*Q9))))))))
1177	ENDIF
1178	RETURN
1179
1180	END FUNCTION bessi0
1181
1182
1183	SUBROUTINE POTTER2(NSTEPS,WW)
1184
1185	! define potter window function.
1186	! modified to fall off over twice the range.
1187
1188	implicit none
1189
1190	integer, intent(in) :: nsteps
1191	real, dimension(0:nsteps),intent(out) :: ww
1192	integer :: n
1193	real :: ck, sum, arg
1194
1195	! local data
1196	real, dimension(0:3) :: d
1197	real :: pi
1198	integer :: ip
1199
1200	d(0) = 0.35577019
1201	d(1) = 0.2436983
1202	d(2) = 0.07211497
1203	d(3) = 0.00630165
1204
1205	PI=4*ATAN(1.)
1206
1207	CK = 1.0
1208	DO N=0,NSTEPS
1209	IF (N.EQ.NSTEPS) CK = 0.5
1210	ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
1211	!min--- modification in next statement
1212	ARG = ARG/2.
1213	!min--- end of modification
1214	SUM = D(0)
1215	DO IP=1,3
1216	SUM = SUM + 2.D(IP)COS(ARG*FLOAT(IP))
1217	END DO
1218	WW(N) = CK*SUM
1219	END DO
1220
1221	RETURN
1222
1223	END SUBROUTINE potter2
1224
1225
1226	SUBROUTINE dolphwin(m, window)
1227
1228	! calculation of dolph-chebyshev window or, for short,
1229	! dolph window, using the expression in the reference:
1230	!
1231	! antoniou, andreas, 1993: digital filters: analysis,
1232	! design and applications. mcgraw-hill, inc., 689pp.
1233	!
1234	! the dolph window is optimal in the following sense:
1235	! for a given main-lobe width, the stop-band attenuation
1236	! is minimal; for a given stop-band level, the main-lobe
1237	! width is minimal.
1238	!
1239	! it is possible to specify either the ripple-ratio r
1240	! or the stop-band edge thetas.
1241
1242	IMPLICIT NONE
1243
1244	! Arguments
1245	INTEGER, INTENT(IN) :: m
1246	REAL, DIMENSION(0:M), INTENT(OUT) :: window
1247
1248	! local data
1249	REAL, DIMENSION(0:2*M) :: t
1250	REAL, DIMENSION(0:M) :: w, time
1251	REAL :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
1252	INTEGER :: n, nm1, nt, i
1253
1254	PI = 4*ATAN(1.D0)
1255	THETAS = 2*PI/M
1256
1257	N = 2*M+1
1258	NM1 = N-1
1259	X0 = 1/COS(THETAS/2)
1260
1261	TERM1 = (X0 + SQRT(X02-1))(FLOAT(N-1))
1262	TERM2 = (X0 - SQRT(X02-1))(FLOAT(N-1))
1263	RR = 0.5*(TERM1+TERM2)
1264	R = 1/RR
1265	DB = 20*LOG10(R)
1266	WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
1267	WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
1268	WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
1269
1270	DO NT=0,M
1271	SUM = RR
1272	DO I=1,M
1273	ARG = X0COS(IPI/N)
1274	CALL CHEBY(T,NM1,ARG)
1275	TERM1 = T(NM1)
1276	TERM2 = COS(2NTPI*I/N)
1277	SUM = SUM + 2TERM1TERM2
1278	ENDDO
1279	W(NT) = SUM/N
1280	TIME(NT) = NT
1281	ENDDO
1282
1283	! fill up the array for return
1284	DO NT=0,M
1285	WINDOW(NT) = W(NT)
1286	ENDDO
1287
1288	RETURN
1289
1290	END SUBROUTINE dolphwin
1291
1292
1293	SUBROUTINE dolph(deltat, taus, m, window)
1294
1295	! calculation of dolph-chebyshev window or, for short,
1296	! dolph window, using the expression in the reference:
1297	!
1298	! antoniou, andreas, 1993: digital filters: analysis,
1299	! design and applications. mcgraw-hill, inc., 689pp.
1300	!
1301	! the dolph window is optimal in the following sense:
1302	! for a given main-lobe width, the stop-band attenuation
1303	! is minimal; for a given stop-band level, the main-lobe
1304	! width is minimal.
1305
1306	IMPLICIT NONE
1307
1308	! Arguments
1309	INTEGER, INTENT(IN) :: m
1310	REAL, DIMENSION(0:M), INTENT(OUT) :: window
1311	REAL, INTENT(IN) :: deltat, taus
1312
1313	! local data
1314	integer, PARAMETER :: NMAX = 5000
1315	REAL, dimension(0:NMAX) :: t, w, time
1316	real, dimension(0:2*nmax) :: w2
1317	INTEGER :: NPRPE=0 ! no of pe
1318	CHARACTER*80 :: MES
1319
1320	real :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
1321	integer :: n, nm1, i, nt
1322
1323	PI = 4*ATAN(1.D0)
1324
1325	WRITE (mes,'(A,F8.2,A,F10.2)') 'In dolph, deltat = ',deltat,' taus = ',taus
1326	CALL wrf_message(TRIM(mes))
1327
1328	N = 2*M+1
1329	NM1 = N-1
1330
1331	THETAS = 2PIABS(DELTAT/TAUS)
1332	X0 = 1/COS(THETAS/2)
1333	TERM1 = (X0 + SQRT(X02-1))(FLOAT(N-1))
1334	TERM2 = (X0 - SQRT(X02-1))(FLOAT(N-1))
1335	RR = 0.5*(TERM1+TERM2)
1336	R = 1/RR
1337	DB = 20*LOG10(R)
1338
1339	WRITE (mes,'(A,2I8)') 'In dolph: M,N = ', M,N
1340	CALL wrf_message(TRIM(mes))
1341	WRITE (mes,'(A,F10.3)') 'In dolph: THETAS (STOP-BAND EDGE) = ', thetas
1342	CALL wrf_message(TRIM(mes))
1343	WRITE (mes,'(A,2F10.3)') 'In dolph: R,DB = ', R,DB
1344	CALL wrf_message(TRIM(mes))
1345
1346	DO NT=0,M
1347	SUM = 1
1348	DO I=1,M
1349	ARG = X0COS(IPI/N)
1350	CALL CHEBY(T,NM1,ARG)
1351	TERM1 = T(NM1)
1352	TERM2 = COS(2NTPI*I/N)
1353	SUM = SUM + R2TERM1*TERM2
1354	ENDDO
1355	W(NT) = SUM/N
1356	TIME(NT) = NT
1357	WRITE (mes,'(A,F10.6,2x,E17.7)') 'In dolph: TIME, W = ', TIME(NT), W(NT)
1358	CALL wrf_message(TRIM(mes))
1359	ENDDO
1360	! fill in the negative-time values by symmetry.
1361	DO NT=0,M
1362	W2(M+NT) = W(NT)
1363	W2(M-NT) = W(NT)
1364	ENDDO
1365
1366	! fill up the array for return
1367	SUMW = 0.
1368	DO NT=0,2*M
1369	SUMW = SUMW + W2(NT)
1370	ENDDO
1371	WRITE (mes,'(A,F10.4)') 'In dolph: SUM OF WEIGHTS W2 = ', sumw
1372	CALL wrf_message(TRIM(mes))
1373
1374	DO NT=0,M
1375	WINDOW(NT) = W2(NT)
1376	ENDDO
1377
1378	RETURN
1379
1380	END SUBROUTINE dolph
1381
1382
1383	SUBROUTINE cheby(t, n, x)
1384
1385	! calculate all chebyshev polynomials up to order n
1386	! for the argument value x.
1387
1388	! reference: numerical recipes, page 184, recurrence
1389	! t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) , n>=2.
1390
1391	IMPLICIT NONE
1392
1393	! Arguments
1394	INTEGER, INTENT(IN) :: n
1395	REAL, INTENT(IN) :: x
1396	REAL, DIMENSION(0:N) :: t
1397
1398	integer :: nn
1399
1400	T(0) = 1
1401	T(1) = X
1402	IF(N.LT.2) RETURN
1403	DO NN=2,N
1404	T(NN) = 2XT(NN-1) - T(NN-2)
1405	ENDDO
1406
1407	RETURN
1408
1409	END SUBROUTINE cheby
1410
1411
1412	SUBROUTINE rhofil(dt, tauc, norder, nstep, ictype, frow, nosize)
1413
1414	! RHO = recurssive high order.
1415	!
1416	! This routine calculates and returns the
1417	! Last Row, FROW, of the FILTER matrix.
1418	!
1419	! Input Parameters:
1420	! DT : Time Step in seconds
1421	! TAUC : Cut-off period (hours)
1422	! NORDER : Order of QS Filter
1423	! NSTEP : Number of step/Size of row.
1424	! ICTYPE : Initial Conditions
1425	! NOSIZE : Max. side of FROW.
1426	!
1427	! Working Fields:
1428	! ACOEF : X-coefficients of filter
1429	! BCOEF : Y-coefficients of filter
1430	! FILTER : Filter Matrix.
1431	!
1432	! Output Parameters:
1433	! FROW : Last Row of Filter Matrix.
1434	!
1435	! Note: Two types of initial conditions are permitted.
1436	! ICTYPE = 1 : Order increasing each row to NORDER.
1437	! ICTYPE = 2 : Order fixed at NORDER throughout.
1438	!
1439	! DOUBLE PRECISION USED THROUGHOUT.
1440
1441	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1442
1443	DOUBLE PRECISION MUC
1444
1445	! N.B. Single Precision for List Parameters.
1446	REAL, intent(in) :: DT,TAUC
1447
1448	! Space for the last row of FILTER.
1449	integer, intent(in) :: norder, nstep, ictype, nosize
1450	REAL , dimension(0:nosize), intent(out):: FROW
1451
1452	! Arrays for rho filter.
1453	integer, PARAMETER :: NOMAX=100
1454	real , dimension(0:NOMAX) :: acoef, bcoef
1455	real , dimension(0:NOMAX,0:NOMAX) :: filter
1456	! Working space.
1457	real , dimension(0:NOMAX) :: alpha, beta
1458
1459	real :: DTT
1460
1461	DTT = ABS(DT)
1462	PI = 2*DASIN(1.D0)
1463	IOTA = CMPLX(0.,1.)
1464
1465	! Filtering Parameters (derived).
1466	THETAC = 2PIDTT/(TAUC)
1467	MUC = tan(THETAC/2)
1468	FC = THETAC/(2*PI)
1469
1470	! Clear the arrays.
1471	DO NC=0,NOMAX
1472	ACOEF(NC) = 0.
1473	BCOEF(NC) = 0.
1474	ALPHA(NC) = 0.
1475	BETA (NC) = 0.
1476	FROW (NC) = 0.
1477	DO NR=0,NOMAX
1478	FILTER(NR,NC) = 0.
1479	ENDDO
1480	ENDDO
1481
1482	! Fill up the Filter Matrix.
1483	FILTER(0,0) = 1.
1484
1485	! Get the coefficients of the Filter.
1486	IF ( ICTYPE.eq.2 ) THEN
1487	CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
1488	ENDIF
1489
1490	DO 100 NROW=1,NSTEP
1491
1492	IF ( ICTYPE.eq.1 ) THEN
1493	NORD = MIN(NROW,NORDER)
1494	IF ( NORD.le.NORDER) THEN
1495	CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
1496	ENDIF
1497	ENDIF
1498
1499	DO K=0,NROW
1500	ALPHA(K) = ACOEF(NROW-K)
1501	IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
1502	ENDDO
1503
1504	! Correction for terms of negative index.
1505	IF ( ICTYPE.eq.2 ) THEN
1506	IF ( NROW.lt.NORDER ) THEN
1507	CN = 0.
1508	DO NN=NROW+1,NORDER
1509	CN = CN + (ACOEF(NN)+BCOEF(NN))
1510	ENDDO
1511	ALPHA(0) = ALPHA(0) + CN
1512	ENDIF
1513	ENDIF
1514
1515	! Check sum of ALPHAs and BETAs = 1
1516	SUMAB = 0.
1517	DO NN=0,NROW
1518	SUMAB = SUMAB + ALPHA(NN)
1519	IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
1520	ENDDO
1521
1522	DO KK=0,NROW-1
1523	SUMBF = 0.
1524	DO LL=0,NROW-1
1525	SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
1526	ENDDO
1527	FILTER(NROW,KK) = ALPHA(KK)+SUMBF
1528	ENDDO
1529	FILTER(NROW,NROW) = ALPHA(NROW)
1530
1531	! Check sum of row elements = 1
1532	SUMROW = 0.
1533	DO NN=0,NROW
1534	SUMROW = SUMROW + FILTER(NROW,NN)
1535	ENDDO
1536
1537	100 CONTINUE
1538
1539	DO NC=0,NSTEP
1540	FROW(NC) = FILTER(NSTEP,NC)
1541	ENDDO
1542
1543	RETURN
1544
1545	END SUBROUTINE rhofil
1546
1547
1548	SUBROUTINE rhocof(nord, nomax, muc, ca, cb)
1549
1550	! Get the coefficients of the RHO Filter
1551
1552	! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1553	IMPLICIT NONE
1554
1555	! Arguments
1556	integer, intent(in) :: nord, nomax
1557	real, dimension(0:nomax) :: ca, cb
1558
1559	! Functions
1560	double precision, external :: cnr
1561
1562	! Local variables
1563	INTEGER :: nn
1564	COMPLEX :: IOTA
1565	DOUBLE PRECISION :: MUC, ZN
1566	DOUBLE PRECISION :: pi, root2, rn, sigma, gain, sumcof
1567
1568	PI = 2*ASIN(1.)
1569	ROOT2 = SQRT(2.)
1570	IOTA = (0.,1.)
1571
1572	RN = 1./FLOAT(NORD)
1573	SIGMA = 1./( SQRT(2.**RN-1.) )
1574
1575	GAIN = (MUCSIGMA/(1+MUCSIGMA))**NORD
1576	ZN = (1-MUCSIGMA)/(1+MUCSIGMA)
1577
1578	DO NN=0,NORD
1579	CA(NN) = CNR(NORD,NN)*GAIN
1580	IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)(-ZN)*NN
1581	ENDDO
1582
1583	! Check sum of coefficients = 1
1584	SUMCOF = 0.
1585	DO NN=0,NORD
1586	SUMCOF = SUMCOF + CA(NN)
1587	IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
1588	ENDDO
1589
1590	RETURN
1591
1592	END SUBROUTINE RHOCOF
1593
1594
1595	DOUBLE PRECISION FUNCTION cnr(n,r)
1596
1597	! Binomial Coefficient (n,r).
1598
1599	! IMPLICIT DOUBLE PRECISION(C,X)
1600	IMPLICIT NONE
1601
1602	! Arguments
1603	INTEGER , intent(in) :: n, R
1604
1605	! Local variables
1606	INTEGER :: k
1607	DOUBLE PRECISION :: coeff, xn, xr, xk
1608
1609	IF ( R.eq.0 ) THEN
1610	CNR = 1.0
1611	RETURN
1612	ENDIF
1613	Coeff = 1.0
1614	XN = DFLOAT(N)
1615	XR = DFLOAT(R)
1616	DO K=1,R
1617	XK = DFLOAT(K)
1618	COEFF = COEFF * ( (XN-XR+XK)/XK )
1619	ENDDO
1620	CNR = COEFF
1621
1622	RETURN
1623
1624	END FUNCTION cnr
1625
1626
1627	SUBROUTINE optfil (grid,NH,DELTAT,NHMAX)
1628	!----------------------------------------------------------------------
1629
1630	! SUBROUTINE optfil (NH,DELTAT,TAUP,TAUS,LPRINT, &
1631	! H,NHMAX)
1632	!
1633	! - Huang and Lynch optimal filter
1634	! Monthly Weather Review, Feb 1993
1635	!----------------------------------------------------------
1636	! Input Parameters in List:
1637	! NH : Half-length of the Filter
1638	! DELTAT : Time-step (in seconds).
1639	! TAUP : Period of pass-band edge (hours).
1640	! TAUS : Period of stop-band edge (hours).
1641	! LPRINT : Logical switch for messages.
1642	! NHMAX : Maximum permitted Half-length.
1643	!
1644	! Output Parameters in List:
1645	! H : Impulse Response.
1646	! DP : Deviation in pass-band (db)
1647	! DS : Deviation in stop-band (db)
1648	!----------------------------------------------------------
1649	!
1650	USE module_domain, ONLY : domain
1651
1652	TYPE(domain) , POINTER :: grid
1653
1654	REAL,DIMENSION( 20) :: EDGE
1655	REAL,DIMENSION( 10) :: FX, WTX, DEVIAT
1656	REAL,DIMENSION(2*NHMAX+1) :: H
1657	logical LPRINT
1658	REAL, INTENT (IN) :: DELTAT
1659	INTEGER, INTENT (IN) :: NH, NHMAX
1660	!
1661	TAUP = 3.
1662	TAUS = 1.5
1663	LPRINT = .true.
1664	!initialize H array
1665
1666	NL=2*NHMAX+1
1667	do 101 n=1,NL
1668	H(n)=0.
1669	101 continue
1670
1671	NFILT = 2*NH+1
1672	print *,' start optfil, NFILT=', nfilt
1673
1674	!
1675	! 930325 PL & XYH : the upper limit is changed from 64 to 128.
1676	IF(NFILT.LE.0 .OR. NFILT.GT.128 ) THEN
1677	WRITE(6,*) 'NH=',NH
1678	CALL wrf_error_fatal (' Sorry, error 1 in call to OPTFIL ')
1679	ENDIF
1680	!
1681	! The following four should always be the same.
1682	JTYPE = 1
1683	NBANDS = 2
1684	!CC JPRINT = 0
1685	LGRID = 16
1686	!
1687	! calculate transition frequencies.
1688	DT = ABS(DELTAT)
1689	FS = DT/(TAUS*3600.)
1690	FP = DT/(TAUP*3600.)
1691	IF(FS.GT.0.5) then
1692	! print *,' FS too large in OPTFIL '
1693	CALL wrf_error_fatal (' FS too large in OPTFIL ')
1694	! return
1695	end if
1696	IF(FP.LT.0.0) then
1697	! print *, ' FP too small in OPTFIL '
1698	CALL wrf_error_fatal (' FP too small in OPTFIL ')
1699	! return
1700	end if
1701	!
1702	! Relative Weights in pass- and stop-bands.
1703	WTP = 1.0
1704	WTS = 1.0
1705	!
1706	!CC NOTE: (FP,FC,FS) is an arithmetic progression, so
1707	!CC (1/FS,1/FC,1/FP) is a harmonic one.
1708	!CC TAUP = 1/( (1/TAUC)-(1/DTAU) )
1709	!CC TAUS = 1/( (1/TAUC)+(1/DTAU) )
1710	!CC TAUC : Cut-off Period (hours).
1711	!CC DTAU : Transition half-width (hours).
1712	!CC FC = 1/TAUC ; DF = 1/DTAU
1713	!CC FP = FC - DF : FS = FC + DF
1714	!
1715	IF ( LPRINT ) THEN
1716	TAUC = 2./((1/TAUS)+(1/TAUP))
1717	DTAU = 2./((1/TAUS)-(1/TAUP))
1718	FC = DT/(TAUC*3600.)
1719	DF = DT/(DTAU*3600.)
1720	WRITE(6,*) ' DT ' , dt
1721	WRITE(6,*) ' TAUS, TAUP ' , TAUS,TAUP
1722	WRITE(6,*) ' TAUC, DTAU ' , TAUC,DTAU
1723	WRITE(6,*) ' FP, FS ' , FP, FS
1724	WRITE(6,*) ' FC, DF ' , FC, DF
1725	WRITE(6,*) ' WTS, WTP ' , WTS, WTP
1726	ENDIF
1727	!
1728	! Fill the control vectors for MCCPAR
1729	EDGE(1) = 0.0
1730	EDGE(2) = FP
1731	EDGE(3) = FS
1732	EDGE(4) = 0.5
1733	FX(1) = 1.0
1734	FX(2) = 0.0
1735	WTX(1) = WTP
1736	WTX(2) = WTS
1737
1738	CALL MCCPAR(NFILT,JTYPE,NBANDS,LPRINT,LGRID, &
1739	EDGE,FX,WTX,DEVIAT, h )
1740	!
1741	! Save the deviations in the pass- and stop-bands.
1742	DP = DEVIAT(1)
1743	DS = DEVIAT(2)
1744	!
1745	! Fill out the array H (only first half filled in MCCPAR).
1746	IF(MOD(NFILT,2).EQ.0) THEN
1747	NHALF = ( NFILT )/2
1748	ELSE
1749	NHALF = (NFILT+1)/2
1750	ENDIF
1751	DO 100 nn=1,NHALF
1752	H(NFILT+1-nn) = h(nn)
1753	100 CONTINUE
1754
1755	! normalize the sums to be unity
1756	sumh = 0
1757	do 150 n=1,NFILT
1758	sumh = sumh + H(n)
1759	150 continue
1760	print *,'SUMH =', sumh
1761
1762	do 200 n=1,NFILT
1763	H(n) = H(n)/sumh
1764	200 continue
1765	do 201 n=1,NFILT
1766	grid%hcoeff(n)=H(n)
1767	201 continue
1768	! print *,'HCOEFF(n) ', grid%hcoeff
1769	!
1770	END SUBROUTINE optfil
1771
1772
1773	SUBROUTINE MCCPAR (NFILT,JTYPE,NBANDS,LPRINT,LGRID, &
1774	EDGE,FX,WTX,DEVIAT,h )
1775
1776	! PROGRAM FOR THE DESIGN OF LINEAR PHASE FINITE IMPULSE
1777	! REPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM
1778	!
1779	!************************************************************
1780	!* Reference: McClellan, J.H., T.W. Parks and L.R.Rabiner, *
1781	!* 1973: A computer program for designing *
1782	!* optimum FIR linear phase digital filters. *
1783	!* IEEE Trans. on Audio and Electroacoustics, *
1784	!* Vol AU-21, No. 6, 506-526. *
1785	!************************************************************
1786	!
1787	! THREE TYPES OF FILTERS ARE INCLUDED -- BANDPASS FILTERS
1788	! DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS
1789	!
1790	!---------------------------------------------------------------
1791	!
1792	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
1793	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
1794	DIMENSION H(66)
1795	DIMENSION DES(1045),GRID(1045),WT(1045)
1796	DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
1797	DOUBLE PRECISION PI2,PI
1798	DOUBLE PRECISION AD,DEV,X,Y
1799	LOGICAL LPRINT
1800
1801	PI = 3.141592653589793
1802	PI2 = 6.283185307179586
1803
1804	! ......
1805
1806	NFMAX = 128
1807	100 CONTINUE
1808
1809	! PROGRAM INPUT SECTION
1810
1811	!CC READ(5,*) NFILT,JTYPE,NBANDS,JPRINT,LGRID
1812
1813	IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) THEN
1814	CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
1815	END IF
1816	IF(NBANDS.LE.0) NBANDS = 1
1817
1818	! ....
1819
1820	IF(LGRID.LE.0) LGRID = 16
1821	JB = 2*NBANDS
1822	!cc READ(5,*) (EDGE(J),J=1,JB)
1823	!cc READ(5,*) (FX(J),J=1,NBANDS)
1824	!cc READ(5,*) (WTX(J),J=1,NBANDS)
1825	IF(JTYPE.EQ.0) THEN
1826	CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
1827	END IF
1828	NEG = 1
1829	IF(JTYPE.EQ.1) NEG = 0
1830	NODD = NFILT/2
1831	NODD = NFILT-2*NODD
1832	NFCNS = NFILT/2
1833	IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS = NFCNS+1
1834
1835	! ...
1836
1837	GRID(1) = EDGE(1)
1838	DELF = LGRID*NFCNS
1839	DELF = 0.5/DELF
1840	IF(NEG.EQ.0) GOTO 135
1841	IF(EDGE(1).LT.DELF) GRID(1) = DELF
1842	135 CONTINUE
1843	J = 1
1844	L = 1
1845	LBAND = 1
1846	140 FUP = EDGE(L+1)
1847	145 TEMP = GRID(J)
1848
1849	! ....
1850
1851	DES(J) = EFF(TEMP,FX,WTX,LBAND,JTYPE)
1852	WT(J) = WATE(TEMP,FX,WTX,LBAND,JTYPE)
1853	J = J+1
1854	GRID(J) = TEMP+DELF
1855	IF(GRID(J).GT.FUP) GOTO 150
1856	GOTO 145
1857	150 GRID(J-1) = FUP
1858	DES(J-1) = EFF(FUP,FX,WTX,LBAND,JTYPE)
1859	WT(J-1) = WATE(FUP,FX,WTX,LBAND,JTYPE)
1860	LBAND = LBAND+1
1861	L = L+2
1862	IF(LBAND.GT.NBANDS) GOTO 160
1863	GRID(J) = EDGE(L)
1864	GOTO 140
1865	160 NGRID = J-1
1866	IF(NEG.NE.NODD) GOTO 165
1867	IF(GRID(NGRID).GT.(0.5-DELF)) NGRID = NGRID-1
1868	165 CONTINUE
1869
1870	! ......
1871
1872	IF(NEG) 170,170,180
1873	170 IF(NODD.EQ.1) GOTO 200
1874	DO 175 J=1,NGRID
1875	CHANGE = DCOS(PI*GRID(J))
1876	DES(J) = DES(J)/CHANGE
1877	WT(J) = WT(J)*CHANGE
1878	175 CONTINUE
1879	GOTO 200
1880	180 IF(NODD.EQ.1) GOTO 190
1881	DO 185 J = 1,NGRID
1882	CHANGE = DSIN(PI*GRID(J))
1883	DES(J) = DES(J)/CHANGE
1884	WT(J) = WT(J)*CHANGE
1885	185 CONTINUE
1886	GOTO 200
1887	190 DO 195 J =1,NGRID
1888	CHANGE = DSIN(PI2*GRID(J))
1889	DES(J) = DES(J)/CHANGE
1890	WT(J) = WT(J)*CHANGE
1891	195 CONTINUE
1892
1893	! ......
1894
1895	200 TEMP = FLOAT(NGRID-1)/FLOAT(NFCNS)
1896	DO 210 J = 1,NFCNS
1897	IEXT(J) = (J-1)*TEMP+1
1898	210 CONTINUE
1899	IEXT(NFCNS+1) = NGRID
1900	NM1 = NFCNS-1
1901	NZ = NFCNS+1
1902
1903	! CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM
1904
1905	CALL REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
1906
1907	! CALCULATE THE IMPULSE RESPONSE
1908
1909	IF(NEG) 300,300,320
1910	300 IF(NODD.EQ.0) GOTO 310
1911	DO 305 J=1,NM1
1912	H(J) = 0.5*ALPHA(NZ-J)
1913	305 CONTINUE
1914	H(NFCNS)=ALPHA(1)
1915	GOTO 350
1916	310 H(1) = 0.25*ALPHA(NFCNS)
1917	DO 315 J = 2,NM1
1918	H(J) = 0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
1919	315 CONTINUE
1920	H(NFCNS) = 0.5ALPHA(1)+0.25ALPHA(2)
1921	GOTO 350
1922	320 IF(NODD.EQ.0) GOTO 330
1923	H(1) = 0.25*ALPHA(NFCNS)
1924	H(2) = 0.25*ALPHA(NM1)
1925	DO 325 J = 3,NM1
1926	H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
1927	325 CONTINUE
1928	H(NFCNS) = 0.5ALPHA(1)-0.25ALPHA(3)
1929	H(NZ) = 0.0
1930	GOTO 350
1931	330 H(1) = 0.25*ALPHA(NFCNS)
1932	DO 335 J =2,NM1
1933	H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
1934	335 CONTINUE
1935	H(NFCNS) = 0.5ALPHA(1)-0.25ALPHA(2)
1936
1937	! PROGRAM OUTPUT SECTION
1938
1939	350 CONTINUE
1940	!
1941	IF(LPRINT) THEN
1942
1943	print , '***************************************************'
1944	print *, 'FINITE IMPULSE RESPONSE (FIR)'
1945	print *, 'LINEAR PHASE DIGITAL FILTER DESIGN'
1946	print *, 'REMEZ EXCHANGE ALGORITHM'
1947
1948	IF(JTYPE.EQ.1) WRITE(6,365)
1949	365 FORMAT(25X,'BANDPASS FILTER'/)
1950
1951	IF(JTYPE.EQ.2) WRITE(6,370)
1952	370 FORMAT(25X,'DIFFERENTIATOR '/)
1953
1954	IF(JTYPE.EQ.3) WRITE(6,375)
1955	375 FORMAT(25X,'HILBERT TRANSFORMER '/)
1956
1957	WRITE(6,378) NFILT
1958	378 FORMAT(15X,'FILTER LENGTH =',I3/)
1959
1960	WRITE(6,380)
1961	380 FORMAT(15X,'*** IMPULSE RESPONSE ***')
1962
1963	DO 381 J = 1,NFCNS
1964	K = NFILT+1-J
1965	IF(NEG.EQ.0) WRITE(6,382) J,H(J),K
1966	IF(NEG.EQ.1) WRITE(6,383) J,H(J),K
1967	381 CONTINUE
1968	382 FORMAT(20X,'H(',I3,') = ',E15.8,' = H(',I4,')')
1969	383 FORMAT(20X,'H(',I3,') = ',E15.8,' = -H(',I4,')')
1970
1971	IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ
1972	384 FORMAT(20X,'H(',I3,') = 0.0')
1973
1974	DO 450 K=1,NBANDS,4
1975	KUP = K+3
1976	IF(KUP.GT.NBANDS) KUP = NBANDS
1977	print *
1978	WRITE(6,385) (J,J=K,KUP)
1979	385 FORMAT(24X,4('BAND',I3,8X))
1980	WRITE(6,390) (EDGE(2*J-1),J=K,KUP)
1981	390 FORMAT(2X,'LOWER BAND EDGE',5F15.8)
1982	WRITE(6,395) (EDGE(2*J),J=K,KUP)
1983	395 FORMAT(2X,'UPPER BAND EDGE',5F15.8)
1984	IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
1985	400 FORMAT(2X,'DESIRED VALUE',2X,5F15.8)
1986	IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
1987	405 FORMAT(2X,'DESIRED SLOPE',2X,5F15.8)
1988	WRITE(6,410) (WTX(J),J=K,KUP)
1989	410 FORMAT(2X,'WEIGHTING',6X,5F15.8)
1990	DO 420 J = K,KUP
1991	DEVIAT(J) = DEV/WTX(J)
1992	420 CONTINUE
1993	WRITE(6,425) (DEVIAT(J),J=K,KUP)
1994	425 FORMAT(2X,'DEVIATION',6X,5F15.8)
1995	IF(JTYPE.NE.1) GOTO 450
1996	DO 430 J = K,KUP
1997	DEVIAT(J) = 20.0*ALOG10(DEVIAT(J))
1998	430 CONTINUE
1999	WRITE(6,435) (DEVIAT(J),J=K,KUP)
2000	435 FORMAT(2X,'DEVIATION IN DB',5F15.8)
2001	450 CONTINUE
2002	print *, 'EXTREMAL FREQUENCIES'
2003	WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
2004	455 FORMAT((2X,5F15.7))
2005	WRITE(6,460)
2006	460 FORMAT(1X,70(1H*))
2007	!
2008	ENDIF
2009	!
2010	!CC IF(NFILT.NE.0) GOTO 100 ! removal of re-run loop.
2011	!
2012	END SUBROUTINE mccpar
2013
2014
2015	FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE)
2016	DIMENSION FX(5),WTX(5)
2017	IF(JTYPE.EQ.2) GOTO 1
2018	EFF = FX(LBAND)
2019	RETURN
2020	1 EFF = FX(LBAND)*TEMP
2021	END FUNCTION eff
2022
2023
2024	FUNCTION WATE(TEMP,FX,WTX,LBAND,JTYPE)
2025	DIMENSION FX(5),WTX(5)
2026	IF(JTYPE.EQ.2) GOTO 1
2027	WATE = WTX(LBAND)
2028	RETURN
2029	1 IF(FX(LBAND).LT.0.0001) GOTO 2
2030	WATE = WTX(LBAND)/TEMP
2031	RETURN
2032	2 WATE = WTX(LBAND)
2033	END FUNCTION wate
2034
2035
2036	! SUBROUTINE ERROR
2037	!! WRITE(6,)' * ERROR IN INPUT DATA **'
2038	! CALL wrf_error_fatal (' ** ERROR IN INPUT DATA **' )
2039	! END SUBROUTINE error
2040
2041
2042	SUBROUTINE REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
2043	! THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
2044	! FOR THE WEIGHTED CHEBCHEV APPROXIMATION OF A CONTINUOUS
2045	! FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE
2046	! ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
2047	! DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
2048	! GRID, THE NUMBER OF COSINES, AND THE INITIAL GUESS OF THE
2049	! EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV
2050	! ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
2051	! FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
2052	! THE COEFFICIENTS OF THE BEST APPROXIMATION.
2053	!
2054	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
2055	DIMENSION EDGE(20)
2056	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
2057	DIMENSION DES(1045),GRID(1045),WT(1045)
2058	DIMENSION A(66),P(65),Q(65)
2059	DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
2060	DOUBLE PRECISION AD,DEV,X,Y
2061	DOUBLE PRECISION, EXTERNAL :: D, GEE
2062	!
2063	! THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
2064	!
2065	ITRMAX=25
2066	DEVL=-1.0
2067	NZ=NFCNS+1
2068	NZZ=NFCNS+2
2069	NITER=0
2070	100 CONTINUE
2071	IEXT(NZZ)=NGRID+1
2072	NITER=NITER+1
2073	IF(NITER.GT.ITRMAX) GO TO 400
2074	DO 110 J=1,NZ
2075	DTEMP=GRID(IEXT(J))
2076	DTEMP=DCOS(DTEMP*PI2)
2077	110 X(J)=DTEMP
2078	JET=(NFCNS-1)/15+1
2079	DO 120 J=1,NZ
2080	120 AD(J)=D(J,NZ,JET,X)
2081	DNUM=0.0
2082	DDEN=0.0
2083	K=1
2084	DO 130 J=1,NZ
2085	L=IEXT(J)
2086	DTEMP=AD(J)*DES(L)
2087	DNUM=DNUM+DTEMP
2088	DTEMP=K*AD(J)/WT(L)
2089	DDEN=DDEN+DTEMP
2090	130 K=-K
2091	DEV=DNUM/DDEN
2092	NU=1
2093	IF(DEV.GT.0.0) NU=-1
2094	DEV=-NU*DEV
2095	K=NU
2096	DO 140 J=1,NZ
2097	L=IEXT(J)
2098	DTEMP=K*DEV/WT(L)
2099	Y(J)=DES(L)+DTEMP
2100	140 K=-K
2101	IF(DEV.GE.DEVL) GO TO 150
2102	WRITE(6,) ' ***** FAILURE TO CONVERGE ********* '
2103	WRITE(6,*) ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR '
2104	WRITE(6,*) ' THE IMPULSE RESPONSE MAY BE CORRECT '
2105	WRITE(6,*) ' CHECK WITH A FREQUENCY RESPONSE '
2106	WRITE(6,) ' *************************************** '
2107	GO TO 400
2108	150 DEVL=DEV
2109	JCHNGE=0
2110	K1=IEXT(1)
2111	KNZ=IEXT(NZ)
2112	KLOW=0
2113	NUT=-NU
2114	J=1
2115	!
2116	! SEARCH FOR THE EXTERMAL FREQUENCIES OF THE BEST
2117	! APPROXIMATION.
2118
2119	200 IF(J.EQ.NZZ) YNZ=COMP
2120	IF(J.GE.NZZ) GO TO 300
2121	KUP=IEXT(J+1)
2122	L=IEXT(J)+1
2123	NUT=-NUT
2124	IF(J.EQ.2) Y1=COMP
2125	COMP=DEV
2126	IF(L.GE.KUP) GO TO 220
2127	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2128	ERR=(ERR-DES(L))*WT(L)
2129	DTEMP=NUT*ERR-COMP
2130	IF(DTEMP.LE.0.0) GO TO 220
2131	COMP=NUT*ERR
2132	210 L=L+1
2133	IF(L.GE.KUP) GO TO 215
2134	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2135	ERR=(ERR-DES(L))*WT(L)
2136	DTEMP=NUT*ERR-COMP
2137	IF(DTEMP.LE.0.0) GO TO 215
2138	COMP=NUT*ERR
2139	GO TO 210
2140	215 IEXT(J)=L-1
2141	J=J+1
2142	KLOW=L-1
2143	JCHNGE=JCHNGE+1
2144	GO TO 200
2145	220 L=L-1
2146	225 L=L-1
2147	IF(L.LE.KLOW) GO TO 250
2148	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2149	ERR=(ERR-DES(L))*WT(L)
2150	DTEMP=NUT*ERR-COMP
2151	IF(DTEMP.GT.0.0) GO TO 230
2152	IF(JCHNGE.LE.0) GO TO 225
2153	GO TO 260
2154	230 COMP=NUT*ERR
2155	235 L=L-1
2156	IF(L.LE.KLOW) GO TO 240
2157	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2158	ERR=(ERR-DES(L))*WT(L)
2159	DTEMP=NUT*ERR-COMP
2160	IF(DTEMP.LE.0.0) GO TO 240
2161	COMP=NUT*ERR
2162	GO TO 235
2163	240 KLOW=IEXT(J)
2164	IEXT(J)=L+1
2165	J=J+1
2166	JCHNGE=JCHNGE+1
2167	GO TO 200
2168	250 L=IEXT(J)+1
2169	IF(JCHNGE.GT.0) GO TO 215
2170	255 L=L+1
2171	IF(L.GE.KUP) GO TO 260
2172	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2173	ERR=(ERR-DES(L))*WT(L)
2174	DTEMP=NUT*ERR-COMP
2175	IF(DTEMP.LE.0.0) GO TO 255
2176	COMP=NUT*ERR
2177	GO TO 210
2178	260 KLOW=IEXT(J)
2179	J=J+1
2180	GO TO 200
2181	300 IF(J.GT.NZZ) GO TO 320
2182	IF(K1.GT.IEXT(1)) K1=IEXT(1)
2183	IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
2184	NUT1=NUT
2185	NUT=-NU
2186	L=0
2187	KUP=K1
2188	COMP=YNZ*(1.00001)
2189	LUCK=1
2190	310 L=L+1
2191	IF(L.GE.KUP) GO TO 315
2192	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2193	ERR=(ERR-DES(L))*WT(L)
2194	DTEMP=NUT*ERR-COMP
2195	IF(DTEMP.LE.0.0) GO TO 310
2196	COMP=NUT*ERR
2197	J=NZZ
2198	GO TO 210
2199	315 LUCK=6
2200	GO TO 325
2201	320 IF(LUCK.GT.9) GO TO 350
2202	IF(COMP.GT.Y1) Y1=COMP
2203	K1=IEXT(NZZ)
2204	325 L=NGRID+1
2205	KLOW=KNZ
2206	NUT=-NUT1
2207	COMP=Y1*(1.00001)
2208	330 L=L-1
2209	IF(L.LE.KLOW) GO TO 340
2210	ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
2211	ERR=(ERR-DES(L))*WT(L)
2212	DTEMP=NUT*ERR-COMP
2213	IF(DTEMP.LE.0.0) GO TO 330
2214	J=NZZ
2215	COMP=NUT*ERR
2216	LUCK=LUCK+10
2217	GO TO 235
2218	340 IF(LUCK.EQ.6) GO TO 370
2219	DO 345 J=1,NFCNS
2220	345 IEXT(NZZ-J)=IEXT(NZ-J)
2221	IEXT(1)=K1
2222	GO TO 100
2223	350 KN=IEXT(NZZ)
2224	DO 360 J=1,NFCNS
2225	360 IEXT(J)=IEXT(J+1)
2226	IEXT(NZ)=KN
2227	GO TO 100
2228	370 IF(JCHNGE.GT.0) GO TO 100
2229	!
2230	! CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
2231	! USING THE INVERSE DISCRETE FOURIER TRANSFORM.
2232	!
2233	400 CONTINUE
2234	NM1=NFCNS-1
2235	FSH=1.0E-06
2236	GTEMP=GRID(1)
2237	X(NZZ)=-2.0
2238	CN=2*NFCNS-1
2239	DELF=1.0/CN
2240	L=1
2241	KKK=0
2242	IF(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1
2243	IF(NFCNS.LE.3) KKK=1
2244	IF(KKK.EQ.1) GO TO 405
2245	DTEMP=DCOS(PI2*GRID(1))
2246	DNUM=DCOS(PI2*GRID(NGRID))
2247	AA=2.0/(DTEMP-DNUM)
2248	BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
2249	405 CONTINUE
2250	DO 430 J=1,NFCNS
2251	FT=(J-1)*DELF
2252	XT=DCOS(PI2*FT)
2253	IF(KKK.EQ.1) GO TO 410
2254	XT=(XT-BB)/AA
2255	! original : FT=ARCOS(XT)/PI2
2256	FT=ACOS(XT)/PI2
2257	410 XE=X(L)
2258	IF(XT.GT.XE) GO TO 420
2259	IF((XE-XT).LT.FSH) GO TO 415
2260	L=L+1
2261	GO TO 410
2262	415 A(J)=Y(L)
2263	GO TO 425
2264	420 IF((XT-XE).LT.FSH) GO TO 415
2265	GRID(1)=FT
2266	A(J)=GEE(1,NZ,GRID,PI2,X,Y,AD)
2267	425 CONTINUE
2268	IF(L.GT.1) L=L-1
2269	430 CONTINUE
2270	GRID(1)=GTEMP
2271	DDEN=PI2/CN
2272	DO 510 J=1,NFCNS
2273	DTEMP=0.0
2274	DNUM=(J-1)*DDEN
2275	IF(NM1.LT.1) GO TO 505
2276	DO 500 K=1,NM1
2277	500 DTEMP=DTEMP+A(K+1)DCOS(DNUMK)
2278	505 DTEMP=2.0*DTEMP+A(1)
2279	510 ALPHA(J)=DTEMP
2280	DO 550 J=2,NFCNS
2281	550 ALPHA(J)=2*ALPHA(J)/CN
2282	ALPHA(1)=ALPHA(1)/CN
2283	IF(KKK.EQ.1) GO TO 545
2284	P(1)=2.0ALPHA(NFCNS)BB+ALPHA(NM1)
2285	P(2)=2.0AAALPHA(NFCNS)
2286	Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
2287	DO 540 J=2,NM1
2288	IF(J.LT.NM1) GO TO 515
2289	AA=0.5*AA
2290	BB=0.5*BB
2291	515 CONTINUE
2292	P(J+1)=0.0
2293	DO 520 K=1,J
2294	A(K)=P(K)
2295	520 P(K)=2.0BBA(K)
2296	P(2)=P(2)+A(1)2.0AA
2297	JM1=J-1
2298	DO 525 K=1,JM1
2299	525 P(K)=P(K)+Q(K)+AA*A(K+1)
2300	JP1=J+1
2301	DO 530 K=3,JP1
2302	530 P(K)=P(K)+AA*A(K-1)
2303	IF(J.EQ.NM1) GO TO 540
2304	DO 535 K=1,J
2305	535 Q(K)=-A(K)
2306	Q(1)=Q(1)+ALPHA(NFCNS-1-J)
2307	540 CONTINUE
2308	DO 543 J=1,NFCNS
2309	543 ALPHA(J)=P(J)
2310	545 CONTINUE
2311	IF(NFCNS.GT.3) RETURN
2312	ALPHA(NFCNS+1)=0.0
2313	ALPHA(NFCNS+2)=0.0
2314	END SUBROUTINE remez
2315
2316	DOUBLE PRECISION FUNCTION D(K,N,M,X)
2317	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
2318	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
2319	DIMENSION DES(1045),GRID(1045),WT(1045)
2320	DOUBLE PRECISION AD,DEV,X,Y
2321	DOUBLE PRECISION Q
2322	DOUBLE PRECISION PI2
2323	D = 1.0
2324	Q = X(K)
2325	DO 3 L = 1,M
2326	DO 2 J = L,N,M
2327	IF(J-K) 1,2,1
2328	1 D = 2.0D(Q-X(J))
2329	2 CONTINUE
2330	3 CONTINUE
2331	D = 1.0/D
2332	END FUNCTION D
2333
2334
2335	DOUBLE PRECISION FUNCTION GEE(K,N,GRID,PI2,X,Y,AD)
2336	! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
2337	DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
2338	DIMENSION DES(1045),GRID(1045),WT(1045)
2339	DOUBLE PRECISION AD,DEV,X,Y
2340	DOUBLE PRECISION P,C,D,XF
2341	DOUBLE PRECISION PI2
2342	P = 0.0
2343	XF = GRID(K)
2344	XF = DCOS(PI2*XF)
2345	D = 0.0
2346	DO 1 J =1,N
2347	C = XF-X(J)
2348	C = AD(J)/C
2349	D = D+C
2350	P = P+C*Y(J)
2351	1 CONTINUE
2352	GEE = P/D
2353	END FUNCTION GEE
2354
2355	REAL FUNCTION RSLF(P,T)
2356
2357	IMPLICIT NONE
2358	REAL, INTENT(IN):: P, T
2359	REAL:: ESL,X
2360	REAL, PARAMETER:: C0= .611583699E03
2361	REAL, PARAMETER:: C1= .444606896E02
2362	REAL, PARAMETER:: C2= .143177157E01
2363	REAL, PARAMETER:: C3= .264224321E-1
2364	REAL, PARAMETER:: C4= .299291081E-3
2365	REAL, PARAMETER:: C5= .203154182E-5
2366	REAL, PARAMETER:: C6= .702620698E-8
2367	REAL, PARAMETER:: C7= .379534310E-11
2368	REAL, PARAMETER:: C8=-.321582393E-13
2369
2370	X=MAX(-80.,T-273.16)
2371
2372	! ESL=612.2EXP(17.67X/(T-29.65))
2373	ESL=C0+X(C1+X(C2+X(C3+X(C4+X(C5+X(C6+X(C7+XC8)))))))
2374	RSLF=.622*ESL/(P-ESL)
2375
2376	END FUNCTION RSLF
2377
2378	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2379	! DFI startfwd group of functions
2380	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2381
2382	SUBROUTINE wrf_dfi_startfwd_init ( )
2383
2384	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
2385	USE module_utility
2386
2387	IMPLICIT NONE
2388
2389	INTERFACE
2390	SUBROUTINE dfi_startfwd_init_recurse(grid)
2391	USE module_domain, ONLY : domain
2392	TYPE (domain), POINTER :: grid
2393	END SUBROUTINE dfi_startfwd_init_recurse
2394	END INTERFACE
2395
2396	! Now, setup all nests
2397
2398	CALL dfi_startfwd_init_recurse(head_grid)
2399
2400	CALL set_current_grid_ptr( head_grid )
2401
2402	END SUBROUTINE wrf_dfi_startfwd_init
2403
2404
2405	RECURSIVE SUBROUTINE dfi_startfwd_init_recurse(grid)
2406
2407	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
2408
2409	IMPLICIT NONE
2410
2411	INTERFACE
2412	SUBROUTINE dfi_startfwd_init(grid)
2413	USE module_domain, ONLY : domain
2414	TYPE (domain), POINTER :: grid
2415	END SUBROUTINE dfi_startfwd_init
2416	END INTERFACE
2417
2418	INTEGER :: kid
2419	TYPE (domain), POINTER :: grid
2420	TYPE (domain), POINTER :: grid_ptr
2421
2422	grid_ptr => grid
2423
2424	DO WHILE ( ASSOCIATED( grid_ptr ) )
2425	!
2426	! Assure that time-step is set back to positive
2427	! for this forward step.
2428	!
2429	grid_ptr%dt = abs(grid_ptr%dt)
2430	grid_ptr%time_step = abs(grid_ptr%time_step)
2431	CALL set_current_grid_ptr( grid_ptr )
2432	CALL dfi_startfwd_init( grid_ptr )
2433	DO kid = 1, max_nests
2434	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2435	CALL dfi_startfwd_init_recurse(grid_ptr%nests(kid)%ptr)
2436	ENDIF
2437	END DO
2438	grid_ptr => grid_ptr%sibling
2439	END DO
2440
2441	END SUBROUTINE dfi_startfwd_init_recurse
2442
2443
2444	SUBROUTINE dfi_startfwd_init ( grid )
2445
2446	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
2447	USE module_utility
2448	USE module_state_description
2449
2450	IMPLICIT NONE
2451
2452	TYPE (domain) , POINTER :: grid
2453	INTEGER rc
2454
2455	INTERFACE
2456	SUBROUTINE Setup_Timekeeping(grid)
2457	USE module_domain, ONLY : domain
2458	TYPE (domain), POINTER :: grid
2459	END SUBROUTINE Setup_Timekeeping
2460
2461	END INTERFACE
2462
2463	grid%dfi_stage = DFI_STARTFWD
2464
2465	#if (EM_CORE == 1)
2466	! No need for adaptive time-step
2467	CALL nl_set_use_adaptive_time_step( grid%id, .false. )
2468	#endif
2469
2470	CALL Setup_Timekeeping (grid)
2471	grid%start_subtime = domain_get_start_time ( head_grid )
2472	grid%stop_subtime = domain_get_stop_time ( head_grid )
2473
2474	CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
2475
2476	END SUBROUTINE dfi_startfwd_init
2477
2478	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2479	! DFI startbck group of functions
2480	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2481
2482	SUBROUTINE wrf_dfi_startbck_init ( )
2483
2484	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
2485	USE module_utility
2486
2487	IMPLICIT NONE
2488
2489	INTERFACE
2490	SUBROUTINE dfi_startbck_init_recurse(grid)
2491	USE module_domain, ONLY : domain
2492	TYPE (domain), POINTER :: grid
2493	END SUBROUTINE dfi_startbck_init_recurse
2494	END INTERFACE
2495
2496	! Now, setup all nests
2497
2498	CALL dfi_startbck_init_recurse(head_grid)
2499
2500	CALL set_current_grid_ptr( head_grid )
2501
2502	END SUBROUTINE wrf_dfi_startbck_init
2503
2504
2505	RECURSIVE SUBROUTINE dfi_startbck_init_recurse(grid)
2506
2507	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
2508
2509	IMPLICIT NONE
2510
2511	INTERFACE
2512	SUBROUTINE dfi_startbck_init(grid)
2513	USE module_domain, ONLY : domain
2514	TYPE (domain), POINTER :: grid
2515	END SUBROUTINE dfi_startbck_init
2516	END INTERFACE
2517
2518	INTEGER :: kid
2519	TYPE (domain), POINTER :: grid
2520	TYPE (domain), POINTER :: grid_ptr
2521
2522	grid_ptr => grid
2523
2524	DO WHILE ( ASSOCIATED( grid_ptr ) )
2525	!
2526	! Assure that time-step is set back to positive
2527	! for this forward step.
2528	!
2529	grid_ptr%dt = abs(grid_ptr%dt)
2530	grid_ptr%time_step = abs(grid_ptr%time_step)
2531	CALL set_current_grid_ptr( grid_ptr )
2532	CALL dfi_startbck_init( grid_ptr )
2533	DO kid = 1, max_nests
2534	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2535	CALL dfi_startbck_init_recurse(grid_ptr%nests(kid)%ptr)
2536	ENDIF
2537	END DO
2538	grid_ptr => grid_ptr%sibling
2539	END DO
2540
2541	END SUBROUTINE dfi_startbck_init_recurse
2542
2543
2544	SUBROUTINE dfi_startbck_init ( grid )
2545
2546	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
2547	USE module_utility
2548	USE module_state_description
2549
2550	IMPLICIT NONE
2551
2552	TYPE (domain) , POINTER :: grid
2553	INTEGER rc
2554
2555	INTERFACE
2556	SUBROUTINE Setup_Timekeeping(grid)
2557	USE module_domain, ONLY : domain
2558	TYPE (domain), POINTER :: grid
2559	END SUBROUTINE Setup_Timekeeping
2560
2561	END INTERFACE
2562
2563	grid%dfi_stage = DFI_STARTBCK
2564
2565	! set physics options to zero
2566	CALL nl_set_mp_physics( grid%id, 0 )
2567	CALL nl_set_ra_lw_physics( grid%id, 0 )
2568	CALL nl_set_ra_sw_physics( grid%id, 0 )
2569	CALL nl_set_sf_surface_physics( grid%id, 0 )
2570	CALL nl_set_sf_sfclay_physics( grid%id, 0 )
2571	CALL nl_set_sf_urban_physics( grid%id, 0 )
2572	CALL nl_set_bl_pbl_physics( grid%id, 0 )
2573	CALL nl_set_cu_physics( grid%id, 0 )
2574	CALL nl_set_damp_opt( grid%id, 0 )
2575	CALL nl_set_sst_update( grid%id, 0 )
2576	CALL nl_set_gwd_opt( grid%id, 0 )
2577	#if (EM_CORE == 1)
2578	CALL nl_set_diff_6th_opt( grid%id, 0 )
2579	CALL nl_set_use_adaptive_time_step( grid%id, .false. )
2580	#endif
2581	CALL nl_set_feedback( grid%id, 0 )
2582	#if (EM_CORE == 1)
2583	! set bc
2584	CALL nl_set_constant_bc( grid%id, head_grid%constant_bc)
2585	#endif
2586
2587	#ifdef WRF_CHEM
2588	! set chemistry option to zero
2589	CALL nl_set_chem_opt (grid%id, 0)
2590	CALL nl_set_aer_ra_feedback (grid%id, 0)
2591	CALL nl_set_io_form_auxinput5 (grid%id, 0)
2592	CALL nl_set_io_form_auxinput7 (grid%id, 0)
2593	CALL nl_set_io_form_auxinput8 (grid%id, 0)
2594	#endif
2595
2596	#if (EM_CORE == 1)
2597	! set diffusion to zero for backward integration
2598	CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
2599	CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
2600	IF ( grid%moist_adv_opt == 2 ) THEN
2601	CALL nl_set_moist_adv_opt( grid%id, 0)
2602	ENDIF
2603	#endif
2604
2605	!tgs need to call start_domain here to reset bc initialization for
2606	! negative dt, but only for outer domain.
2607	if (grid%id == 1) then
2608	CALL start_domain ( grid , .TRUE. )
2609	endif
2610
2611	! Call wrf_run to advance forward 1 step
2612
2613	! Negate time step
2614	CALL nl_set_time_step ( grid%id, -grid%time_step )
2615
2616	CALL Setup_Timekeeping (grid)
2617
2618	grid%start_subtime = domain_get_start_time ( grid )
2619	grid%stop_subtime = domain_get_stop_time ( grid )
2620
2621	CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
2622
2623	END SUBROUTINE dfi_startbck_init
2624
2625
2626	SUBROUTINE wrf_dfi_bck_init ( )
2627
2628	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
2629	USE module_utility
2630	USE module_state_description
2631
2632	IMPLICIT NONE
2633
2634	INTERFACE
2635	SUBROUTINE dfi_bck_init_recurse(grid)
2636	USE module_domain, ONLY : domain
2637	TYPE (domain), POINTER :: grid
2638	END SUBROUTINE dfi_bck_init_recurse
2639	END INTERFACE
2640
2641	! We can only call dfi_bck_init for the head_grid
2642	! since nests have not been instantiated at this point,
2643	! so, dfi_bck_init will need to be called for each
2644	! nest from integrate.
2645	CALL dfi_bck_init_recurse(head_grid)
2646
2647	END SUBROUTINE wrf_dfi_bck_init
2648
2649	RECURSIVE SUBROUTINE dfi_bck_init_recurse(grid)
2650
2651	USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
2652
2653	IMPLICIT NONE
2654
2655	INTERFACE
2656	SUBROUTINE dfi_bck_init(grid)
2657	USE module_domain, ONLY : domain
2658	TYPE (domain), POINTER :: grid
2659	END SUBROUTINE dfi_bck_init
2660	END INTERFACE
2661
2662	INTEGER :: kid
2663	TYPE (domain), POINTER :: grid
2664	TYPE (domain), POINTER :: grid_ptr
2665
2666	grid_ptr => grid
2667
2668	DO WHILE ( ASSOCIATED( grid_ptr ) )
2669	!
2670	! Assure that time-step is set back to positive
2671	! for this forward step.
2672	!
2673	grid_ptr%dt = abs(grid_ptr%dt)
2674	grid_ptr%time_step = abs(grid_ptr%time_step)
2675	CALL set_current_grid_ptr( grid_ptr )
2676	CALL dfi_bck_init( grid_ptr )
2677	DO kid = 1, max_nests
2678	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2679	CALL dfi_bck_init_recurse(grid_ptr%nests(kid)%ptr)
2680	ENDIF
2681	END DO
2682	grid_ptr => grid_ptr%sibling
2683	END DO
2684
2685	END SUBROUTINE dfi_bck_init_recurse
2686
2687	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2688	! DFI forward initialization group of functions
2689	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2690
2691	SUBROUTINE wrf_dfi_fwd_init ( )
2692
2693	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
2694	USE module_utility
2695
2696	IMPLICIT NONE
2697
2698	INTERFACE
2699	SUBROUTINE dfi_fwd_init_recurse(grid)
2700	USE module_domain, ONLY : domain
2701	TYPE (domain), POINTER :: grid
2702	END SUBROUTINE dfi_fwd_init_recurse
2703	END INTERFACE
2704
2705	! Now, setup all nests
2706
2707	CALL dfi_fwd_init_recurse(head_grid)
2708
2709	CALL set_current_grid_ptr( head_grid )
2710
2711	END SUBROUTINE wrf_dfi_fwd_init
2712
2713	RECURSIVE SUBROUTINE dfi_fwd_init_recurse(grid)
2714
2715	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
2716
2717	IMPLICIT NONE
2718
2719	INTERFACE
2720	SUBROUTINE dfi_fwd_init(grid)
2721	USE module_domain, ONLY : domain
2722	TYPE (domain), POINTER :: grid
2723	END SUBROUTINE dfi_fwd_init
2724	END INTERFACE
2725
2726	INTEGER :: kid
2727	TYPE (domain), POINTER :: grid
2728	TYPE (domain), POINTER :: grid_ptr
2729
2730	grid_ptr => grid
2731
2732	DO WHILE ( ASSOCIATED( grid_ptr ) )
2733	!
2734	! Assure that time-step is set back to positive
2735	! for this forward step.
2736	!
2737	grid_ptr%dt = abs(grid_ptr%dt)
2738	grid_ptr%time_step = abs(grid_ptr%time_step)
2739	CALL set_current_grid_ptr( grid_ptr )
2740	CALL dfi_fwd_init( grid_ptr )
2741	DO kid = 1, max_nests
2742	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2743	CALL dfi_fwd_init_recurse(grid_ptr%nests(kid)%ptr)
2744	ENDIF
2745	END DO
2746	grid_ptr => grid_ptr%sibling
2747	END DO
2748
2749	END SUBROUTINE dfi_fwd_init_recurse
2750
2751	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2752	! DFI forecast initialization group of functions
2753	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2754
2755	SUBROUTINE wrf_dfi_fst_init ( )
2756
2757	USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
2758	USE module_utility
2759
2760	IMPLICIT NONE
2761
2762	INTERFACE
2763	SUBROUTINE dfi_fst_init_recurse(grid)
2764	USE module_domain, ONLY : domain
2765	TYPE (domain), POINTER :: grid
2766	END SUBROUTINE dfi_fst_init_recurse
2767	END INTERFACE
2768
2769	! Now, setup all nests
2770
2771	CALL dfi_fst_init_recurse(head_grid)
2772
2773	CALL set_current_grid_ptr( head_grid )
2774
2775	END SUBROUTINE wrf_dfi_fst_init
2776
2777	RECURSIVE SUBROUTINE dfi_fst_init_recurse ( grid )
2778
2779	USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
2780
2781	IMPLICIT NONE
2782
2783	INTERFACE
2784	SUBROUTINE dfi_fst_init(grid)
2785	USE module_domain, ONLY : domain
2786	TYPE (domain), POINTER :: grid
2787	END SUBROUTINE dfi_fst_init
2788	END INTERFACE
2789
2790	INTEGER :: kid
2791	TYPE (domain), POINTER :: grid
2792	TYPE (domain), POINTER :: grid_ptr
2793
2794	grid_ptr => grid
2795
2796	DO WHILE ( ASSOCIATED( grid_ptr ) )
2797	!
2798	! Assure that time-step is set back to positive
2799	! for this forward step.
2800	!
2801	grid_ptr%dt = abs(grid_ptr%dt)
2802	grid_ptr%time_step = abs(grid_ptr%time_step)
2803	CALL set_current_grid_ptr( grid_ptr )
2804	CALL dfi_fst_init( grid_ptr )
2805	DO kid = 1, max_nests
2806	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2807	CALL dfi_fst_init_recurse(grid_ptr%nests(kid)%ptr)
2808	ENDIF
2809	END DO
2810	grid_ptr => grid_ptr%sibling
2811	END DO
2812
2813	END SUBROUTINE dfi_fst_init_recurse
2814
2815	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2816	! DFI write initialization group of functions
2817	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2818
2819	SUBROUTINE wrf_dfi_write_initialized_state( )
2820
2821	USE module_domain, ONLY : domain, head_grid
2822
2823	INTERFACE
2824	SUBROUTINE dfi_write_initialized_state_recurse(grid)
2825	USE module_domain, ONLY : domain
2826	TYPE (domain), POINTER :: grid
2827	END SUBROUTINE dfi_write_initialized_state_recurse
2828	END INTERFACE
2829
2830	! Now, setup all nests
2831
2832	CALL dfi_write_initialized_state_recurse(head_grid)
2833
2834	END SUBROUTINE wrf_dfi_write_initialized_state
2835
2836	RECURSIVE SUBROUTINE dfi_write_initialized_state_recurse( grid )
2837
2838	USE module_domain, ONLY : domain, max_nests
2839
2840	IMPLICIT NONE
2841
2842	INTERFACE
2843	SUBROUTINE dfi_write_initialized_state( grid )
2844	USE module_domain, ONLY : domain
2845	TYPE (domain), POINTER :: grid
2846	END SUBROUTINE dfi_write_initialized_state
2847	END INTERFACE
2848
2849	INTEGER :: kid
2850	TYPE (domain), POINTER :: grid
2851	TYPE (domain), POINTER :: grid_ptr
2852
2853	grid_ptr => grid
2854
2855	DO WHILE ( ASSOCIATED( grid_ptr ) )
2856	!
2857	! Assure that time-step is set back to positive
2858	! for this forward step.
2859	!
2860	CALL dfi_write_initialized_state( grid_ptr )
2861	DO kid = 1, max_nests
2862	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2863	CALL dfi_write_initialized_state_recurse(grid_ptr%nests(kid)%ptr)
2864	ENDIF
2865	END DO
2866	grid_ptr => grid_ptr%sibling
2867	END DO
2868
2869	END SUBROUTINE dfi_write_initialized_state_recurse
2870
2871
2872	RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid)
2873
2874	USE module_domain, ONLY : domain, max_nests, set_current_grid_ptr
2875
2876	IMPLICIT NONE
2877
2878	INTERFACE
2879	SUBROUTINE dfi_array_reset(grid)
2880	USE module_domain, ONLY : domain
2881	TYPE (domain), POINTER :: grid
2882	END SUBROUTINE dfi_array_reset
2883	END INTERFACE
2884
2885	INTEGER :: kid
2886	TYPE (domain), POINTER :: grid
2887	TYPE (domain), POINTER :: grid_ptr
2888
2889	grid_ptr => grid
2890
2891	DO WHILE ( ASSOCIATED( grid_ptr ) )
2892	CALL set_current_grid_ptr( grid_ptr )
2893	CALL dfi_array_reset( grid_ptr )
2894	DO kid = 1, max_nests
2895	IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
2896	CALL dfi_array_reset_recurse(grid_ptr%nests(kid)%ptr)
2897	ENDIF
2898	END DO
2899	grid_ptr => grid_ptr%sibling
2900	END DO
2901
2902	END SUBROUTINE dfi_array_reset_recurse
2903

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

Download in other formats:

Original Format