PageRenderTime 69ms CodeModel.GetById 24ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/share/dfi.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2909 lines | 1731 code | 567 blank | 611 comment | 17 complexity | cc24b016fea48a0ad33b12ada0258f2f MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1.    SUBROUTINE dfi_accumulate( grid )
    2. 

  2. 

  3.       USE module_domain, ONLY : domain
    4. 

  4. !      USE module_configure
    5. 

  5.       USE module_driver_constants
    6. 

  6.       USE module_machine
    7. 

  7. !      USE module_dm
    8. 

  8.       USE module_model_constants
    9. 

  9.       USE module_state_description
    10. 

  10. 

  11.       IMPLICIT NONE
    12. 

  12. 

  13.       REAL hn
    14. 

  14.       CHARACTER*80 mess
    15. 

  15. 

  16.       !  Input data.
    17. 

  17. 

  18.       TYPE(domain) , POINTER          :: grid
    19. 

  19. 

  20.       IF ( grid%dfi_opt .EQ. DFI_NODFI .OR. grid%dfi_stage .EQ. DFI_FST ) RETURN
    21. 

  21. 

  22. #if (EM_CORE == 1)
    23. 

  23. 

  24. 

  25.       hn = grid%hcoeff(grid%itimestep+1)
    26. 

  26. 

  27.       ! accumulate dynamic variables
    28. 

  28.        grid%dfi_mu(:,:)    = grid%dfi_mu(:,:)    + grid%mu_2(:,:)    * hn
    29. 

  29.        grid%dfi_u(:,:,:)   = grid%dfi_u(:,:,:)   + grid%u_2(:,:,:)   * hn
    30. 

  30.        grid%dfi_v(:,:,:)   = grid%dfi_v(:,:,:)   + grid%v_2(:,:,:)   * hn
    31. 

  31.        grid%dfi_w(:,:,:)   = grid%dfi_w(:,:,:)   + grid%w_2(:,:,:)   * hn
    32. 

  32.        grid%dfi_ww(:,:,:)  = grid%dfi_ww(:,:,:)  + grid%ww(:,:,:)    * hn
    33. 

  33.        grid%dfi_t(:,:,:)   = grid%dfi_t(:,:,:)   + grid%t_2(:,:,:)   * hn
    34. 

  34.        grid%dfi_phb(:,:,:) = grid%dfi_phb(:,:,:) + grid%phb(:,:,:)   * hn
    35. 

  35.        grid%dfi_ph0(:,:,:) = grid%dfi_ph0(:,:,:) + grid%ph0(:,:,:)   * hn
    36. 

  36.        grid%dfi_php(:,:,:) = grid%dfi_php(:,:,:) + grid%php(:,:,:)   * hn
    37. 

  37.        grid%dfi_p(:,:,:)   = grid%dfi_p(:,:,:)   + grid%p(:,:,:)     * hn
    38. 

  38.        grid%dfi_ph(:,:,:)  = grid%dfi_ph(:,:,:)  + grid%ph_2(:,:,:)  * hn
    39. 

  39.        grid%dfi_tke(:,:,:) = grid%dfi_tke(:,:,:) + grid%tke_2(:,:,:) * hn
    40. 

  40.        grid%dfi_al(:,:,:)  = grid%dfi_al(:,:,:)  + grid%al(:,:,:)    * hn
    41. 

  41.        grid%dfi_alt(:,:,:) = grid%dfi_alt(:,:,:) + grid%alt(:,:,:)   * hn
    42. 

  42.        grid%dfi_pb(:,:,:)  = grid%dfi_pb(:,:,:)  + grid%pb(:,:,:)    * hn
    43. 

  43.       ! neg. check on hydrometeor and scalar variables
    44. 

  44.        grid%moist(:,:,:,:) = max(0.,grid%moist(:,:,:,:))
    45. 

  45.        grid%dfi_scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:))
    46. 

  46. #if ( WRF_DFI_RADAR == 1 )
    47. 

  47.      IF ( grid%dfi_radar .EQ. 0 ) then  
    48. 

  48.        grid%dfi_moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  + grid%moist(:,:,:,:)  * hn
    49. 

  49.      ELSE
    50. 

  50.        grid%dfi_moist(:,:,:,P_QV) = grid%dfi_moist(:,:,:,P_QV) + grid%moist(:,:,:,P_QV) * hn
    51. 

  51.      ENDIF
    52. 

  52. #else
    53. 

  53.        grid%dfi_moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  + grid%moist(:,:,:,:)  * hn
    54. 

  54. #endif
    55. 

  55.        grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
    56. 

  56. 

  57.       ! accumulate DFI coefficient
    58. 

  58.       grid%hcoeff_tot = grid%hcoeff_tot + hn
    59. 

  59. #endif
    60. 

  60. 

  61. #if (NMM_CORE == 1)
    62. 

  62.       hn = grid%hcoeff(grid%ntsd+1)
    63. 

  63.       grid%dfi_pd(:,:)    = grid%dfi_pd(:,:)    + grid%pd(:,:)    * hn
    64. 

  64.       grid%dfi_pint(:,:,:)   = grid%dfi_pint(:,:,:)   + grid%pint(:,:,:)   * hn
    65. 

  65.       grid%dfi_dwdt(:,:,:)   = grid%dfi_dwdt(:,:,:)   + grid%dwdt(:,:,:)   * hn
    66. 

  66.       grid%dfi_t(:,:,:)   = grid%dfi_t(:,:,:)   + grid%t(:,:,:)   * hn
    67. 

  67.       grid%dfi_q(:,:,:)   = grid%dfi_q(:,:,:)   + grid%q(:,:,:)   * hn
    68. 

  68.       grid%dfi_q2(:,:,:)  = grid%dfi_q2(:,:,:)  + grid%q2(:,:,:)  * hn
    69. 

  69.       grid%dfi_cwm(:,:,:) = grid%dfi_cwm(:,:,:) + grid%cwm(:,:,:) * hn
    70. 

  70.       grid%dfi_u(:,:,:)   = grid%dfi_u(:,:,:)   + grid%u(:,:,:)   * hn
    71. 

  71.       grid%dfi_v(:,:,:)   = grid%dfi_v(:,:,:)   + grid%v(:,:,:)   * hn
    72. 

  72.       grid%dfi_moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  + grid%moist(:,:,:,:)  * hn
    73. 

  73.       grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
    74. 

  74.       ! accumulate DFI coefficient
    75. 

  75.       grid%hcoeff_tot = grid%hcoeff_tot + hn
    76. 

  76.       write(mess,*) 'grid%hcoeff_tot: ', grid%hcoeff_tot
    77. 

  77.       call wrf_message(mess)
    78. 

  78. #endif
    79. 

  79. 

  80.    END SUBROUTINE dfi_accumulate
    81. 

  81. 

  82.    SUBROUTINE dfi_bck_init ( grid )
    83. 

  83. 

  84.       USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
    85. 

  85.       USE module_utility
    86. 

  86.       USE module_state_description
    87. 

  87. 

  88.       IMPLICIT NONE
    89. 

  89. 

  90.       TYPE (domain) , POINTER                 :: grid
    91. 

  91.       INTEGER rc
    92. 

  92.       CHARACTER*80 mess
    93. 

  93. 

  94.       INTERFACE
    95. 

  95.          SUBROUTINE Setup_Timekeeping(grid)
    96. 

  96.             USE module_domain, ONLY : domain
    97. 

  97.             TYPE (domain), POINTER :: grid
    98. 

  98.          END SUBROUTINE Setup_Timekeeping
    99. 

  99. 

  100.          SUBROUTINE dfi_save_arrays(grid)
    101. 

  101.             USE module_domain, ONLY : domain
    102. 

  102.             TYPE (domain), POINTER :: grid
    103. 

  103.          END SUBROUTINE dfi_save_arrays
    104. 

  104. 

  105.          SUBROUTINE dfi_clear_accumulation(grid)
    106. 

  106.             USE module_domain, ONLY : domain
    107. 

  107.             TYPE (domain), POINTER :: grid
    108. 

  108.          END SUBROUTINE dfi_clear_accumulation
    109. 

  109. 

  110.          SUBROUTINE optfil_driver(grid)
    111. 

  111.             USE module_domain, ONLY : domain
    112. 

  112.             TYPE (domain), POINTER :: grid
    113. 

  113.          END SUBROUTINE optfil_driver
    114. 

  114. 

  115.          SUBROUTINE start_domain(grid, allowed_to_read)
    116. 

  116.             USE module_domain, ONLY : domain
    117. 

  117.             TYPE (domain)          :: grid
    118. 

  118.             LOGICAL, INTENT(IN)    :: allowed_to_read
    119. 

  119.          END SUBROUTINE start_domain
    120. 

  120.       END INTERFACE
    121. 

  121. 

  122.       grid%dfi_stage = DFI_BCK
    123. 

  123. 

  124.       ! Negate time step
    125. 

  125.       IF ( grid%time_step_dfi .gt. 0 ) THEN
    126. 

  126.         CALL nl_set_time_step ( 1, -grid%time_step_dfi )
    127. 

  127.       ELSE
    128. 

  128.         CALL nl_set_time_step ( 1, -grid%time_step )
    129. 

  129.         CALL nl_set_time_step_fract_num ( 1, -grid%time_step_fract_num )
    130. 

  130.       ENDIF
    131. 

  131. 

  132.       CALL Setup_Timekeeping (grid)
    133. 

  133. 

  134.       !tgs 7apr11 - need to call start_domain here to reset bc initialization for negative dt
    135. 

  135.       CALL start_domain ( grid , .TRUE. )
    136. 

  136.       !tgs 7apr11 - save arrays should be done after start_domain to get correct grid%p field
    137. 

  137.       CALL dfi_save_arrays ( grid )
    138. 

  138. 

  139.       ! set physics options to zero
    140. 

  140.       CALL nl_set_mp_physics( grid%id, 0 )
    141. 

  141.       CALL nl_set_ra_lw_physics( grid%id, 0 )
    142. 

  142.       CALL nl_set_ra_sw_physics( grid%id, 0 )
    143. 

  143.       CALL nl_set_sf_surface_physics( grid%id, 0 )
    144. 

  144.       CALL nl_set_sf_sfclay_physics( grid%id, 0 )
    145. 

  145.       CALL nl_set_sf_urban_physics( grid%id, 0 )
    146. 

  146.       CALL nl_set_bl_pbl_physics( grid%id, 0 )
    147. 

  147.       CALL nl_set_cu_physics( grid%id, 0 )
    148. 

  148.       CALL nl_set_damp_opt( grid%id, 0 )
    149. 

  149.       CALL nl_set_sst_update( grid%id, 0 )
    150. 

  150.       CALL nl_set_fractional_seaice( grid%id, 0 )
    151. 

  151.       CALL nl_set_gwd_opt( grid%id, 0 )
    152. 

  152.       CALL nl_set_feedback( grid%id, 0 )
    153. 

  153.       ! set bc
    154. 

  154. #if (EM_CORE == 1)
    155. 

  155.       CALL nl_set_diff_6th_opt( grid%id, 0 )
    156. 

  156.       CALL nl_set_constant_bc(1, grid%constant_bc)
    157. 

  157.       CALL nl_set_use_adaptive_time_step( grid%id, .false. )
    158. 

  158. #endif
    159. 

  159. 

  160. #ifdef WRF_CHEM
    161. 

  161.       ! set chemistry option to zero
    162. 

  162.       CALL nl_set_chem_opt (grid%id, 0)
    163. 

  163.       CALL nl_set_aer_ra_feedback (grid%id, 0)
    164. 

  164.       CALL nl_set_io_form_auxinput5 (grid%id, 0)
    165. 

  165.       CALL nl_set_io_form_auxinput7 (grid%id, 0)
    166. 

  166.       CALL nl_set_io_form_auxinput8 (grid%id, 0)
    167. 

  167. #endif
    168. 

  168. 

  169.       ! set diffusion to zero for backward integration
    170. 

  170. 

  171. #if (EM_CORE == 1)
    172. 

  172.       CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
    173. 

  173.       CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
    174. 

  174.       IF ( grid%moist_adv_opt == 2 ) THEN
    175. 

  175.          CALL nl_set_moist_adv_opt( grid%id, 0)
    176. 

  176.       ENDIF
    177. 

  177. #endif
    178. 

  178. 

  179. #if (NMM_CORE == 1)
    180. 

  180. 

  181. !
    182. 

  182. ! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
    183. 

  183. !
    184. 

  184.       CALL nl_get_time_step( grid%id, grid%time_step )
    185. 

  185.       if (grid%time_step .lt. 0) then
    186. 

  186. !       DT   =-DT
    187. 

  187. 

  188.         write(mess,*) 'changing signs for backward integration'
    189. 

  189.         call wrf_message(mess)
    190. 

  190.         grid%CPGFV = -grid%CPGFV
    191. 

  191.         grid%EN    = -grid%EN
    192. 

  192.         grid%ENT   = -grid%ENT
    193. 

  193.         grid%F4D   = -grid%F4D
    194. 

  194.         grid%F4Q   = -grid%F4Q
    195. 

  195.         grid%EF4T  = -grid%EF4T
    196. 

  196.   
    197. 

  197.         grid%EM(:) = -grid%EM(:)
    198. 

  198.         grid%EMT(:) = -grid%EMT(:)
    199. 

  199.         grid%F4Q2(:) = -grid%F4Q2(:)
    200. 

  200.   
    201. 

  201.         grid%WPDAR(:,:)= -grid%WPDAR(:,:)
    202. 

  202.         grid%CPGFU(:,:)= -grid%CPGFU(:,:)
    203. 

  203.         grid%CURV(:,:)= -grid%CURV(:,:)
    204. 

  204.         grid%FCP(:,:)= -grid%FCP(:,:)
    205. 

  205.         grid%FAD(:,:)= -grid%FAD(:,:)
    206. 

  206.         grid%F(:,:)= -grid%F(:,:)
    207. 

  207.       endif
    208. 

  208. #endif
    209. 

  209. 

  210.       grid%start_subtime = domain_get_start_time ( grid )
    211. 

  211.       grid%stop_subtime = domain_get_stop_time ( grid )
    212. 

  212. 

  213.       CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
    214. 

  214. 

  215. !tgs 7apr11 - this call has been moved up
    216. 

  216. !     CALL dfi_save_arrays ( grid )
    217. 

  217.       CALL dfi_clear_accumulation( grid )
    218. 

  218.       CALL optfil_driver(grid)
    219. 

  219. 

  220.       !tgs need to call start_domain here to reset bc initialization for negative dt
    221. 

  221.       CALL start_domain ( grid , .TRUE. )
    222. 

  222. 

  223.    END SUBROUTINE dfi_bck_init
    224. 

  224. 

  225. 

  226.    SUBROUTINE dfi_fwd_init ( grid )
    227. 

  227. 

  228.       USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
    229. 

  229.       USE module_utility, ONLY : WRFU_ClockSet
    230. 

  230.       USE module_state_description
    231. 

  231. 

  232.       IMPLICIT NONE
    233. 

  233. 

  234.       TYPE (domain) , POINTER                 :: grid
    235. 

  235.       INTEGER rc
    236. 

  236.       CHARACTER*80 mess
    237. 

  237. 

  238.       INTERFACE
    239. 

  239.          SUBROUTINE Setup_Timekeeping(grid)
    240. 

  240.             USE module_domain, ONLY : domain
    241. 

  241.             TYPE (domain), POINTER :: grid
    242. 

  242.          END SUBROUTINE Setup_Timekeeping
    243. 

  243. 

  244.          SUBROUTINE dfi_save_arrays(grid)
    245. 

  245.             USE module_domain, ONLY : domain
    246. 

  246.             TYPE (domain), POINTER :: grid
    247. 

  247.          END SUBROUTINE dfi_save_arrays
    248. 

  248. 

  249.          SUBROUTINE dfi_clear_accumulation(grid)
    250. 

  250.             USE module_domain, ONLY : domain
    251. 

  251.             TYPE (domain), POINTER :: grid
    252. 

  252.          END SUBROUTINE dfi_clear_accumulation
    253. 

  253. 

  254.          SUBROUTINE optfil_driver(grid)
    255. 

  255.             USE module_domain, ONLY : domain
    256. 

  256.             TYPE (domain), POINTER :: grid
    257. 

  257.          END SUBROUTINE optfil_driver
    258. 

  258. 

  259.          SUBROUTINE start_domain(grid, allowed_to_read)
    260. 

  260.             USE module_domain, ONLY : domain
    261. 

  261.             TYPE (domain)          :: grid
    262. 

  262.             LOGICAL, INTENT(IN)    :: allowed_to_read
    263. 

  263.          END SUBROUTINE start_domain
    264. 

  264.       END INTERFACE
    265. 

  265. 

  266.       grid%dfi_stage = DFI_FWD
    267. 

  267. 

  268.       ! for Setup_Timekeeping to use when setting the clock
    269. 

  269. 

  270.       IF ( grid%time_step_dfi .gt. 0 ) THEN
    271. 

  271.         CALL nl_set_time_step ( grid%id, grid%time_step_dfi )
    272. 

  272.       ELSE
    273. 

  273. 

  274.       CALL nl_get_time_step( grid%id, grid%time_step )
    275. 

  275.         CALL nl_get_time_step_fract_num( grid%id, grid%time_step_fract_num )
    276. 

  276. 

  277.       grid%time_step = abs(grid%time_step)
    278. 

  278.       CALL nl_set_time_step( grid%id, grid%time_step )
    279. 

  279. 

  280.         grid%time_step_fract_num = abs(grid%time_step_fract_num)
    281. 

  281.         CALL nl_set_time_step_fract_num( grid%id, grid%time_step_fract_num )
    282. 

  282. 

  283.       ENDIF
    284. 

  284. 

  285.       grid%itimestep=0
    286. 

  286.       grid%xtime=0.
    287. 

  287. 

  288.       ! reset physics options to normal
    289. 

  289.       CALL nl_set_mp_physics( grid%id, grid%mp_physics)
    290. 

  290.       CALL nl_set_ra_lw_physics( grid%id, grid%ra_lw_physics)
    291. 

  291.       CALL nl_set_ra_sw_physics( grid%id, grid%ra_sw_physics)
    292. 

  292.       CALL nl_set_sf_surface_physics( grid%id, grid%sf_surface_physics)
    293. 

  293.       CALL nl_set_sf_sfclay_physics( grid%id, grid%sf_sfclay_physics)
    294. 

  294.       CALL nl_set_sf_urban_physics( grid%id, grid%sf_urban_physics)
    295. 

  295.       CALL nl_set_bl_pbl_physics( grid%id, grid%bl_pbl_physics)
    296. 

  296.       CALL nl_set_cu_physics( grid%id, grid%cu_physics)
    297. 

  297.       CALL nl_set_damp_opt( grid%id, grid%damp_opt)
    298. 

  298.       CALL nl_set_sst_update( grid%id, 0)
    299. 

  299.       CALL nl_set_fractional_seaice( grid%id, grid%fractional_seaice)
    300. 

  300.       CALL nl_set_gwd_opt( grid%id, grid%gwd_opt)
    301. 

  301.       CALL nl_set_feedback( grid%id, 0 )
    302. 

  302. #if (EM_CORE == 1) 
    303. 

  303.       CALL nl_set_diff_6th_opt( grid%id, grid%diff_6th_opt)
    304. 

  304.       CALL nl_set_use_adaptive_time_step( grid%id, .false. )
    305. 

  305.       ! set bc
    306. 

  306.       CALL nl_set_constant_bc(grid%id, grid%constant_bc)
    307. 

  307. #endif
    308. 

  308. 

  309. #if (NMM_CORE == 1)
    310. 

  310. !
    311. 

  311. ! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
    312. 

  312. !
    313. 

  313. !mptest      CALL nl_get_time_step( grid%id, grid%time_step )
    314. 

  314. 

  315. 

  316. !!! here need to key off something else being the "wrong" sign
    317. 

  317.       if (grid%cpgfv .gt. 0) then
    318. 

  318. 

  319.         write(mess,*) 'changing signs for forward integration'
    320. 

  320.         call wrf_message(mess)
    321. 

  321.         grid%CPGFV = -grid%CPGFV
    322. 

  322.         grid%EN    = -grid%EN
    323. 

  323.         grid%ENT   = -grid%ENT
    324. 

  324.         grid%F4D   = -grid%F4D
    325. 

  325.         grid%F4Q   = -grid%F4Q
    326. 

  326.         grid%EF4T  = -grid%EF4T
    327. 

  327.   
    328. 

  328.         grid%EM(:) = -grid%EM(:)
    329. 

  329.         grid%EMT(:) = -grid%EMT(:)
    330. 

  330.         grid%F4Q2(:) = -grid%F4Q2(:)
    331. 

  331.   
    332. 

  332.         grid%WPDAR(:,:)= -grid%WPDAR(:,:)
    333. 

  333.         grid%CPGFU(:,:)= -grid%CPGFU(:,:)
    334. 

  334.         grid%CURV(:,:)= -grid%CURV(:,:)
    335. 

  335.         grid%FCP(:,:)= -grid%FCP(:,:)
    336. 

  336.         grid%FAD(:,:)= -grid%FAD(:,:)
    337. 

  337.         grid%F(:,:)= -grid%F(:,:)
    338. 

  338.       endif
    339. 

  339. #endif
    340. 

  340. 

  341. 

  342. !#ifdef WRF_CHEM
    343. 

  343. !      ! reset chem option to normal
    344. 

  344. !      CALL nl_set_chem_opt( grid%id, grid%chem_opt)
    345. 

  345. !      CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
    346. 

  346. !#endif
    347. 

  347. 

  348. #if (EM_CORE == 1)
    349. 

  349.       ! reset km_opt to norlmal
    350. 

  350.       CALL nl_set_km_opt( grid%id, grid%km_opt)
    351. 

  351.       CALL nl_set_moist_adv_opt( grid%id, grid%moist_adv_opt)
    352. 

  352. #endif
    353. 

  353. 

  354. 

  355.       CALL Setup_Timekeeping (grid)
    356. 

  356.       grid%start_subtime = domain_get_start_time ( head_grid )
    357. 

  357.       grid%stop_subtime = domain_get_stop_time ( head_grid )
    358. 

  358. 

  359.       CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
    360. 

  360. 

  361.       IF ( grid%dfi_opt .EQ. DFI_DFL ) THEN
    362. 

  362.          CALL dfi_save_arrays ( grid )
    363. 

  363.       END IF
    364. 

  364.       CALL dfi_clear_accumulation( grid )
    365. 

  365.       CALL optfil_driver(grid)
    366. 

  366. 

  367.       !tgs need to call it here to reset bc initialization for positive time_step
    368. 

  368.       CALL start_domain ( grid , .TRUE. )
    369. 

  369. 

  370.    END SUBROUTINE dfi_fwd_init
    371. 

  371. 

  372.    SUBROUTINE dfi_fst_init ( grid )
    373. 

  373. 

  374.       USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time
    375. 

  375.       USE module_state_description
    376. 

  376. 

  377.       IMPLICIT NONE
    378. 

  378. 

  379.       TYPE (domain) , POINTER :: grid
    380. 

  380.       CHARACTER (LEN=80)      :: wrf_error_message
    381. 

  381. 

  382.       INTERFACE
    383. 

  383.          SUBROUTINE Setup_Timekeeping(grid)
    384. 

  384.             USE module_domain, ONLY : domain
    385. 

  385.             TYPE (domain), POINTER :: grid
    386. 

  386.          END SUBROUTINE Setup_Timekeeping
    387. 

  387. 

  388.          SUBROUTINE dfi_save_arrays(grid)
    389. 

  389.             USE module_domain, ONLY : domain
    390. 

  390.             TYPE (domain), POINTER :: grid
    391. 

  391.          END SUBROUTINE dfi_save_arrays
    392. 

  392. 

  393.          SUBROUTINE dfi_clear_accumulation(grid)
    394. 

  394.             USE module_domain, ONLY : domain
    395. 

  395.             TYPE (domain), POINTER :: grid
    396. 

  396.          END SUBROUTINE dfi_clear_accumulation
    397. 

  397. 

  398.          SUBROUTINE optfil_driver(grid)
    399. 

  399.             USE module_domain, ONLY : domain
    400. 

  400.             TYPE (domain), POINTER :: grid
    401. 

  401.          END SUBROUTINE optfil_driver
    402. 

  402. 

  403.          SUBROUTINE start_domain(grid, allowed_to_read)
    404. 

  404.             USE module_domain, ONLY : domain
    405. 

  405.             TYPE (domain)          :: grid
    406. 

  406.             LOGICAL, INTENT(IN)    :: allowed_to_read
    407. 

  407.          END SUBROUTINE start_domain
    408. 

  408.       END INTERFACE
    409. 

  409. 

  410.       grid%dfi_stage = DFI_FST 
    411. 

  411. 

  412.      ! reset time_step to normal and adaptive_time_step
    413. 

  413.       CALL nl_set_time_step( grid%id, grid%time_step )
    414. 

  414. 

  415.       grid%itimestep=0
    416. 

  416.       grid%xtime=0.         ! BUG: This will probably not work for all DFI options
    417. 

  417. ! only use adaptive time stepping for forecast
    418. 

  418. #if (EM_CORE == 1)
    419. 

  419.       if (grid%id == 1) then
    420. 

  420.          CALL nl_set_use_adaptive_time_step( 1, grid%use_adaptive_time_step )
    421. 

  421.       endif
    422. 

  422.       CALL nl_set_sst_update( grid%id, grid%sst_update)
    423. 

  423.       ! reset to normal bc
    424. 

  424.       CALL nl_set_constant_bc(grid%id, .false.)
    425. 

  425. #endif
    426. 

  426.       CALL nl_set_feedback( grid%id, grid%feedback )
    427. 

  427. 

  428. #ifdef WRF_CHEM
    429. 

  429.       ! reset chem option to normal
    430. 

  430.       CALL nl_set_chem_opt( grid%id, grid%chem_opt)
    431. 

  431.       CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
    432. 

  432.       CALL nl_set_io_form_auxinput5 (grid%id, grid%io_form_auxinput5)
    433. 

  433.       CALL nl_set_io_form_auxinput7 (grid%id, grid%io_form_auxinput7)
    434. 

  434.       CALL nl_set_io_form_auxinput8 (grid%id, grid%io_form_auxinput8)
    435. 

  435. #endif
    436. 

  436. 

  437. 

  438.       CALL Setup_Timekeeping (grid)
    439. 

  439.       grid%start_subtime = domain_get_start_time ( grid )
    440. 

  440.       grid%stop_subtime = domain_get_stop_time ( grid )
    441. 

  441. 

  442.       CALL start_domain ( grid , .TRUE. )
    443. 

  443. 

  444.    END SUBROUTINE dfi_fst_init
    445. 

  445. 

  446. 

  447.    SUBROUTINE dfi_write_initialized_state( grid )
    448. 

  448. 

  449.       ! Driver layer
    450. 

  450.       USE module_domain, ONLY : domain, head_grid
    451. 

  451.       USE module_io_domain
    452. 

  452.       USE module_configure, ONLY : grid_config_rec_type, model_config_rec
    453. 

  453. 

  454.       IMPLICIT NONE
    455. 

  455. 

  456.       TYPE (domain) , POINTER :: grid
    457. 

  457.       INTEGER             :: fid, ierr
    458. 

  458.       CHARACTER (LEN=80)  :: wrf_error_message
    459. 

  459.       CHARACTER (LEN=80)  :: rstname
    460. 

  460.       CHARACTER (LEN=132) :: message
    461. 

  461. 

  462.       TYPE (grid_config_rec_type) :: config_flags
    463. 

  463. 

  464.       CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
    465. 

  465. 

  466.       WRITE (wrf_err_message,'(A,I4)') 'Writing out initialized model state'
    467. 

  467.       CALL wrf_message(TRIM(wrf_err_message))
    468. 

  468. 

  469.       WRITE (rstname,'(A,I2.2)')'wrfinput_initialized_d',grid%id
    470. 

  470.       CALL open_w_dataset ( fid, TRIM(rstname), grid, config_flags, output_input, "DATASET=INPUT", ierr )
    471. 

  471.       IF ( ierr .NE. 0 ) THEN
    472. 

  472.          WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
    473. 

  473.          CALL WRF_ERROR_FATAL ( message )
    474. 

  474.       END IF
    475. 

  475.       CALL output_input ( fid,   grid, config_flags, ierr )
    476. 

  476.       CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
    477. 

  477. 

  478.    END SUBROUTINE dfi_write_initialized_state
    479. 

  479. 

  480. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    481. 

  481. ! DFI array reset group of functions
    482. 

  482. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    483. 

  483. 

  484.    SUBROUTINE wrf_dfi_array_reset ( )
    485. 

  485. 

  486.       USE module_domain, ONLY : domain, head_grid, set_current_grid_ptr
    487. 

  487. 

  488.       IMPLICIT NONE
    489. 

  489. 

  490.       INTERFACE
    491. 

  491.          RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid)
    492. 

  492.             USE module_domain, ONLY : domain
    493. 

  493.             TYPE (domain), POINTER :: grid
    494. 

  494.          END SUBROUTINE dfi_array_reset_recurse
    495. 

  495.       END INTERFACE
    496. 

  496. 

  497.       ! Copy filtered arrays back into state arrays in grid structure, and 
    498. 

  498.       !   restore original land surface fields
    499. 

  499. 

  500.       CALL dfi_array_reset_recurse(head_grid)
    501. 

  501. 

  502.       CALL set_current_grid_ptr( head_grid )
    503. 

  503. 

  504.    END SUBROUTINE wrf_dfi_array_reset
    505. 

  505. 

  506.    SUBROUTINE dfi_array_reset( grid )
    507. 

  507. 

  508.       USE module_domain, ONLY : domain
    509. 

  509. !      USE module_configure
    510. 

  510. !      USE module_driver_constants
    511. 

  511. !      USE module_machine
    512. 

  512. !      USE module_dm
    513. 

  513.       USE module_model_constants
    514. 

  514.       USE module_state_description
    515. 

  515. 

  516.       IMPLICIT NONE
    517. 

  517. 

  518.       INTEGER :: its, ite, jts, jte, kts, kte, &
    519. 

  519.                  i, j, k
    520. 

  520. 

  521.       !  Input data.
    522. 

  522.       TYPE(domain) , POINTER          :: grid
    523. 

  523. 

  524.       ! local
    525. 

  525. !      real p1000mb,eps,svp1,svp2,svp3,svpt0
    526. 

  526.       real eps
    527. 

  527. !      parameter (p1000mb = 1.e+05, eps=0.622,svp1=0.6112,svp3=29.65,svpt0=273.15)
    528. 

  528.       parameter (eps=0.622)
    529. 

  529.       REAL es,qs,pol,tx,temp,pres,rslf             
    530. 

  530.       CHARACTER*80 mess
    531. 

  531. 

  532.       IF ( grid%dfi_opt .EQ. DFI_NODFI ) RETURN
    533. 

  533. 

  534. 

  535.       ! Set dynamic variables
    536. 

  536.       ! divide by total DFI coefficient
    537. 

  537. 

  538. #if (EM_CORE == 1)
    539. 

  539.       grid%mu_2(:,:)    = grid%dfi_mu(:,:)      / grid%hcoeff_tot
    540. 

  540.       grid%u_2(:,:,:)   = grid%dfi_u(:,:,:)     / grid%hcoeff_tot
    541. 

  541.       grid%v_2(:,:,:)   = grid%dfi_v(:,:,:)     / grid%hcoeff_tot
    542. 

  542.       grid%w_2(:,:,:)   = grid%dfi_w(:,:,:)     / grid%hcoeff_tot
    543. 

  543.       grid%ww(:,:,:)    = grid%dfi_ww(:,:,:)    / grid%hcoeff_tot
    544. 

  544.       grid%t_2(:,:,:)   = grid%dfi_t(:,:,:)     / grid%hcoeff_tot
    545. 

  545.       grid%phb(:,:,:)   = grid%dfi_phb(:,:,:)   / grid%hcoeff_tot
    546. 

  546.       grid%ph0(:,:,:)   = grid%dfi_ph0(:,:,:)   / grid%hcoeff_tot
    547. 

  547.       grid%php(:,:,:)   = grid%dfi_php(:,:,:)   / grid%hcoeff_tot
    548. 

  548.       grid%p(:,:,:)     = grid%dfi_p(:,:,:)     / grid%hcoeff_tot
    549. 

  549.       grid%ph_2(:,:,:)  = grid%dfi_ph(:,:,:)    / grid%hcoeff_tot
    550. 

  550.       grid%tke_2(:,:,:) = grid%dfi_tke(:,:,:)   / grid%hcoeff_tot
    551. 

  551.       grid%al(:,:,:)    = grid%dfi_al(:,:,:)    / grid%hcoeff_tot
    552. 

  552.       grid%alt(:,:,:)   = grid%dfi_alt(:,:,:)   / grid%hcoeff_tot
    553. 

  553.       grid%pb(:,:,:)    = grid%dfi_pb(:,:,:)    / grid%hcoeff_tot
    554. 

  554. #if ( WRF_DFI_RADAR == 1 )
    555. 

  555.     IF ( grid%dfi_radar .EQ. 0 ) then   ! tgs no radar assimilation
    556. 

  556.       grid%moist(:,:,:,:)  = max(0.,grid%dfi_moist(:,:,:,:)  / grid%hcoeff_tot)
    557. 

  557.     ELSE
    558. 

  558.       grid%moist(:,:,:,P_QV)  = max(0.,grid%dfi_moist(:,:,:,P_QV)  / grid%hcoeff_tot)
    559. 

  559.     ENDIF
    560. 

  560. #else
    561. 

  561.       grid%moist(:,:,:,:)  = max(0.,grid%dfi_moist(:,:,:,:)  / grid%hcoeff_tot)
    562. 

  562. #endif
    563. 

  563.       grid%scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot)
    564. 

  564. 

  565.       ! restore initial fields
    566. 

  566.       grid%SNOW  (:,:) = grid%dfi_SNOW  (:,:)
    567. 

  567.       grid%SNOWH (:,:) = grid%dfi_SNOWH (:,:)
    568. 

  568.       grid%SNOWC (:,:) = grid%dfi_SNOWC (:,:)
    569. 

  569.       grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
    570. 

  570.       grid%TSK   (:,:) = grid%dfi_TSK   (:,:)
    571. 

  571. 

  572.       grid%TSLB  (:,:,:)       = grid%dfi_TSLB   (:,:,:)
    573. 

  573.       grid%SMOIS (:,:,:)       = grid%dfi_SMOIS  (:,:,:)
    574. 

  574.       IF ( grid%sf_surface_physics .EQ. 3 ) then
    575. 

  575.          grid%QVG   (:,:) = grid%dfi_QVG   (:,:)
    576. 

  576.          grid%TSNAV (:,:) = grid%dfi_TSNAV (:,:)  
    577. 

  577.          grid%SOILT1(:,:) = grid%dfi_SOILT1(:,:)   
    578. 

  578.          grid%SMFR3D(:,:,:)       = grid%dfi_SMFR3D (:,:,:)
    579. 

  579.          grid%KEEPFR3DFLAG(:,:,:) = grid%dfi_KEEPFR3DFLAG(:,:,:)
    580. 

  580.       ENDIF
    581. 

  581. 

  582.       ! restore analized hydrometeor fileds 
    583. 

  583. #if ( WRF_DFI_RADAR == 1 )
    584. 

  584.     IF ( grid%dfi_radar .EQ. 1 ) then
    585. 

  585. !      grid%moist(:,:,:,:)      = grid%dfi_moist(:,:,:,:)    !tgs
    586. 

  586.       grid%moist(:,:,:,P_QC)      = grid%dfi_moist(:,:,:,P_QC)    !tgs
    587. 

  587.       grid%moist(:,:,:,P_QR)      = grid%dfi_moist(:,:,:,P_QR)    !tgs
    588. 

  588.       grid%moist(:,:,:,P_QI)      = grid%dfi_moist(:,:,:,P_QI)    !tgs
    589. 

  589.       grid%moist(:,:,:,P_QS)      = grid%dfi_moist(:,:,:,P_QS)    !tgs
    590. 

  590.       grid%moist(:,:,:,P_QG)      = grid%dfi_moist(:,:,:,P_QG)    !tgs
    591. 

  591. 

  592.        if(grid%dfi_stage .EQ. DFI_FWD) then
    593. 

  593. !tgs change QV to restore initial RH field after the diabatic DFI
    594. 

  594.             its = grid%sp31 ; ite = grid%ep31 ;  
    595. 

  595.             kts = grid%sp32 ; kte = grid%ep32 ; 
    596. 

  596.             jts = grid%sp33 ; jte = grid%ep33 ;
    597. 

  597.    DO j=jts,jte
    598. 

  598.       DO i=its,ite
    599. 

  599.          do k = kts , kte
    600. 

  600.           temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
    601. 

  601.                      ** (r_d / Cp)
    602. 

  602.           pres = grid%p(i,k,j)+grid%pb(i,k,j)
    603. 

  603. !tgs rslf - function to compute qs from Thompson microphysics
    604. 

  604.           qs = rslf(pres, temp)
    605. 

  605. 

  606. !      if(i.eq. 178 .and. j.eq. 148 .and. k.eq.11) then
    607. 

  607. !       print *,'temp,pres,qs-thomp',temp,pres,qs
    608. 

  608. !      endif
    609. 

  609. 

  610.       IF(grid%moist(i,k,j,P_QC) .GT. 1.e-6 .or.                    &
    611. 

  611.          grid%moist(i,k,j,P_QI) .GT. 1.e-6) THEN
    612. 

  612.          grid%moist (i,k,j,P_QV) = MAX(0.,grid%dfi_rh(i,k,j)*QS)
    613. 

  613.       ENDIF
    614. 

  614. 

  615. !      if(i.eq. 178 .and. j.eq. 148 .and. k.eq.11) then
    616. 

  616. !       print *,'temp,pres,qs,grid%moist (i,k,j,P_QV)',temp,pres,qs,  &
    617. 

  617. !                grid%moist(i,k,j,P_QV)
    618. 

  618. !      endif
    619. 

  619.          enddo
    620. 

  620.       ENDDO
    621. 

  621.    ENDDO
    622. 

  622.        endif
    623. 

  623. 

  624.     ENDIF
    625. 

  625. #endif
    626. 

  626. #endif
    627. 

  627. 

  628. #if (NMM_CORE == 1)
    629. 

  629.         write(mess,*) ' divide by grid%hcoeff_tot: ', grid%hcoeff_tot
    630. 

  630.         call wrf_message(mess)
    631. 

  631. 	if (grid%hcoeff_tot .lt. 0.0001) then
    632. 

  632. 	call wrf_error_fatal("bad grid%hcoeff_tot")
    633. 

  633. 	endif
    634. 

  634.        grid%pd(:,:)       = grid%dfi_pd(:,:)  / grid%hcoeff_tot
    635. 

  635.        grid%pint(:,:,:)      = grid%dfi_pint(:,:,:) / grid%hcoeff_tot
    636. 

  636. !       grid%dwdt(:,:,:)      = grid%dfi_dwdt(:,:,:) / grid%hcoeff_tot
    637. 

  637.        grid%t(:,:,:)      = grid%dfi_t(:,:,:) / grid%hcoeff_tot
    638. 

  638.        grid%q(:,:,:)      = grid%dfi_q(:,:,:) / grid%hcoeff_tot
    639. 

  639.        grid%q2(:,:,:)     = grid%dfi_q2(:,:,:) / grid%hcoeff_tot
    640. 

  640.        grid%cwm(:,:,:)    = grid%dfi_cwm(:,:,:) / grid%hcoeff_tot
    641. 

  641.        grid%u(:,:,:)      = grid%dfi_u(:,:,:) / grid%hcoeff_tot
    642. 

  642.        grid%v(:,:,:)      = grid%dfi_v(:,:,:) / grid%hcoeff_tot
    643. 

  643.       grid%moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  / grid%hcoeff_tot
    644. 

  644.       grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot
    645. 

  645. 

  646.       ! restore initial fields
    647. 

  647.       grid%SNOW(:,:)   = grid%dfi_SNOW(:,:)
    648. 

  648.       grid%SNOWH(:,:)  = grid%dfi_SNOWH(:,:)
    649. 

  649. !      grid%SNOWC(:,:)  = grid%dfi_SNOWC(:,:)
    650. 

  650.       grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
    651. 

  651.       grid%NMM_TSK(:,:)    = grid%dfi_NMM_TSK(:,:)
    652. 

  652.       ! save soil fields
    653. 

  653.       grid%STC(:,:,:)         = grid%dfi_STC(:,:,:)
    654. 

  654.       grid%SMC(:,:,:)        = grid%dfi_SMC(:,:,:)
    655. 

  655.       grid%SH2O(:,:,:)       = grid%dfi_SH2O(:,:,:)
    656. 

  656. #endif
    657. 

  657. 

  658. 

  659.    END SUBROUTINE dfi_array_reset
    660. 

  660. 

  661.    SUBROUTINE optfil_driver( grid )
    662. 

  662. 

  663.       USE module_domain, ONLY : domain
    664. 

  664.       USE module_utility
    665. 

  665. !      USE module_wrf_error
    666. 

  666. !      USE module_timing
    667. 

  667. !      USE module_date_time
    668. 

  668. !      USE module_configure
    669. 

  669.       USE module_state_description
    670. 

  670.       USE module_model_constants
    671. 

  671. 

  672.       IMPLICIT NONE
    673. 

  673. 

  674.       TYPE (domain) , POINTER                 :: grid
    675. 

  675. 

  676.       ! Local variables
    677. 

  677.       integer :: nstep2, nstepmax, rundfi, i, rc,hr,min,sec
    678. 

  678.       integer :: yr,jday
    679. 

  679.       real :: timestep, tauc
    680. 

  680.       TYPE(WRFU_TimeInterval) :: run_interval
    681. 

  681.       CHARACTER*80 mess
    682. 

  682. 

  683.       timestep=abs(grid%dt)
    684. 

  684.       run_interval = grid%stop_subtime - grid%start_subtime
    685. 

  685.       CALL WRFU_TimeGet(grid%start_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
    686. 

  686.       CALL WRFU_TimeGet(grid%stop_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
    687. 

  687. 

  688.       CALL WRFU_TimeIntervalGet( run_interval, S=rundfi, rc=rc )
    689. 

  689.       rundfi = abs(rundfi)
    690. 

  690. 

  691.       nstep2= ceiling((1.0 + real(rundfi)/timestep) / 2.0)
    692. 

  692. 

  693.       ! nstep2 is equal to a half of timesteps per initialization period,
    694. 

  694.       ! should not exceed nstepmax
    695. 

  695. 

  696.       tauc = real(grid%dfi_cutoff_seconds)
    697. 

  697. 

  698.       ! Get DFI coefficient
    699. 

  699.       grid%hcoeff(:) = 0.0
    700. 

  700.       IF ( grid%dfi_nfilter < 0 .OR. grid%dfi_nfilter > 8 ) THEN
    701. 

  701. write(mess,*) 'Invalid filter specified in namelist.'
    702. 

  702. call wrf_message(mess)
    703. 

  703.       ELSE
    704. 

  704.          call dfcoef(nstep2-1, grid%dt, tauc, grid%dfi_nfilter, grid%hcoeff)
    705. 

  705.       END IF
    706. 

  706. 

  707.       IF ( MOD(int(1.0 + real(rundfi)/timestep),2) /= 0 ) THEN
    708. 

  708.          DO i=1,nstep2-1 
    709. 

  709.             grid%hcoeff(2*nstep2-i) = grid%hcoeff(i)
    710. 

  710.          END DO
    711. 

  711.       ELSE
    712. 

  712.          DO i=1,nstep2 
    713. 

  713.             grid%hcoeff(2*nstep2-i+1) = grid%hcoeff(i)
    714. 

  714.          END DO
    715. 

  715.       END IF
    716. 

  716. 

  717.    END SUBROUTINE optfil_driver
    718. 

  718. 

  719. 

  720.    SUBROUTINE dfi_clear_accumulation( grid )
    721. 

  721. 

  722.       USE module_domain, ONLY : domain
    723. 

  723. !      USE module_configure
    724. 

  724. !      USE module_driver_constants
    725. 

  725. !      USE module_machine
    726. 

  726. !      USE module_dm
    727. 

  727. !      USE module_model_constants
    728. 

  728.       USE module_state_description
    729. 

  729. 

  730.       IMPLICIT NONE
    731. 

  731. 

  732.       !  Input data.
    733. 

  733.       TYPE(domain) , POINTER          :: grid
    734. 

  734. 

  735. #if (EM_CORE == 1)
    736. 

  736.       grid%dfi_mu(:,:)       = 0.
    737. 

  737.       grid%dfi_u(:,:,:)      = 0.
    738. 

  738.       grid%dfi_v(:,:,:)      = 0.
    739. 

  739.       grid%dfi_w(:,:,:)      = 0.
    740. 

  740.       grid%dfi_ww(:,:,:)     = 0.
    741. 

  741.       grid%dfi_t(:,:,:)      = 0.
    742. 

  742.       grid%dfi_phb(:,:,:)    = 0.
    743. 

  743.       grid%dfi_ph0(:,:,:)    = 0.
    744. 

  744.       grid%dfi_php(:,:,:)    = 0.
    745. 

  745.       grid%dfi_p(:,:,:)      = 0.
    746. 

  746.       grid%dfi_ph(:,:,:)     = 0.
    747. 

  747.       grid%dfi_tke(:,:,:)    = 0.
    748. 

  748.       grid%dfi_al(:,:,:)     = 0.
    749. 

  749.       grid%dfi_alt(:,:,:)    = 0.
    750. 

  750.       grid%dfi_pb(:,:,:)     = 0.
    751. 

  751. #if ( WRF_DFI_RADAR == 1 )
    752. 

  752.     IF ( grid%dfi_radar .EQ. 0 ) then  
    753. 

  753.       grid%dfi_moist(:,:,:,:)  = 0.
    754. 

  754.     ELSE
    755. 

  755.       grid%dfi_moist(:,:,:,P_QV) = 0.
    756. 

  756.     ENDIF
    757. 

  757. #else
    758. 

  758.       grid%dfi_moist(:,:,:,:)  = 0.
    759. 

  759. #endif
    760. 

  760.       grid%dfi_scalar(:,:,:,:) = 0.
    761. 

  761. #endif
    762. 

  762. 

  763. #if (NMM_CORE == 1)
    764. 

  764.        grid%dfi_pd(:,:)    = 0.
    765. 

  765.        grid%dfi_pint(:,:,:)   = 0.
    766. 

  766.        grid%dfi_dwdt(:,:,:)   = 0.
    767. 

  767.        grid%dfi_t(:,:,:)   = 0.
    768. 

  768.        grid%dfi_q(:,:,:)   = 0.
    769. 

  769.        grid%dfi_q2(:,:,:)  = 0.
    770. 

  770.        grid%dfi_cwm(:,:,:) = 0.
    771. 

  771.        grid%dfi_u(:,:,:)   = 0.
    772. 

  772.        grid%dfi_v(:,:,:)   = 0.
    773. 

  773.       grid%dfi_moist(:,:,:,:)  = 0.
    774. 

  774.       grid%dfi_scalar(:,:,:,:) = 0.
    775. 

  775. #endif
    776. 

  776. 

  777.       grid%hcoeff_tot  = 0.0
    778. 

  778. 

  779.    END SUBROUTINE dfi_clear_accumulation
    780. 

  780. 

  781. 

  782.    SUBROUTINE dfi_save_arrays( grid )
    783. 

  783. 

  784.       USE module_domain, ONLY : domain
    785. 

  785. !      USE module_configure
    786. 

  786. !      USE module_driver_constants
    787. 

  787. !      USE module_machine
    788. 

  788. !      USE module_dm
    789. 

  789.       USE module_model_constants
    790. 

  790.       USE module_state_description
    791. 

  791. 

  792.       IMPLICIT NONE
    793. 

  793. 

  794.       INTEGER :: its, ite, jts, jte, kts, kte, &
    795. 

  795.                  i, j, k
    796. 

  796. 

  797.       !  Input data.
    798. 

  798.       TYPE(domain) , POINTER          :: grid
    799. 

  799.       ! local
    800. 

  800. 

  801.       REAL es,qs,pol,tx,temp,pres,rslf
    802. 

  802. 

  803. #if (EM_CORE == 1)
    804. 

  804.       ! save surface 2-D fields
    805. 

  805.       grid%dfi_SNOW(:,:)   = grid%SNOW(:,:) 
    806. 

  806.       grid%dfi_SNOWH(:,:)  = grid%SNOWH(:,:) 
    807. 

  807.       grid%dfi_SNOWC(:,:)  = grid%SNOWC(:,:) 
    808. 

  808.       grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:) 
    809. 

  809.       grid%dfi_TSK(:,:)    = grid%TSK(:,:) 
    810. 

  810. 

  811.       ! save soil fields
    812. 

  812.       grid%dfi_TSLB(:,:,:)         = grid%TSLB(:,:,:) 
    813. 

  813.       grid%dfi_SMOIS(:,:,:)        = grid%SMOIS(:,:,:) 
    814. 

  814.       ! RUC LSM only, need conditional
    815. 

  815.       IF ( grid%sf_surface_physics .EQ. 3 ) then
    816. 

  816.          grid%dfi_QVG(:,:)            = grid%QVG(:,:) 
    817. 

  817.          grid%dfi_SOILT1(:,:)         = grid%SOILT1(:,:) 
    818. 

  818.          grid%dfi_TSNAV(:,:)          = grid%TSNAV(:,:) 
    819. 

  819.          grid%dfi_SMFR3D(:,:,:)       = grid%SMFR3D(:,:,:) 
    820. 

  820.          grid%dfi_KEEPFR3DFLAG(:,:,:) = grid%KEEPFR3DFLAG(:,:,:) 
    821. 

  821.       ENDIF
    822. 

  822. #endif
    823. 

  823. 

  824. #if (NMM_CORE == 1)
    825. 

  825.       ! save surface 2-D fields
    826. 

  826.       grid%dfi_SNOW(:,:)   = grid%SNOW(:,:)
    827. 

  827.       grid%dfi_SNOWH(:,:)  = grid%SNOWH(:,:)
    828. 

  828. !      grid%dfi_SNOWC(:,:)  = grid%SNOWC(:,:)
    829. 

  829.       grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
    830. 

  830.       grid%dfi_NMM_TSK(:,:)    = grid%NMM_TSK(:,:)
    831. 

  831.       ! save soil fields
    832. 

  832.       grid%dfi_STC(:,:,:)         = grid%STC(:,:,:)
    833. 

  833.       grid%dfi_SMC(:,:,:)        = grid%SMC(:,:,:)
    834. 

  834.       grid%dfi_SH2O(:,:,:)       = grid%SH2O(:,:,:)
    835. 

  835. #endif
    836. 

  836. 

  837. 

  838.       ! save hydrometeor fields 
    839. 

  839. #if (EM_CORE == 1)
    840. 

  840. #if ( WRF_DFI_RADAR == 1 )
    841. 

  841.     IF ( grid%dfi_radar .EQ. 1 ) then    !tgs
    842. 

  842. !      grid%dfi_moist(:,:,:,:)  =   grid%moist(:,:,:,:)
    843. 

  843.       grid%dfi_moist(:,:,:,P_QC)  =   grid%moist(:,:,:,P_QC)
    844. 

  844.       grid%dfi_moist(:,:,:,P_QR)  =   grid%moist(:,:,:,P_QR)
    845. 

  845.       grid%dfi_moist(:,:,:,P_QI)  =   grid%moist(:,:,:,P_QI)
    846. 

  846.       grid%dfi_moist(:,:,:,P_QS)  =   grid%moist(:,:,:,P_QS)
    847. 

  847.       grid%dfi_moist(:,:,:,P_QG)  =   grid%moist(:,:,:,P_QG)
    848. 

  848. 

  849.        if(grid%dfi_stage .EQ. DFI_BCK) then
    850. 

  850. ! compute initial RH field to be reintroduced after the diabatic DFI
    851. 

  851.             its = grid%sp31 ; ite = grid%ep31 ;  
    852. 

  852.             kts = grid%sp32 ; kte = grid%ep32 ; 
    853. 

  853.             jts = grid%sp33 ; jte = grid%ep33 ;
    854. 

  854.    DO j=jts,jte
    855. 

  855.       DO i=its,ite
    856. 

  856.          do k = kts , kte
    857. 

  857.           temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
    858. 

  858.                      ** (r_d / Cp)
    859. 

  859.           pres = grid%p(i,k,j)+grid%pb(i,k,j)
    860. 

  860. !tgs rslf - function to compute qs from Thompson microphysics
    861. 

  861.           qs = rslf(pres, temp)
    862. 

  862.           grid%dfi_rh (i,k,j) = MIN(1.,MAX(0.,grid%moist(i,k,j,P_QV)/qs))
    863. 

  863. 

  864. !tgs saturation check for points with water or ice clouds
    865. 

  865.       IF ((grid%moist (i,k,j,P_QC) .GT. 1.e-6  .or.        &
    866. 

  866.           grid%moist (i,k,j,P_QI) .GT. 1.e-6)  .and.       &
    867. 

  867.           grid%dfi_rh (i,k,j) .lt. 1.) THEN 
    868. 

  868.                grid%dfi_rh (i,k,j)=1.
    869. 

  869.       ENDIF
    870. 

  870.          
    871. 

  871.         end do
    872. 

  872.       END DO
    873. 

  873.     ENDDO
    874. 

  874.        endif
    875. 

  875. 

  876.     ENDIF
    877. 

  877. #endif
    878. 

  878. #endif
    879. 

  879. 

  880.    END SUBROUTINE dfi_save_arrays
    881. 

  881. 

  882. 

  883.    SUBROUTINE dfcoef (NSTEPS,DT,TAUC,NFILT,H)
    884. 

  884. !
    885. 

  885. !     calculate filter weights with selected window.
    886. 

  886. !
    887. 

  887. !       peter lynch and xiang-yu huang
    888. 

  888. !
    889. 

  889. !       ref: see hamming, r.w., 1989: digital filters,
    890. 

  890. !                prentice-hall international. 3rd edition.
    891. 

  891. !
    892. 

  892. !       input:      nsteps  -  number of timesteps
    893. 

  893. !                              forward or backward.
    894. 

  894. !                   dt      -  time step in seconds.
    895. 

  895. !                   tauc    -  cut-off period in seconds.
    896. 

  896. !                   nfilt   -  indicator for selected filter.
    897. 

  897. !
    898. 

  898. !       output:     h       -  array(0:nsteps) with the
    899. 

  899. !                              required filter weights
    900. 

  900. !
    901. 

  901. !------------------------------------------------------------
    902. 

  902.  
    903. 

  903.       implicit none
    904. 

  904. 

  905.       integer, intent(in)   :: nsteps, nfilt
    906. 

  906.       real   , intent(in)   :: dt, tauc
    907. 

  907.       real, intent(out)  :: h(1:nsteps+1)
    908. 

  908.       
    909. 

  909.       ! Local data
    910. 

  910. 

  911.       integer               :: n
    912. 

  912.       real                  :: pi, omegac, x, window, deltat
    913. 

  913.       real                  :: hh(0:nsteps)
    914. 

  914. 

  915.       pi=4*ATAN(1.)
    916. 

  916.       deltat=ABS(dt)
    917. 

  917. 

  918.       ! windows are defined by a call and held in hh.
    919. 

  919. 

  920.       if ( nfilt .eq. -1) call debug    (nsteps,h)
    921. 

  921. 

  922.       IF ( NFILT .EQ. 0 ) CALL UNIFORM  (NSTEPS,HH)
    923. 

  923.       IF ( NFILT .EQ. 1 ) CALL LANCZOS  (NSTEPS,HH)
    924. 

  924.       IF ( NFILT .EQ. 2 ) CALL HAMMING  (NSTEPS,HH)
    925. 

  925.       IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
    926. 

  926.       IF ( NFILT .EQ. 4 ) CALL KAISER   (NSTEPS,HH)
    927. 

  927.       IF ( NFILT .EQ. 5 ) CALL POTTER2  (NSTEPS,HH)
    928. 

  928.       IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
    929. 

  929. 

  930.       IF ( NFILT .LE. 6 ) THEN     ! sinc-windowed filters
    931. 

  931. 

  932.          ! calculate the cutoff frequency
    933. 

  933.          OMEGAC = 2.*PI/TAUC
    934. 

  934. 

  935.          DO N=0,NSTEPS
    936. 

  936.             WINDOW = HH(N)
    937. 

  937.             IF ( N .EQ. 0 ) THEN
    938. 

  938.               X = (OMEGAC*DELTAT/PI)
    939. 

  939.             ELSE
    940. 

  940.               X = SIN(N*OMEGAC*DELTAT)/(N*PI)
    941. 

  941.             END IF
    942. 

  942.             HH(N) = X*WINDOW
    943. 

  943.          END DO
    944. 

  944. 

  945.          ! normalize the sums to be unity
    946. 

  946.          CALL NORMLZ(HH,NSTEPS)
    947. 

  947. 

  948.          DO N=0,NSTEPS
    949. 

  949.             H(N+1) = HH(NSTEPS-N)
    950. 

  950.          END DO
    951. 

  951. 

  952.       ELSE IF ( NFILT .EQ. 7 ) THEN     ! dolph filter
    953. 

  953. 

  954.          CALL DOLPH(DT,TAUC,NSTEPS,H)
    955. 

  955. 

  956.       ELSE IF ( NFILT .EQ. 8 ) THEN     ! 2nd order, 2nd type quick start filter (RHO)
    957. 

  957. 

  958.          CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS)
    959. 

  959. 

  960.       END IF
    961. 

  961. 

  962.       RETURN
    963. 

  963. 

  964.    END SUBROUTINE dfcoef
    965. 

  965. 

  966. 

  967.    SUBROUTINE NORMLZ(HH,NMAX)
    968. 

  968. 

  969.    ! normalize the sum of hh to be unity
    970. 

  970. 

  971.       implicit none
    972. 

  972.   
    973. 

  973.       integer, intent(in)                       :: nmax
    974. 

  974.       real   , dimension(0:nmax), intent(out)   :: hh
    975. 

  975. 

  976.       ! local data
    977. 

  977.       real     :: sumhh
    978. 

  978.       integer  :: n
    979. 

  979. 

  980.       SUMHH = HH(0)
    981. 

  981.       DO N=1,NMAX
    982. 

  982.         SUMHH = SUMHH + 2*HH(N)
    983. 

  983.       ENDDO
    984. 

  984.       DO N=0,NMAX
    985. 

  985.         HH(N)  = HH(N)/SUMHH
    986. 

  986.       ENDDO
    987. 

  987. 

  988.       RETURN
    989. 

  989. 

  990.    END subroutine normlz
    991. 

  991. 

  992. 

  993.    subroutine debug(nsteps, ww)
    994. 

  994. 

  995.       implicit none
    996. 

  996. 

  997.       integer, intent(in)                        :: nsteps
    998. 

  998.       real   , dimension(0:nsteps), intent(out)  :: ww
    999. 

  999.       integer :: n
    1000. 

  1000. 

  1001.       do n=0,nsteps
    1002. 

  1002.          ww(n)=0
    1003. 

  1003.       end do
    1004. 

  1004. 

  1005.       ww(int(nsteps/2))=1
    1006. 

  1006. 

  1007.       return
    1008. 

  1008. 

  1009.    end subroutine debug
    1010. 

  1010. 

  1011. 

  1012.    SUBROUTINE UNIFORM(NSTEPS,WW)
    1013. 

  1013. 

  1014.    !  define uniform or rectangular window function.
    1015. 

  1015. 

  1016.       implicit none
    1017. 

  1017. 

  1018.       integer, intent(in)                        :: nsteps
    1019. 

  1019.       real   , dimension(0:nsteps), intent(out)  :: ww
    1020. 

  1020.        
    1021. 

  1021.       integer          :: n
    1022. 

  1022. 

  1023.       DO N=0,NSTEPS
    1024. 

  1024.         WW(N) = 1.
    1025. 

  1025.       ENDDO
    1026. 

  1026. 

  1027.       RETURN
    1028. 

  1028. 

  1029.    END subroutine uniform
    1030. 

  1030. 

  1031. 

  1032.    SUBROUTINE LANCZOS(NSTEPS,WW)
    1033. 

  1033. 

  1034.    ! define (genaralised) lanczos window function.
    1035. 

  1035. 

  1036.       implicit none 
    1037. 

  1037. 

  1038.       integer,  parameter                      :: nmax = 1000
    1039. 

  1039.       integer,  intent(in)                     :: nsteps
    1040. 

  1040.       real   ,  dimension(0:nmax), intent(out) :: ww
    1041. 

  1041.       integer  ::  n
    1042. 

  1042.       real     :: power, pi, w
    1043. 

  1043. 

  1044.       ! (for the usual lanczos window, power = 1 )
    1045. 

  1045.       POWER = 1
    1046. 

  1046. 

  1047.       PI=4*ATAN(1.)
    1048. 

  1048.       DO N=0,NSTEPS
    1049. 

  1049.          IF ( N .EQ. 0 ) THEN
    1050. 

  1050.             W = 1.0
    1051. 

  1051.          ELSE
    1052. 

  1052.             W = SIN(N*PI/(NSTEPS+1)) / ( N*PI/(NSTEPS+1))
    1053. 

  1053.          ENDIF
    1054. 

  1054.          WW(N) = W**POWER
    1055. 

  1055.       ENDDO
    1056. 

  1056. 

  1057.       RETURN
    1058. 

  1058. 

  1059.    END SUBROUTINE lanczos
    1060. 

  1060. 

  1061. 

  1062.    SUBROUTINE HAMMING(NSTEPS,WW)
    1063. 

  1063. 

  1064.    ! define (genaralised) hamming window function.
    1065. 

  1065. 

  1066.       implicit none
    1067. 

  1067. 

  1068.       integer, intent(in)           :: nsteps
    1069. 

  1069.       real, dimension(0:nsteps)    :: ww
    1070. 

  1070.       integer   ::   n
    1071. 

  1071.       real      :: alpha, pi, w
    1072. 

  1072. 

  1073.       ! (for the usual hamming window, alpha=0.54,
    1074. 

  1074.       !      for the hann window, alpha=0.50).
    1075. 

  1075.       ALPHA=0.54
    1076. 

  1076. 

  1077.       PI=4*ATAN(1.)
    1078. 

  1078.       DO N=0,NSTEPS
    1079. 

  1079.          IF ( N .EQ. 0 ) THEN
    1080. 

  1080.             W = 1.0
    1081. 

  1081.          ELSE
    1082. 

  1082.             W = ALPHA + (1-ALPHA)*COS(N*PI/(NSTEPS))
    1083. 

  1083.          ENDIF
    1084. 

  1084.          WW(N) = W
    1085. 

  1085.       ENDDO
    1086. 

  1086. 

  1087.       RETURN
    1088. 

  1088. 

  1089.    END SUBROUTINE hamming
    1090. 

  1090. 

  1091. 

  1092.    SUBROUTINE BLACKMAN(NSTEPS,WW)
    1093. 

  1093. 

  1094.    ! define blackman window function.
    1095. 

  1095. 

  1096.       implicit none
    1097. 

  1097. 

  1098.       integer, intent(in)           :: nsteps
    1099. 

  1099.       real, dimension(0:nsteps)    :: ww
    1100. 

  1100.       integer  :: n
    1101. 

  1101. 

  1102.       real :: pi, w
    1103. 

  1103. 

  1104.       PI=4*ATAN(1.)
    1105. 

  1105.       DO N=0,NSTEPS
    1106. 

  1106.          IF ( N .EQ. 0 ) THEN
    1107. 

  1107.             W = 1.0
    1108. 

  1108.          ELSE
    1109. 

  1109.             W = 0.42 + 0.50*COS(  N*PI/(NSTEPS))   &
    1110. 

  1110.                      + 0.08*COS(2*N*PI/(NSTEPS))
    1111. 

  1111.          ENDIF
    1112. 

  1112.          WW(N) = W
    1113. 

  1113.       ENDDO
    1114. 

  1114. 

  1115.       RETURN
    1116. 

  1116. 

  1117.    END SUBROUTINE blackman
    1118. 

  1118. 

  1119. 

  1120.    SUBROUTINE KAISER(NSTEPS,WW)
    1121. 

  1121. 

  1122.    ! define kaiser window function.
    1123. 

  1123. 

  1124.       implicit none
    1125. 

  1125. 

  1126.       real, external :: bessi0
    1127. 

  1127. 

  1128.       integer, intent(in)           :: nsteps
    1129. 

  1129.       real, dimension(0:nsteps)    :: ww
    1130. 

  1130.       integer   :: n
    1131. 

  1131.       real      :: alpha, xi0a, xn, as
    1132. 

  1132. 

  1133.       ALPHA = 1
    1134. 

  1134. 

  1135.       XI0A =  BESSI0(ALPHA)
    1136. 

  1136.       DO N=0,NSTEPS
    1137. 

  1137.          XN = N
    1138. 

  1138.          AS = ALPHA*SQRT(1.-(XN/NSTEPS)**2)
    1139. 

  1139.          WW(N) = BESSI0(AS) / XI0A
    1140. 

  1140.       ENDDO
    1141. 

  1141. 

  1142.       RETURN
    1143. 

  1143. 

  1144.    END SUBROUTINE kaiser
    1145. 

  1145. 

  1146. 

  1147.    REAL FUNCTION BESSI0(X)
    1148. 

  1148. 

  1149.    ! from numerical recipes (press, et al.)
    1150. 

  1150. 

  1151.       implicit none
    1152. 

  1152. 

  1153.       real(8)   :: Y
    1154. 

  1154.       real(8)   :: P1 = 1.0d0
    1155. 

  1155.       real(8)   :: P2 = 3.5156230D0
    1156. 

  1156.       real(8)   :: P3 = 3.0899424D0
    1157. 

  1157.       real(8)   :: P4 = 1.2067492D0
    1158. 

  1158.       real(8)   :: P5 = 0.2659732D0
    1159. 

  1159.       real(8)   :: P6 = 0.360768D-1
    1160. 

  1160.       real(8)   :: P7 = 0.45813D-2
    1161. 

  1161. 

  1162.       real*8    :: Q1 = 0.39894228D0
    1163. 

  1163.       real*8    :: Q2 = 0.1328592D-1
    1164. 

  1164.       real*8    :: Q3 = 0.225319D-2
    1165. 

  1165.       real*8    :: Q4 = -0.157565D-2
    1166. 

  1166.       real*8    :: Q5 = 0.916281D-2
    1167. 

  1167.       real*8    :: Q6 = -0.2057706D-1
    1168. 

  1168.       real*8    :: Q7 = 0.2635537D-1
    1169. 

  1169.       real*8    :: Q8 = -0.1647633D-1
    1170. 

  1170.       real*8    :: Q9 = 0.392377D-2
    1171. 

  1171. 

  1172.       real            :: x, ax
    1173. 

  1173. 

  1174. 

  1175.       IF (ABS(X).LT.3.75) THEN
    1176. 

  1176.         Y=(X/3.75)**2
    1177. 

  1177.         BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
    1178. 

  1178.       ELSE
    1179. 

  1179.         AX=ABS(X)
    1180. 

  1180.         Y=3.75/AX
    1181. 

  1181.         BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4    &
    1182. 

  1182.            +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
    1183. 

  1183.       ENDIF
    1184. 

  1184.       RETURN
    1185. 

  1185. 

  1186.    END FUNCTION bessi0
    1187. 

  1187. 

  1188. 

  1189.    SUBROUTINE POTTER2(NSTEPS,WW)
    1190. 

  1190. 

  1191.       ! define potter window function.
    1192. 

  1192.       ! modified to fall off over twice the range.
    1193. 

  1193. 

  1194.       implicit none
    1195. 

  1195. 

  1196.       integer, intent(in)                       :: nsteps
    1197. 

  1197.       real, dimension(0:nsteps),intent(out)     ::  ww
    1198. 

  1198.       integer  :: n
    1199. 

  1199.       real     :: ck, sum, arg
    1200. 

  1200. 

  1201.       ! local data
    1202. 

  1202.       real, dimension(0:3)   :: d 
    1203. 

  1203.       real                   :: pi
    1204. 

  1204.       integer                :: ip
    1205. 

  1205. 

  1206.       d(0) = 0.35577019
    1207. 

  1207.       d(1) = 0.2436983
    1208. 

  1208.       d(2) = 0.07211497
    1209. 

  1209.       d(3) = 0.00630165
    1210. 

  1210. 

  1211.       PI=4*ATAN(1.)
    1212. 

  1212. 

  1213.       CK = 1.0
    1214. 

  1214.       DO N=0,NSTEPS
    1215. 

  1215.          IF (N.EQ.NSTEPS) CK = 0.5
    1216. 

  1216.          ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
    1217. 

  1217. !min---        modification in next statement
    1218. 

  1218.          ARG = ARG/2.
    1219. 

  1219. !min---        end of modification
    1220. 

  1220.          SUM = D(0)
    1221. 

  1221.          DO IP=1,3
    1222. 

  1222.             SUM = SUM + 2.*D(IP)*COS(ARG*FLOAT(IP))
    1223. 

  1223.          END DO
    1224. 

  1224.          WW(N) = CK*SUM
    1225. 

  1225.       END DO
    1226. 

  1226. 

  1227.       RETURN
    1228. 

  1228. 

  1229.    END SUBROUTINE potter2
    1230. 

  1230. 

  1231. 

  1232.    SUBROUTINE dolphwin(m, window)
    1233. 

  1233. 

  1234. !     calculation of dolph-chebyshev window or, for short,
    1235. 

  1235. !     dolph window, using the expression in the reference:
    1236. 

  1236. !
    1237. 

  1237. !     antoniou, andreas, 1993: digital filters: analysis,
    1238. 

  1238. !     design and applications. mcgraw-hill, inc., 689pp.
    1239. 

  1239. !
    1240. 

  1240. !     the dolph window is optimal in the following sense:
    1241. 

  1241. !     for a given main-lobe width, the stop-band attenuation
    1242. 

  1242. !     is minimal; for a given stop-band level, the main-lobe
    1243. 

  1243. !     width is minimal.
    1244. 

  1244. !
    1245. 

  1245. !     it is possible to specify either the ripple-ratio r
    1246. 

  1246. !     or the stop-band edge thetas.
    1247. 

  1247. 

  1248.       IMPLICIT NONE
    1249. 

  1249. 

  1250.       ! Arguments
    1251. 

  1251.       INTEGER, INTENT(IN)                  ::  m
    1252. 

  1252.       REAL, DIMENSION(0:M), INTENT(OUT)    ::  window
    1253. 

  1253. 

  1254.       ! local data
    1255. 

  1255.       REAL, DIMENSION(0:2*M)               :: t
    1256. 

  1256.       REAL, DIMENSION(0:M)                 :: w, time
    1257. 

  1257.       REAL    :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
    1258. 

  1258.       INTEGER :: n, nm1, nt, i
    1259. 

  1259. 

  1260.       PI = 4*ATAN(1.D0)
    1261. 

  1261.       THETAS = 2*PI/M
    1262. 

  1262. 

  1263.       N = 2*M+1
    1264. 

  1264.       NM1 = N-1
    1265. 

  1265.       X0 = 1/COS(THETAS/2)
    1266. 

  1266. 

  1267.       TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
    1268. 

  1268.       TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
    1269. 

  1269.       RR = 0.5*(TERM1+TERM2)
    1270. 

  1270.       R = 1/RR
    1271. 

  1271.       DB = 20*LOG10(R)
    1272. 

  1272.       WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
    1273. 

  1273.       WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
    1274. 

  1274.       WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
    1275. 

  1275. 

  1276.       DO NT=0,M
    1277. 

  1277.         SUM = RR
    1278. 

  1278.         DO I=1,M
    1279. 

  1279.           ARG = X0*COS(I*PI/N)
    1280. 

  1280.           CALL CHEBY(T,NM1,ARG)
    1281. 

  1281.           TERM1 = T(NM1)
    1282. 

  1282.           TERM2 = COS(2*NT*PI*I/N)
    1283. 

  1283.           SUM = SUM + 2*TERM1*TERM2
    1284. 

  1284.         ENDDO
    1285. 

  1285.         W(NT) = SUM/N
    1286. 

  1286.         TIME(NT) = NT
    1287. 

  1287.       ENDDO
    1288. 

  1288. 

  1289. !     fill up the array for return
    1290. 

  1290.       DO NT=0,M
    1291. 

  1291.         WINDOW(NT) = W(NT)
    1292. 

  1292.       ENDDO
    1293. 

  1293. 

  1294.       RETURN
    1295. 

  1295. 

  1296.    END SUBROUTINE dolphwin
    1297. 

  1297. 

  1298. 

  1299.    SUBROUTINE dolph(deltat, taus, m, window)
    1300. 

  1300. 

  1301. !     calculation of dolph-chebyshev window or, for short,
    1302. 

  1302. !     dolph window, using the expression in the reference:
    1303. 

  1303. !
    1304. 

  1304. !     antoniou, andreas, 1993: digital filters: analysis,
    1305. 

  1305. !     design and applications. mcgraw-hill, inc., 689pp.
    1306. 

  1306. !
    1307. 

  1307. !     the dolph window is optimal in the following sense:
    1308. 

  1308. !     for a given main-lobe width, the stop-band attenuation
    1309. 

  1309. !     is minimal; for a given stop-band level, the main-lobe
    1310. 

  1310. !     width is minimal.
    1311. 

  1311. 

  1312.       IMPLICIT NONE
    1313. 

  1313. 

  1314.       ! Arguments
    1315. 

  1315.       INTEGER, INTENT(IN)                  ::  m
    1316. 

  1316.       REAL, DIMENSION(0:M), INTENT(OUT)    ::  window
    1317. 

  1317.       REAL, INTENT(IN)                     :: deltat, taus
    1318. 

  1318. 

  1319.       ! local data
    1320. 

  1320.       integer, PARAMETER        :: NMAX = 5000
    1321. 

  1321.       REAL, dimension(0:NMAX)   :: t, w, time
    1322. 

  1322.       real, dimension(0:2*nmax) :: w2
    1323. 

  1323.       INTEGER                   :: NPRPE=0        ! no of pe
    1324. 

  1324.       CHARACTER*80              :: MES
    1325. 

  1325. 

  1326.       real    :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
    1327. 

  1327.       integer :: n, nm1, i, nt
    1328. 

  1328.       
    1329. 

  1329.       PI = 4*ATAN(1.D0)
    1330. 

  1330. 

  1331.       WRITE (mes,'(A,F8.2,A,F10.2)') 'In dolph, deltat = ',deltat,' taus = ',taus
    1332. 

  1332.       CALL wrf_message(TRIM(mes))
    1333. 

  1333. 

  1334.       N = 2*M+1
    1335. 

  1335.       NM1 = N-1
    1336. 

  1336. 

  1337.       THETAS = 2*PI*ABS(DELTAT/TAUS)
    1338. 

  1338.       X0 = 1/COS(THETAS/2)
    1339. 

  1339.       TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
    1340. 

  1340.       TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
    1341. 

  1341.       RR = 0.5*(TERM1+TERM2)
    1342. 

  1342.       R = 1/RR
    1343. 

  1343.       DB = 20*LOG10(R)
    1344. 

  1344. 

  1345.       WRITE (mes,'(A,2I8)') 'In dolph: M,N = ', M,N
    1346. 

  1346.       CALL wrf_message(TRIM(mes))
    1347. 

  1347.       WRITE (mes,'(A,F10.3)') 'In dolph: THETAS (STOP-BAND EDGE) = ', thetas
    1348. 

  1348.       CALL wrf_message(TRIM(mes))
    1349. 

  1349.       WRITE (mes,'(A,2F10.3)') 'In dolph: R,DB = ', R,DB
    1350. 

  1350.       CALL wrf_message(TRIM(mes))
    1351. 

  1351. 

  1352.       DO NT=0,M
    1353. 

  1353.          SUM = 1
    1354. 

  1354.          DO I=1,M
    1355. 

  1355.             ARG = X0*COS(I*PI/N)
    1356. 

  1356.             CALL CHEBY(T,NM1,ARG)
    1357. 

  1357.             TERM1 = T(NM1)
    1358. 

  1358.             TERM2 = COS(2*NT*PI*I/N)
    1359. 

  1359.             SUM = SUM + R*2*TERM1*TERM2
    1360. 

  1360.          ENDDO
    1361. 

  1361.          W(NT) = SUM/N
    1362. 

  1362.          TIME(NT) = NT
    1363. 

  1363.          WRITE (mes,'(A,F10.6,2x,E17.7)') 'In dolph: TIME, W = ', TIME(NT), W(NT)
    1364. 

  1364.          CALL wrf_message(TRIM(mes))
    1365. 

  1365.       ENDDO
    1366. 

  1366. !     fill in the negative-time values by symmetry.
    1367. 

  1367.       DO NT=0,M
    1368. 

  1368.          W2(M+NT) = W(NT)
    1369. 

  1369.          W2(M-NT) = W(NT)
    1370. 

  1370.       ENDDO
    1371. 

  1371. 

  1372. !     fill up the array for return
    1373. 

  1373.       SUMW = 0.
    1374. 

  1374.       DO NT=0,2*M
    1375. 

  1375.          SUMW = SUMW + W2(NT)
    1376. 

  1376.       ENDDO
    1377. 

  1377.       WRITE (mes,'(A,F10.4)') 'In dolph: SUM OF WEIGHTS W2 = ', sumw
    1378. 

  1378.       CALL wrf_message(TRIM(mes))
    1379. 

  1379. 

  1380.       DO NT=0,M
    1381. 

  1381.          WINDOW(NT) = W2(NT)
    1382. 

  1382.       ENDDO
    1383. 

  1383. 

  1384.       RETURN
    1385. 

  1385. 

  1386.    END SUBROUTINE dolph
    1387. 

  1387. 

  1388. 

  1389.    SUBROUTINE cheby(t, n, x)
    1390. 

  1390. 

  1391. !     calculate all chebyshev polynomials up to order n
    1392. 

  1392. !     for the argument value x.
    1393. 

  1393. 

  1394. !     reference: numerical recipes, page 184, recurrence
    1395. 

  1395. !         t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) ,  n>=2.
    1396. 

  1396. 

  1397.       IMPLICIT NONE
    1398. 

  1398. 

  1399.       ! Arguments
    1400. 

  1400.       INTEGER, INTENT(IN)  :: n
    1401. 

  1401.       REAL, INTENT(IN)     :: x
    1402. 

  1402.       REAL, DIMENSION(0:N) :: t
    1403. 

  1403.  
    1404. 

  1404.       integer  :: nn
    1405. 

  1405. 

  1406.       T(0) = 1
    1407. 

  1407.       T(1) = X
    1408. 

  1408.       IF(N.LT.2) RETURN
    1409. 

  1409.       DO NN=2,N
    1410. 

  1410.          T(NN) = 2*X*T(NN-1) - T(NN-2)
    1411. 

  1411.       ENDDO
    1412. 

  1412. 

  1413.       RETURN
    1414. 

  1414. 

  1415.    END SUBROUTINE cheby
    1416. 

  1416. 

  1417. 

  1418.    SUBROUTINE rhofil(dt, tauc, norder, nstep, ictype, frow, nosize)
    1419. 

  1419. 

  1420. !       RHO = recurssive high order.
    1421. 

  1421. !
    1422. 

  1422. !       This routine calculates and returns the 
    1423. 

  1423. !       Last Row, FROW, of the FILTER matrix.
    1424. 

  1424. !
    1425. 

  1425. !       Input Parameters:
    1426. 

  1426. !              DT  :  Time Step in seconds
    1427. 

  1427. !            TAUC  :  Cut-off period (hours)
    1428. 

  1428. !          NORDER  :  Order of QS Filter
    1429. 

  1429. !           NSTEP  :  Number of step/Size of row.
    1430. 

  1430. !          ICTYPE  :  Initial Conditions
    1431. 

  1431. !          NOSIZE  :  Max. side of FROW.
    1432. 

  1432. !
    1433. 

  1433. !       Working Fields:
    1434. 

  1434. !           ACOEF  :  X-coefficients of filter
    1435. 

  1435. !           BCOEF  :  Y-coefficients of filter
    1436. 

  1436. !          FILTER  :  Filter Matrix.
    1437. 

  1437. !
    1438. 

  1438. !       Output Parameters:
    1439. 

  1439. !            FROW  : Last Row of Filter Matrix.
    1440. 

  1440. !
    1441. 

  1441. !       Note: Two types of initial conditions are permitted.
    1442. 

  1442. !          ICTYPE = 1 : Order increasing each row to NORDER.
    1443. 

  1443. !          ICTYPE = 2 : Order fixed at NORDER throughout.
    1444. 

  1444. !
    1445. 

  1445. !       DOUBLE PRECISION USED THROUGHOUT.
    1446. 

  1446. 

  1447.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    1448. 

  1448. 

  1449.       DOUBLE PRECISION MUC
    1450. 

  1450. 

  1451. !       N.B. Single Precision for List Parameters.
    1452. 

  1452.       REAL, intent(in)    ::  DT,TAUC
    1453. 

  1453. 

  1454. !       Space for the last row of FILTER.
    1455. 

  1455.       integer, intent(in) ::  norder, nstep, ictype, nosize
    1456. 

  1456.       REAL   , dimension(0:nosize), intent(out)::  FROW
    1457. 

  1457. 

  1458. !       Arrays for rho filter.
    1459. 

  1459.       integer, PARAMETER  :: NOMAX=100
    1460. 

  1460.       real   , dimension(0:NOMAX) :: acoef, bcoef
    1461. 

  1461.       real   , dimension(0:NOMAX,0:NOMAX) :: filter
    1462. 

  1462. !       Working space.
    1463. 

  1463.       real   , dimension(0:NOMAX) :: alpha, beta
    1464. 

  1464. 

  1465.       real   :: DTT
    1466. 

  1466. 

  1467.       DTT = ABS(DT)
    1468. 

  1468.       PI = 2*DASIN(1.D0)
    1469. 

  1469.       IOTA = CMPLX(0.,1.)
    1470. 

  1470. 

  1471. !       Filtering Parameters (derived).
    1472. 

  1472.       THETAC = 2*PI*DTT/(TAUC)
    1473. 

  1473.       MUC = tan(THETAC/2) 
    1474. 

  1474.       FC = THETAC/(2*PI)
    1475. 

  1475. 

  1476. !     Clear the arrays.
    1477. 

  1477.     DO NC=0,NOMAX
    1478. 

  1478.        ACOEF(NC) = 0.
    1479. 

  1479.        BCOEF(NC) = 0.
    1480. 

  1480.        ALPHA(NC) = 0.
    1481. 

  1481.        BETA (NC) = 0.
    1482. 

  1482.        FROW (NC) = 0.
    1483. 

  1483.        DO NR=0,NOMAX
    1484. 

  1484.           FILTER(NR,NC) = 0.
    1485. 

  1485.        ENDDO
    1486. 

  1486.     ENDDO
    1487. 

  1487. 

  1488. !     Fill up the Filter Matrix.
    1489. 

  1489.     FILTER(0,0) = 1.
    1490. 

  1490. 

  1491. !     Get the coefficients of the Filter.
    1492. 

  1492.     IF ( ICTYPE.eq.2 ) THEN
    1493. 

  1493.        CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
    1494. 

  1494.     ENDIF
    1495. 

  1495. 

  1496.     DO 100 NROW=1,NSTEP
    1497. 

  1497. 

  1498.        IF ( ICTYPE.eq.1 ) THEN
    1499. 

  1499.           NORD = MIN(NROW,NORDER)
    1500. 

  1500.           IF ( NORD.le.NORDER) THEN
    1501. 

  1501.              CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
    1502. 

  1502.           ENDIF
    1503. 

  1503.        ENDIF
    1504. 

  1504. 

  1505.        DO K=0,NROW
    1506. 

  1506.           ALPHA(K) = ACOEF(NROW-K)
    1507. 

  1507.           IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
    1508. 

  1508.        ENDDO
    1509. 

  1509. 

  1510. !        Correction for terms of negative index.
    1511. 

  1511.        IF ( ICTYPE.eq.2 ) THEN
    1512. 

  1512.           IF ( NROW.lt.NORDER ) THEN
    1513. 

  1513.              CN = 0.
    1514. 

  1514.              DO NN=NROW+1,NORDER
    1515. 

  1515.                 CN = CN + (ACOEF(NN)+BCOEF(NN))
    1516. 

  1516.              ENDDO
    1517. 

  1517.              ALPHA(0) = ALPHA(0) + CN
    1518. 

  1518.           ENDIF
    1519. 

  1519.        ENDIF
    1520. 

  1520. 

  1521. !       Check sum of ALPHAs and BETAs = 1
    1522. 

  1522.       SUMAB = 0.
    1523. 

  1523.       DO NN=0,NROW
    1524. 

  1524.         SUMAB = SUMAB + ALPHA(NN)
    1525. 

  1525.           IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
    1526. 

  1526.       ENDDO
    1527. 

  1527. 

  1528.       DO KK=0,NROW-1
    1529. 

  1529.          SUMBF = 0.
    1530. 

  1530.          DO LL=0,NROW-1
    1531. 

  1531.             SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
    1532. 

  1532.          ENDDO
    1533. 

  1533.          FILTER(NROW,KK) = ALPHA(KK)+SUMBF
    1534. 

  1534.       ENDDO
    1535. 

  1535.       FILTER(NROW,NROW) = ALPHA(NROW)
    1536. 

  1536. 

  1537. !       Check sum of row elements = 1
    1538. 

  1538.       SUMROW = 0.
    1539. 

  1539.       DO NN=0,NROW
    1540. 

  1540.         SUMROW = SUMROW + FILTER(NROW,NN)
    1541. 

  1541.       ENDDO
    1542. 

  1542. 

  1543. 100 CONTINUE
    1544. 

  1544. 

  1545.       DO NC=0,NSTEP
    1546. 

  1546.         FROW(NC) = FILTER(NSTEP,NC)
    1547. 

  1547.       ENDDO
    1548. 

  1548. 

  1549.       RETURN
    1550. 

  1550. 

  1551.    END SUBROUTINE rhofil
    1552. 

  1552. 

  1553. 

  1554.    SUBROUTINE rhocof(nord, nomax, muc, ca, cb)
    1555. 

  1555. 

  1556. !     Get the coefficients of the RHO Filter 
    1557. 

  1557. 

  1558. !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    1559. 

  1559.       IMPLICIT NONE 
    1560. 

  1560. 

  1561.       ! Arguments
    1562. 

  1562.       integer, intent(in)      :: nord, nomax
    1563. 

  1563.       real, dimension(0:nomax) :: ca, cb
    1564. 

  1564. 

  1565.       ! Functions
    1566. 

  1566.       double precision, external :: cnr
    1567. 

  1567. 

  1568.       ! Local variables
    1569. 

  1569.       INTEGER                  :: nn
    1570. 

  1570.       COMPLEX                  :: IOTA
    1571. 

  1571.       DOUBLE PRECISION         :: MUC, ZN
    1572. 

  1572.       DOUBLE PRECISION         :: pi, root2, rn, sigma, gain, sumcof
    1573. 

  1573. 

  1574.       PI = 2*ASIN(1.)
    1575. 

  1575.       ROOT2 = SQRT(2.)
    1576. 

  1576.       IOTA = (0.,1.)
    1577. 

  1577. 

  1578.       RN = 1./FLOAT(NORD)
    1579. 

  1579.       SIGMA = 1./( SQRT(2.**RN-1.) )
    1580. 

  1580. 

  1581.       GAIN  = (MUC*SIGMA/(1+MUC*SIGMA))**NORD
    1582. 

  1582.       ZN = (1-MUC*SIGMA)/(1+MUC*SIGMA)
    1583. 

  1583. 

  1584.       DO NN=0,NORD
    1585. 

  1585.         CA(NN) = CNR(NORD,NN)*GAIN
    1586. 

  1586.         IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)*(-ZN)**NN
    1587. 

  1587.       ENDDO
    1588. 

  1588. 

  1589. !     Check sum of coefficients = 1
    1590. 

  1590.       SUMCOF = 0.
    1591. 

  1591.       DO NN=0,NORD
    1592. 

  1592.         SUMCOF = SUMCOF + CA(NN)
    1593. 

  1593.         IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
    1594. 

  1594.       ENDDO
    1595. 

  1595. 

  1596.       RETURN
    1597. 

  1597. 

  1598.    END SUBROUTINE RHOCOF
    1599. 

  1599. 

  1600. 

  1601.    DOUBLE PRECISION FUNCTION cnr(n,r)
    1602. 

  1602. 

  1603.    ! Binomial Coefficient (n,r).
    1604. 

  1604. 

  1605. !      IMPLICIT DOUBLE PRECISION(C,X)
    1606. 

  1606.       IMPLICIT NONE
    1607. 

  1607. 

  1608.       ! Arguments
    1609. 

  1609.       INTEGER , intent(in)    :: n, R
    1610. 

  1610. 

  1611.       ! Local variables
    1612. 

  1612.       INTEGER          :: k
    1613. 

  1613.       DOUBLE PRECISION :: coeff, xn, xr, xk
    1614. 

  1614. 

  1615.       IF ( R.eq.0 ) THEN
    1616. 

  1616.          CNR = 1.0
    1617. 

  1617.          RETURN
    1618. 

  1618.       ENDIF
    1619. 

  1619.       Coeff = 1.0
    1620. 

  1620.       XN = DFLOAT(N)
    1621. 

  1621.       XR = DFLOAT(R)
    1622. 

  1622.       DO K=1,R
    1623. 

  1623.         XK = DFLOAT(K)
    1624. 

  1624.         COEFF = COEFF * ( (XN-XR+XK)/XK )
    1625. 

  1625.       ENDDO
    1626. 

  1626.       CNR = COEFF
    1627. 

  1627. 

  1628.       RETURN
    1629. 

  1629. 

  1630.    END FUNCTION cnr
    1631. 

  1631. 

  1632. 

  1633.     SUBROUTINE optfil (grid,NH,DELTAT,NHMAX)
    1634. 

  1634. !----------------------------------------------------------------------
    1635. 

  1635. 

  1636. !    SUBROUTINE optfil (NH,DELTAT,TAUP,TAUS,LPRINT,         &
    1637. 

  1637. !                         H,NHMAX)
    1638. 

  1638. !
    1639. 

  1639. ! - Huang and Lynch optimal filter
    1640. 

  1640. !   Monthly Weather Review, Feb 1993
    1641. 

  1641. !----------------------------------------------------------
    1642. 

  1642. !       Input Parameters in List:
    1643. 

  1643. !             NH     :  Half-length of the Filter
    1644. 

  1644. !             DELTAT :  Time-step (in seconds).      
    1645. 

  1645. !             TAUP   :  Period of pass-band edge (hours).
    1646. 

  1646. !             TAUS   :  Period of stop-band edge (hours).
    1647. 

  1647. !             LPRINT :  Logical switch for messages.
    1648. 

  1648. !             NHMAX  :  Maximum permitted Half-length.
    1649. 

  1649. !
    1650. 

  1650. !       Output Parameters in List:
    1651. 

  1651. !             H      :  Impulse Response.
    1652. 

  1652. !             DP     :  Deviation in pass-band (db)
    1653. 

  1653. !             DS     :  Deviation in stop-band (db)
    1654. 

  1654. !----------------------------------------------------------
    1655. 

  1655. !
    1656. 

  1656.      USE module_domain, ONLY : domain
    1657. 

  1657.     
    1658. 

  1658.      TYPE(domain) , POINTER :: grid
    1659. 

  1659. 

  1660.      REAL,DIMENSION( 20) :: EDGE
    1661. 

  1661.      REAL,DIMENSION( 10) :: FX, WTX, DEVIAT
    1662. 

  1662.      REAL,DIMENSION(2*NHMAX+1) :: H
    1663. 

  1663.      logical LPRINT
    1664. 

  1664.      REAL, INTENT (IN) :: DELTAT
    1665. 

  1665.      INTEGER, INTENT (IN) :: NH, NHMAX
    1666. 

  1666. !
    1667. 

  1667.          TAUP = 3.
    1668. 

  1668.          TAUS = 1.5
    1669. 

  1669.          LPRINT = .true.
    1670. 

  1670. !initialize H array
    1671. 

  1671. 

  1672.          NL=2*NHMAX+1
    1673. 

  1673.         do 101 n=1,NL
    1674. 

  1674.           H(n)=0.
    1675. 

  1675.  101    continue
    1676. 

  1676. 

  1677.         NFILT = 2*NH+1
    1678. 

  1678.         print *,' start optfil, NFILT=', nfilt
    1679. 

  1679. 

  1680. !
    1681. 

  1681. ! 930325 PL & XYH : the upper limit is changed from 64 to 128.
    1682. 

  1682.         IF(NFILT.LE.0 .OR. NFILT.GT.128 ) THEN
    1683. 

  1683. 	   WRITE(6,*) 'NH=',NH
    1684. 

  1684.        CALL wrf_error_fatal (' Sorry, error 1 in call to OPTFIL ')
    1685. 

  1685. 	ENDIF
    1686. 

  1686. !
    1687. 

  1687. !       The following four should always be the same.
    1688. 

  1688.         JTYPE = 1
    1689. 

  1689.         NBANDS = 2
    1690. 

  1690. !CC     JPRINT = 0
    1691. 

  1691.         LGRID = 16
    1692. 

  1692. !
    1693. 

  1693. !       calculate transition frequencies.
    1694. 

  1694. 	  DT = ABS(DELTAT)
    1695. 

  1695.           FS = DT/(TAUS*3600.)
    1696. 

  1696.           FP = DT/(TAUP*3600.)
    1697. 

  1697.           IF(FS.GT.0.5) then
    1698. 

  1698. !            print *,' FS too large in OPTFIL '
    1699. 

  1699.        CALL wrf_error_fatal (' FS too large in OPTFIL ')
    1700. 

  1700. !            return
    1701. 

  1701.           end if
    1702. 

  1702.           IF(FP.LT.0.0) then
    1703. 

  1703. !            print *, ' FP too small in OPTFIL '
    1704. 

  1704.        CALL wrf_error_fatal (' FP too small in OPTFIL ')
    1705. 

  1705. !            return
    1706. 

  1706.           end if
    1707. 

  1707. !
    1708. 

  1708. !       Relative Weights in pass- and stop-bands.
    1709. 

  1709.           WTP = 1.0
    1710. 

  1710.           WTS = 1.0
    1711. 

  1711. !
    1712. 

  1712. !CC     NOTE: (FP,FC,FS) is an arithmetic progression, so
    1713. 

  1713. !CC           (1/FS,1/FC,1/FP) is a harmonic one.
    1714. 

  1714. !CC           TAUP = 1/( (1/TAUC)-(1/DTAU) )
    1715. 

  1715. !CC           TAUS = 1/( (1/TAUC)+(1/DTAU) )
    1716. 

  1716. !CC           TAUC   :  Cut-off Period (hours).
    1717. 

  1717. !CC           DTAU   :  Transition half-width (hours).
    1718. 

  1718. !CC           FC = 1/TAUC  ;  DF = 1/DTAU
    1719. 

  1719. !CC           FP = FC - DF :  FS = FC + DF
    1720. 

  1720. !
    1721. 

  1721.           IF ( LPRINT ) THEN
    1722. 

  1722.                TAUC = 2./((1/TAUS)+(1/TAUP))
    1723. 

  1723.                DTAU = 2./((1/TAUS)-(1/TAUP))
    1724. 

  1724.                FC = DT/(TAUC*3600.)
    1725. 

  1725.                DF = DT/(DTAU*3600.)
    1726. 

  1726.                WRITE(6,*) ' DT ' , dt
    1727. 

  1727.                WRITE(6,*) ' TAUS, TAUP ' , TAUS,TAUP
    1728. 

  1728.                WRITE(6,*) ' TAUC, DTAU ' , TAUC,DTAU
    1729. 

  1729.                WRITE(6,*) '   FP,   FS ' ,   FP,  FS   
    1730. 

  1730.                WRITE(6,*) '   FC,   DF ' ,   FC,  DF
    1731. 

  1731.                WRITE(6,*) '  WTS,  WTP ' ,  WTS, WTP
    1732. 

  1732.           ENDIF
    1733. 

  1733. !
    1734. 

  1734. !       Fill the control vectors for MCCPAR
    1735. 

  1735.         EDGE(1) = 0.0
    1736. 

  1736.         EDGE(2) = FP
    1737. 

  1737.         EDGE(3) = FS
    1738. 

  1738.         EDGE(4) = 0.5
    1739. 

  1739.         FX(1) = 1.0
    1740. 

  1740.         FX(2) = 0.0
    1741. 

  1741.         WTX(1) = WTP
    1742. 

  1742.         WTX(2) = WTS
    1743. 

  1743. 

  1744.         CALL MCCPAR(NFILT,JTYPE,NBANDS,LPRINT,LGRID,       &
    1745. 

  1745.                    EDGE,FX,WTX,DEVIAT, h )
    1746. 

  1746. !
    1747. 

  1747. !       Save the deviations in the pass- and stop-bands.
    1748. 

  1748.         DP = DEVIAT(1)
    1749. 

  1749.         DS = DEVIAT(2)
    1750. 

  1750. !
    1751. 

  1751. !     Fill out the array H (only first half filled in MCCPAR).
    1752. 

  1752.       IF(MOD(NFILT,2).EQ.0) THEN
    1753. 

  1753.          NHALF = ( NFILT )/2
    1754. 

  1754.       ELSE
    1755. 

  1755.          NHALF = (NFILT+1)/2
    1756. 

  1756.       ENDIF
    1757. 

  1757.       DO 100 nn=1,NHALF
    1758. 

  1758.          H(NFILT+1-nn) = h(nn)
    1759. 

  1759.   100 CONTINUE
    1760. 

  1760. 

  1761. !       normalize the sums to be unity
    1762. 

  1762.         sumh = 0 
    1763. 

  1763.         do 150 n=1,NFILT
    1764. 

  1764.           sumh = sumh + H(n)
    1765. 

  1765.   150   continue
    1766. 

  1766.   print *,'SUMH =', sumh
    1767. 

  1767. 

  1768.         do 200 n=1,NFILT 
    1769. 

  1769.           H(n)  = H(n)/sumh
    1770. 

  1770.   200   continue
    1771. 

  1771.         do 201 n=1,NFILT
    1772. 

  1772.         grid%hcoeff(n)=H(n)
    1773. 

  1773.   201   continue
    1774. 

  1774. !       print *,'HCOEFF(n) ', grid%hcoeff
    1775. 

  1775. ! 
    1776. 

  1776.    END SUBROUTINE optfil
    1777. 

  1777. 

  1778. 

  1779.       SUBROUTINE MCCPAR (NFILT,JTYPE,NBANDS,LPRINT,LGRID,    &
    1780. 

  1780.                    EDGE,FX,WTX,DEVIAT,h )
    1781. 

  1781. 

  1782. !      PROGRAM FOR THE DESIGN OF LINEAR PHASE FINITE IMPULSE
    1783. 

  1783. !      REPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM
    1784. 

  1784. !
    1785. 

  1785. !************************************************************
    1786. 

  1786. !*  Reference: McClellan, J.H., T.W. Parks and L.R.Rabiner, *
    1787. 

  1787. !*             1973:  A computer program for designing      *
    1788. 

  1788. !*             optimum FIR linear phase digital filters.    *
    1789. 

  1789. !*             IEEE Trans. on Audio and Electroacoustics,   *
    1790. 

  1790. !*             Vol AU-21, No. 6, 506-526.                   *
    1791. 

  1791. !************************************************************
    1792. 

  1792. !
    1793. 

  1793. !      THREE TYPES OF FILTERS ARE INCLUDED -- BANDPASS FILTERS
    1794. 

  1794. !      DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS
    1795. 

  1795. !
    1796. 

  1796. !---------------------------------------------------------------
    1797. 

  1797. !
    1798. 

  1798. !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
    1799. 

  1799.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
    1800. 

  1800.       DIMENSION H(66)
    1801. 

  1801.       DIMENSION DES(1045),GRID(1045),WT(1045)
    1802. 

  1802.       DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
    1803. 

  1803.       DOUBLE PRECISION PI2,PI
    1804. 

  1804.       DOUBLE PRECISION AD,DEV,X,Y
    1805. 

  1805.       LOGICAL LPRINT
    1806. 

  1806.       
    1807. 

  1807.       PI  = 3.141592653589793
    1808. 

  1808.       PI2 = 6.283185307179586
    1809. 

  1809. 

  1810. !      ......
    1811. 

  1811. 

  1812.       NFMAX = 128
    1813. 

  1813. 100      CONTINUE
    1814. 

  1814. 

  1815. !      PROGRAM INPUT SECTION
    1816. 

  1816. 

  1817. !CC     READ(5,*) NFILT,JTYPE,NBANDS,JPRINT,LGRID
    1818. 

  1818. 

  1819.       IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) THEN
    1820. 

  1820.          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
    1821. 

  1821.       END IF
    1822. 

  1822.       IF(NBANDS.LE.0) NBANDS = 1
    1823. 

  1823. 

  1824. !      ....
    1825. 

  1825. 

  1826.       IF(LGRID.LE.0) LGRID = 16
    1827. 

  1827.       JB = 2*NBANDS
    1828. 

  1828. !cc       READ(5,*) (EDGE(J),J=1,JB)
    1829. 

  1829. !cc       READ(5,*) (FX(J),J=1,NBANDS)
    1830. 

  1830. !cc       READ(5,*) (WTX(J),J=1,NBANDS)
    1831. 

  1831.       IF(JTYPE.EQ.0) THEN
    1832. 

  1832.          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
    1833. 

  1833.       END IF
    1834. 

  1834.       NEG = 1
    1835. 

  1835.       IF(JTYPE.EQ.1) NEG = 0
    1836. 

  1836.       NODD = NFILT/2
    1837. 

  1837.       NODD = NFILT-2*NODD
    1838. 

  1838.       NFCNS = NFILT/2
    1839. 

  1839.       IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS = NFCNS+1
    1840. 

  1840. 

  1841. !      ...
    1842. 

  1842. 

  1843.       GRID(1) = EDGE(1)
    1844. 

  1844.       DELF = LGRID*NFCNS
    1845. 

  1845.       DELF = 0.5/DELF
    1846. 

  1846.       IF(NEG.EQ.0) GOTO 135
    1847. 

  1847.       IF(EDGE(1).LT.DELF) GRID(1) = DELF
    1848. 

  1848. 135      CONTINUE
    1849. 

  1849.       J = 1
    1850. 

  1850.       L = 1
    1851. 

  1851.       LBAND = 1
    1852. 

  1852. 140      FUP = EDGE(L+1)
    1853. 

  1853. 145      TEMP = GRID(J)
    1854. 

  1854. 

  1855. !      ....
    1856. 

  1856. 

  1857.       DES(J) = EFF(TEMP,FX,WTX,LBAND,JTYPE)
    1858. 

  1858.       WT(J) = WATE(TEMP,FX,WTX,LBAND,JTYPE)
    1859. 

  1859.       J = J+1
    1860. 

  1860.       GRID(J) = TEMP+DELF
    1861. 

  1861.       IF(GRID(J).GT.FUP) GOTO 150
    1862. 

  1862.       GOTO 145
    1863. 

  1863. 150      GRID(J-1) = FUP
    1864. 

  1864.       DES(J-1) = EFF(FUP,FX,WTX,LBAND,JTYPE)
    1865. 

  1865.       WT(J-1) = WATE(FUP,FX,WTX,LBAND,JTYPE)
    1866. 

  1866.       LBAND = LBAND+1
    1867. 

  1867.       L = L+2
    1868. 

  1868.       IF(LBAND.GT.NBANDS) GOTO 160
    1869. 

  1869.       GRID(J) = EDGE(L)
    1870. 

  1870.       GOTO 140
    1871. 

  1871. 160      NGRID = J-1
    1872. 

  1872.       IF(NEG.NE.NODD) GOTO 165
    1873. 

  1873.       IF(GRID(NGRID).GT.(0.5-DELF)) NGRID = NGRID-1
    1874. 

  1874. 165      CONTINUE
    1875. 

  1875. 

  1876. !      ......
    1877. 

  1877. 

  1878.       IF(NEG) 170,170,180
    1879. 

  1879. 170      IF(NODD.EQ.1) GOTO 200
    1880. 

  1880.       DO 175 J=1,NGRID
    1881. 

  1881.             CHANGE = DCOS(PI*GRID(J))
    1882. 

  1882.             DES(J) = DES(J)/CHANGE
    1883. 

  1883.             WT(J) = WT(J)*CHANGE
    1884. 

  1884. 175      CONTINUE
    1885. 

  1885.       GOTO 200
    1886. 

  1886. 180      IF(NODD.EQ.1) GOTO 190
    1887. 

  1887.       DO 185 J = 1,NGRID
    1888. 

  1888.             CHANGE = DSIN(PI*GRID(J))
    1889. 

  1889.             DES(J) = DES(J)/CHANGE
    1890. 

  1890.             WT(J)  = WT(J)*CHANGE
    1891. 

  1891. 185      CONTINUE
    1892. 

  1892.       GOTO 200
    1893. 

  1893. 190      DO 195 J =1,NGRID
    1894. 

  1894.             CHANGE = DSIN(PI2*GRID(J))
    1895. 

  1895.             DES(J) = DES(J)/CHANGE
    1896. 

  1896.             WT(J)  = WT(J)*CHANGE
    1897. 

  1897. 195      CONTINUE
    1898. 

  1898. 

  1899. !      ......
    1900. 

  1900. 

  1901. 200      TEMP = FLOAT(NGRID-1)/FLOAT(NFCNS)
    1902. 

  1902.       DO 210 J = 1,NFCNS
    1903. 

  1903.             IEXT(J) = (J-1)*TEMP+1
    1904. 

  1904. 210      CONTINUE
    1905. 

  1905.       IEXT(NFCNS+1) = NGRID
    1906. 

  1906.       NM1 = NFCNS-1
    1907. 

  1907.       NZ  = NFCNS+1
    1908. 

  1908. 

  1909. ! CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM
    1910. 

  1910. 

  1911.       CALL REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
    1912. 

  1912. 

  1913. ! CALCULATE THE IMPULSE RESPONSE
    1914. 

  1914. 

  1915.       IF(NEG) 300,300,320
    1916. 

  1916. 300      IF(NODD.EQ.0) GOTO 310
    1917. 

  1917.       DO 305 J=1,NM1
    1918. 

  1918.             H(J) = 0.5*ALPHA(NZ-J)
    1919. 

  1919. 305      CONTINUE
    1920. 

  1920.       H(NFCNS)=ALPHA(1)
    1921. 

  1921.       GOTO 350
    1922. 

  1922. 310      H(1) = 0.25*ALPHA(NFCNS)
    1923. 

  1923.       DO 315 J = 2,NM1
    1924. 

  1924.             H(J) = 0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
    1925. 

  1925. 315      CONTINUE
    1926. 

  1926.       H(NFCNS) = 0.5*ALPHA(1)+0.25*ALPHA(2)
    1927. 

  1927.       GOTO 350
    1928. 

  1928. 320      IF(NODD.EQ.0) GOTO 330
    1929. 

  1929.       H(1) = 0.25*ALPHA(NFCNS)
    1930. 

  1930.       H(2) = 0.25*ALPHA(NM1)
    1931. 

  1931.       DO 325 J = 3,NM1
    1932. 

  1932.             H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
    1933. 

  1933. 325      CONTINUE
    1934. 

  1934.       H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(3)
    1935. 

  1935.       H(NZ) = 0.0
    1936. 

  1936.       GOTO 350
    1937. 

  1937. 330      H(1) = 0.25*ALPHA(NFCNS)
    1938. 

  1938.       DO 335 J =2,NM1
    1939. 

  1939.             H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
    1940. 

  1940. 335      CONTINUE
    1941. 

  1941.       H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(2)
    1942. 

  1942. 

  1943. !   PROGRAM OUTPUT SECTION
    1944. 

  1944. 

  1945. 350   CONTINUE
    1946. 

  1946. !
    1947. 

  1947.       IF(LPRINT) THEN
    1948. 

  1948. 

  1949.          print *, '****************************************************'
    1950. 

  1950.          print *, 'FINITE IMPULSE RESPONSE (FIR)'
    1951. 

  1951.          print *, 'LINEAR PHASE DIGITAL FILTER DESIGN'
    1952. 

  1952.          print *, 'REMEZ EXCHANGE ALGORITHM'
    1953. 

  1953.       
    1954. 

  1954.       IF(JTYPE.EQ.1) WRITE(6,365)
    1955. 

  1955. 365      FORMAT(25X,'BANDPASS FILTER'/)
    1956. 

  1956. 

  1957.       IF(JTYPE.EQ.2) WRITE(6,370)
    1958. 

  1958. 370      FORMAT(25X,'DIFFERENTIATOR '/)
    1959. 

  1959. 

  1960.       IF(JTYPE.EQ.3) WRITE(6,375)
    1961. 

  1961. 375      FORMAT(25X,'HILBERT TRANSFORMER '/)
    1962. 

  1962. 

  1963.       WRITE(6,378) NFILT
    1964. 

  1964. 378      FORMAT(15X,'FILTER LENGTH =',I3/)
    1965. 

  1965. 

  1966.       WRITE(6,380)
    1967. 

  1967. 380      FORMAT(15X,'***** IMPULSE RESPONSE *****')
    1968. 

  1968. 

  1969.       DO 381 J = 1,NFCNS
    1970. 

  1970.             K = NFILT+1-J
    1971. 

  1971.             IF(NEG.EQ.0) WRITE(6,382) J,H(J),K
    1972. 

  1972.             IF(NEG.EQ.1) WRITE(6,383) J,H(J),K
    1973. 

  1973. 381      CONTINUE
    1974. 

  1974. 382      FORMAT(20X,'H(',I3,') = ',E15.8,' = H(',I4,')')
    1975. 

  1975. 383      FORMAT(20X,'H(',I3,') = ',E15.8,' = -H(',I4,')')
    1976. 

  1976. 

  1977.       IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ
    1978. 

  1978. 384      FORMAT(20X,'H(',I3,') = 0.0')
    1979. 

  1979. 

  1980.       DO 450 K=1,NBANDS,4
    1981. 

  1981.             KUP = K+3
    1982. 

  1982.             IF(KUP.GT.NBANDS) KUP = NBANDS
    1983. 

  1983. 	    print *
    1984. 

  1984.             WRITE(6,385) (J,J=K,KUP)
    1985. 

  1985. 385            FORMAT(24X,4('BAND',I3,8X))
    1986. 

  1986.             WRITE(6,390) (EDGE(2*J-1),J=K,KUP)
    1987. 

  1987. 390            FORMAT(2X,'LOWER BAND EDGE',5F15.8)
    1988. 

  1988.             WRITE(6,395) (EDGE(2*J),J=K,KUP)
    1989. 

  1989. 395            FORMAT(2X,'UPPER BAND EDGE',5F15.8)
    1990. 

  1990.             IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
    1991. 

  1991. 400            FORMAT(2X,'DESIRED VALUE',2X,5F15.8)
    1992. 

  1992.             IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
    1993. 

  1993. 405            FORMAT(2X,'DESIRED SLOPE',2X,5F15.8)
    1994. 

  1994.             WRITE(6,410) (WTX(J),J=K,KUP)
    1995. 

  1995. 410            FORMAT(2X,'WEIGHTING',6X,5F15.8)
    1996. 

  1996.             DO 420 J = K,KUP
    1997. 

  1997.                   DEVIAT(J) = DEV/WTX(J)
    1998. 

  1998. 420            CONTINUE
    1999. 

  1999.             WRITE(6,425) (DEVIAT(J),J=K,KUP)
    2000. 

  2000. 425            FORMAT(2X,'DEVIATION',6X,5F15.8)
    2001. 

  2001.             IF(JTYPE.NE.1) GOTO 450
    2002. 

  2002.             DO 430 J = K,KUP
    2003. 

  2003.                   DEVIAT(J) = 20.0*ALOG10(DEVIAT(J))
    2004. 

  2004. 430            CONTINUE
    2005. 

  2005.             WRITE(6,435) (DEVIAT(J),J=K,KUP)
    2006. 

  2006. 435            FORMAT(2X,'DEVIATION IN DB',5F15.8)
    2007. 

  2007. 450      CONTINUE
    2008. 

  2008.       print *, 'EXTREMAL FREQUENCIES'
    2009. 

  2009.       WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
    2010. 

  2010. 455      FORMAT((2X,5F15.7))
    2011. 

  2011.       WRITE(6,460)
    2012. 

  2012. 460      FORMAT(1X,70(1H*))
    2013. 

  2013. !
    2014. 

  2014.       ENDIF
    2015. 

  2015. !
    2016. 

  2016. !CC     IF(NFILT.NE.0) GOTO 100  ! removal of re-run loop.
    2017. 

  2017. !
    2018. 

  2018.    END SUBROUTINE mccpar
    2019. 

  2019. 

  2020. 

  2021.       FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE)
    2022. 

  2022.          DIMENSION FX(5),WTX(5)
    2023. 

  2023.          IF(JTYPE.EQ.2) GOTO 1
    2024. 

  2024.          EFF = FX(LBAND)
    2025. 

  2025.          RETURN
    2026. 

  2026. 1        EFF = FX(LBAND)*TEMP
    2027. 

  2027.       END FUNCTION eff
    2028. 

  2028. 

  2029. 

  2030.       FUNCTION WATE(TEMP,FX,WTX,LBAND,JTYPE)
    2031. 

  2031.          DIMENSION FX(5),WTX(5)
    2032. 

  2032.          IF(JTYPE.EQ.2) GOTO 1
    2033. 

  2033.          WATE = WTX(LBAND)
    2034. 

  2034.          RETURN
    2035. 

  2035. 1        IF(FX(LBAND).LT.0.0001) GOTO 2
    2036. 

  2036.          WATE = WTX(LBAND)/TEMP
    2037. 

  2037.          RETURN
    2038. 

  2038. 2        WATE = WTX(LBAND)
    2039. 

  2039.       END FUNCTION wate
    2040. 

  2040. 

  2041. 

  2042. !      SUBROUTINE ERROR
    2043. 

  2043. !!         WRITE(6,*)' **** ERROR IN INPUT DATA ****'
    2044. 

  2044. !          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
    2045. 

  2045. !      END SUBROUTINE error
    2046. 

  2046. 

  2047.       
    2048. 

  2048.       SUBROUTINE REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
    2049. 

  2049. !  THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM 
    2050. 

  2050. !  FOR THE WEIGHTED CHEBCHEV APPROXIMATION OF A CONTINUOUS 
    2051. 

  2051. !  FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE 
    2052. 

  2052. !  ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
    2053. 

  2053. !  DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
    2054. 

  2054. !  GRID, THE NUMBER OF COSINES, AND THE INITIAL GUESS OF THE
    2055. 

  2055. !  EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV 
    2056. 

  2056. !  ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
    2057. 

  2057. !  FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
    2058. 

  2058. !  THE COEFFICIENTS OF THE BEST APPROXIMATION.
    2059. 

  2059. !
    2060. 

  2060. !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
    2061. 

  2061.       DIMENSION EDGE(20)
    2062. 

  2062.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
    2063. 

  2063.       DIMENSION DES(1045),GRID(1045),WT(1045)
    2064. 

  2064.       DIMENSION A(66),P(65),Q(65)
    2065. 

  2065.       DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
    2066. 

  2066.       DOUBLE PRECISION AD,DEV,X,Y
    2067. 

  2067.       DOUBLE PRECISION, EXTERNAL :: D, GEE
    2068. 

  2068. !
    2069. 

  2069. !  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
    2070. 

  2070. !
    2071. 

  2071.       ITRMAX=25
    2072. 

  2072.       DEVL=-1.0
    2073. 

  2073.       NZ=NFCNS+1
    2074. 

  2074.       NZZ=NFCNS+2
    2075. 

  2075.       NITER=0
    2076. 

  2076.   100 CONTINUE
    2077. 

  2077.       IEXT(NZZ)=NGRID+1
    2078. 

  2078.       NITER=NITER+1
    2079. 

  2079.       IF(NITER.GT.ITRMAX) GO TO 400
    2080. 

  2080.       DO 110 J=1,NZ
    2081. 

  2081.       DTEMP=GRID(IEXT(J))
    2082. 

  2082.       DTEMP=DCOS(DTEMP*PI2)
    2083. 

  2083.   110 X(J)=DTEMP
    2084. 

  2084.       JET=(NFCNS-1)/15+1
    2085. 

  2085.       DO 120 J=1,NZ
    2086. 

  2086.   120 AD(J)=D(J,NZ,JET,X)
    2087. 

  2087.       DNUM=0.0
    2088. 

  2088.       DDEN=0.0
    2089. 

  2089.       K=1
    2090. 

  2090.       DO 130 J=1,NZ
    2091. 

  2091.       L=IEXT(J)
    2092. 

  2092.       DTEMP=AD(J)*DES(L)
    2093. 

  2093.       DNUM=DNUM+DTEMP
    2094. 

  2094.       DTEMP=K*AD(J)/WT(L)
    2095. 

  2095.       DDEN=DDEN+DTEMP
    2096. 

  2096.   130 K=-K
    2097. 

  2097.       DEV=DNUM/DDEN
    2098. 

  2098.       NU=1
    2099. 

  2099.       IF(DEV.GT.0.0) NU=-1
    2100. 

  2100.       DEV=-NU*DEV
    2101. 

  2101.       K=NU
    2102. 

  2102.       DO 140 J=1,NZ
    2103. 

  2103.       L=IEXT(J)
    2104. 

  2104.       DTEMP=K*DEV/WT(L)
    2105. 

  2105.       Y(J)=DES(L)+DTEMP
    2106. 

  2106.   140 K=-K
    2107. 

  2107.       IF(DEV.GE.DEVL) GO TO 150
    2108. 

  2108.       WRITE(6,*) ' ******** FAILURE TO CONVERGE *********** '
    2109. 

  2109.       WRITE(6,*) ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR '
    2110. 

  2110.       WRITE(6,*) ' THE IMPULSE RESPONSE MAY BE CORRECT '
    2111. 

  2111.       WRITE(6,*) ' CHECK WITH A FREQUENCY RESPONSE '
    2112. 

  2112.       WRITE(6,*) ' **************************************** '
    2113. 

  2113.       GO TO 400
    2114. 

  2114.   150 DEVL=DEV
    2115. 

  2115.       JCHNGE=0
    2116. 

  2116.       K1=IEXT(1)
    2117. 

  2117.       KNZ=IEXT(NZ)
    2118. 

  2118.       KLOW=0
    2119. 

  2119.       NUT=-NU
    2120. 

  2120.       J=1
    2121. 

  2121. !
    2122. 

  2122. !  SEARCH FOR THE EXTERMAL FREQUENCIES OF THE BEST
    2123. 

  2123. !  APPROXIMATION.
    2124. 

  2124. 

  2125.   200 IF(J.EQ.NZZ) YNZ=COMP
    2126. 

  2126.       IF(J.GE.NZZ) GO TO 300
    2127. 

  2127.       KUP=IEXT(J+1)
    2128. 

  2128.       L=IEXT(J)+1
    2129. 

  2129.       NUT=-NUT
    2130. 

  2130.       IF(J.EQ.2) Y1=COMP
    2131. 

  2131.       COMP=DEV
    2132. 

  2132.       IF(L.GE.KUP) GO TO 220
    2133. 

  2133.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2134. 

  2134.       ERR=(ERR-DES(L))*WT(L)
    2135. 

  2135.       DTEMP=NUT*ERR-COMP
    2136. 

  2136.       IF(DTEMP.LE.0.0) GO TO 220
    2137. 

  2137.       COMP=NUT*ERR
    2138. 

  2138.   210 L=L+1
    2139. 

  2139.       IF(L.GE.KUP) GO TO 215
    2140. 

  2140.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2141. 

  2141.       ERR=(ERR-DES(L))*WT(L)
    2142. 

  2142.       DTEMP=NUT*ERR-COMP
    2143. 

  2143.       IF(DTEMP.LE.0.0) GO TO 215
    2144. 

  2144.       COMP=NUT*ERR
    2145. 

  2145.       GO TO 210
    2146. 

  2146.   215 IEXT(J)=L-1
    2147. 

  2147.       J=J+1
    2148. 

  2148.       KLOW=L-1
    2149. 

  2149.       JCHNGE=JCHNGE+1
    2150. 

  2150.       GO TO 200
    2151. 

  2151.   220 L=L-1
    2152. 

  2152.   225 L=L-1
    2153. 

  2153.       IF(L.LE.KLOW) GO TO 250
    2154. 

  2154.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2155. 

  2155.       ERR=(ERR-DES(L))*WT(L)
    2156. 

  2156.       DTEMP=NUT*ERR-COMP
    2157. 

  2157.       IF(DTEMP.GT.0.0) GO TO 230
    2158. 

  2158.       IF(JCHNGE.LE.0) GO TO 225
    2159. 

  2159.       GO TO 260
    2160. 

  2160.   230 COMP=NUT*ERR
    2161. 

  2161.   235 L=L-1
    2162. 

  2162.       IF(L.LE.KLOW) GO TO 240
    2163. 

  2163.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2164. 

  2164.       ERR=(ERR-DES(L))*WT(L)
    2165. 

  2165.       DTEMP=NUT*ERR-COMP
    2166. 

  2166.       IF(DTEMP.LE.0.0) GO TO 240
    2167. 

  2167.       COMP=NUT*ERR
    2168. 

  2168.       GO TO 235
    2169. 

  2169.   240 KLOW=IEXT(J)
    2170. 

  2170.       IEXT(J)=L+1
    2171. 

  2171.       J=J+1
    2172. 

  2172.       JCHNGE=JCHNGE+1
    2173. 

  2173.       GO TO 200
    2174. 

  2174.   250 L=IEXT(J)+1
    2175. 

  2175.       IF(JCHNGE.GT.0) GO TO 215
    2176. 

  2176.   255 L=L+1
    2177. 

  2177.       IF(L.GE.KUP) GO TO 260
    2178. 

  2178.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2179. 

  2179.       ERR=(ERR-DES(L))*WT(L)
    2180. 

  2180.       DTEMP=NUT*ERR-COMP
    2181. 

  2181.       IF(DTEMP.LE.0.0) GO TO 255
    2182. 

  2182.       COMP=NUT*ERR
    2183. 

  2183.       GO TO 210
    2184. 

  2184.   260 KLOW=IEXT(J)
    2185. 

  2185.       J=J+1
    2186. 

  2186.       GO TO 200
    2187. 

  2187.   300 IF(J.GT.NZZ) GO TO 320
    2188. 

  2188.       IF(K1.GT.IEXT(1)) K1=IEXT(1)
    2189. 

  2189.       IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
    2190. 

  2190.       NUT1=NUT
    2191. 

  2191.       NUT=-NU
    2192. 

  2192.       L=0
    2193. 

  2193.       KUP=K1
    2194. 

  2194.       COMP=YNZ*(1.00001)
    2195. 

  2195.       LUCK=1
    2196. 

  2196.   310 L=L+1
    2197. 

  2197.       IF(L.GE.KUP) GO TO 315
    2198. 

  2198.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2199. 

  2199.       ERR=(ERR-DES(L))*WT(L)
    2200. 

  2200.       DTEMP=NUT*ERR-COMP
    2201. 

  2201.       IF(DTEMP.LE.0.0) GO TO 310
    2202. 

  2202.       COMP=NUT*ERR
    2203. 

  2203.       J=NZZ
    2204. 

  2204.       GO TO 210
    2205. 

  2205.   315 LUCK=6
    2206. 

  2206.       GO TO 325
    2207. 

  2207.   320 IF(LUCK.GT.9) GO TO 350
    2208. 

  2208.       IF(COMP.GT.Y1) Y1=COMP
    2209. 

  2209.       K1=IEXT(NZZ)
    2210. 

  2210.   325 L=NGRID+1
    2211. 

  2211.       KLOW=KNZ
    2212. 

  2212.       NUT=-NUT1
    2213. 

  2213.       COMP=Y1*(1.00001)
    2214. 

  2214.   330 L=L-1
    2215. 

  2215.       IF(L.LE.KLOW) GO TO 340
    2216. 

  2216.       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
    2217. 

  2217.       ERR=(ERR-DES(L))*WT(L)
    2218. 

  2218.       DTEMP=NUT*ERR-COMP
    2219. 

  2219.       IF(DTEMP.LE.0.0) GO TO 330
    2220. 

  2220.       J=NZZ
    2221. 

  2221.       COMP=NUT*ERR
    2222. 

  2222.       LUCK=LUCK+10
    2223. 

  2223.       GO TO 235
    2224. 

  2224.   340 IF(LUCK.EQ.6) GO TO 370
    2225. 

  2225.       DO 345 J=1,NFCNS
    2226. 

  2226.   345 IEXT(NZZ-J)=IEXT(NZ-J)
    2227. 

  2227.       IEXT(1)=K1
    2228. 

  2228.       GO TO 100
    2229. 

  2229.   350 KN=IEXT(NZZ)
    2230. 

  2230.       DO 360 J=1,NFCNS
    2231. 

  2231.   360 IEXT(J)=IEXT(J+1)
    2232. 

  2232.       IEXT(NZ)=KN
    2233. 

  2233.       GO TO 100
    2234. 

  2234.   370 IF(JCHNGE.GT.0) GO TO 100
    2235. 

  2235. !
    2236. 

  2236. !  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
    2237. 

  2237. !  USING THE INVERSE DISCRETE FOURIER TRANSFORM.
    2238. 

  2238. !      
    2239. 

  2239.   400 CONTINUE
    2240. 

  2240.       NM1=NFCNS-1
    2241. 

  2241.       FSH=1.0E-06
    2242. 

  2242.       GTEMP=GRID(1)
    2243. 

  2243.       X(NZZ)=-2.0
    2244. 

  2244.       CN=2*NFCNS-1
    2245. 

  2245.       DELF=1.0/CN
    2246. 

  2246.       L=1
    2247. 

  2247.       KKK=0
    2248. 

  2248.       IF(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1
    2249. 

  2249.       IF(NFCNS.LE.3) KKK=1
    2250. 

  2250.       IF(KKK.EQ.1) GO TO 405
    2251. 

  2251.       DTEMP=DCOS(PI2*GRID(1))
    2252. 

  2252.       DNUM=DCOS(PI2*GRID(NGRID))
    2253. 

  2253.       AA=2.0/(DTEMP-DNUM)
    2254. 

  2254.       BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
    2255. 

  2255.   405 CONTINUE
    2256. 

  2256.       DO 430 J=1,NFCNS
    2257. 

  2257.       FT=(J-1)*DELF
    2258. 

  2258.       XT=DCOS(PI2*FT)
    2259. 

  2259.       IF(KKK.EQ.1) GO TO 410
    2260. 

  2260.       XT=(XT-BB)/AA
    2261. 

  2261. ! original :      FT=ARCOS(XT)/PI2
    2262. 

  2262.       FT=ACOS(XT)/PI2
    2263. 

  2263.   410 XE=X(L)
    2264. 

  2264.       IF(XT.GT.XE) GO TO 420
    2265. 

  2265.       IF((XE-XT).LT.FSH) GO TO 415
    2266. 

  2266.       L=L+1
    2267. 

  2267.       GO TO 410
    2268. 

  2268.   415 A(J)=Y(L)
    2269. 

  2269.       GO TO 425
    2270. 

  2270.   420 IF((XT-XE).LT.FSH) GO TO 415
    2271. 

  2271.       GRID(1)=FT
    2272. 

  2272.       A(J)=GEE(1,NZ,GRID,PI2,X,Y,AD)
    2273. 

  2273.   425 CONTINUE
    2274. 

  2274.       IF(L.GT.1) L=L-1
    2275. 

  2275.   430 CONTINUE
    2276. 

  2276.       GRID(1)=GTEMP
    2277. 

  2277.       DDEN=PI2/CN
    2278. 

  2278.       DO 510 J=1,NFCNS
    2279. 

  2279.       DTEMP=0.0
    2280. 

  2280.       DNUM=(J-1)*DDEN
    2281. 

  2281.       IF(NM1.LT.1) GO TO 505
    2282. 

  2282.       DO 500 K=1,NM1
    2283. 

  2283.   500 DTEMP=DTEMP+A(K+1)*DCOS(DNUM*K)
    2284. 

  2284.   505 DTEMP=2.0*DTEMP+A(1)
    2285. 

  2285.   510 ALPHA(J)=DTEMP
    2286. 

  2286.       DO 550 J=2,NFCNS
    2287. 

  2287.   550 ALPHA(J)=2*ALPHA(J)/CN
    2288. 

  2288.       ALPHA(1)=ALPHA(1)/CN
    2289. 

  2289.       IF(KKK.EQ.1) GO TO 545
    2290. 

  2290.       P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
    2291. 

  2291.       P(2)=2.0*AA*ALPHA(NFCNS)
    2292. 

  2292.       Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
    2293. 

  2293.       DO 540 J=2,NM1
    2294. 

  2294.       IF(J.LT.NM1) GO TO 515
    2295. 

  2295.       AA=0.5*AA
    2296. 

  2296.       BB=0.5*BB
    2297. 

  2297.   515 CONTINUE
    2298. 

  2298.       P(J+1)=0.0
    2299. 

  2299.       DO 520 K=1,J
    2300. 

  2300.       A(K)=P(K)
    2301. 

  2301.   520 P(K)=2.0*BB*A(K)
    2302. 

  2302.       P(2)=P(2)+A(1)*2.0*AA
    2303. 

  2303.       JM1=J-1
    2304. 

  2304.       DO 525 K=1,JM1
    2305. 

  2305.   525 P(K)=P(K)+Q(K)+AA*A(K+1)
    2306. 

  2306.       JP1=J+1
    2307. 

  2307.       DO 530 K=3,JP1
    2308. 

  2308.   530 P(K)=P(K)+AA*A(K-1)
    2309. 

  2309.       IF(J.EQ.NM1) GO TO 540
    2310. 

  2310.       DO 535 K=1,J
    2311. 

  2311.   535 Q(K)=-A(K)
    2312. 

  2312.       Q(1)=Q(1)+ALPHA(NFCNS-1-J)
    2313. 

  2313.   540 CONTINUE
    2314. 

  2314.       DO 543 J=1,NFCNS
    2315. 

  2315.   543 ALPHA(J)=P(J)
    2316. 

  2316.   545 CONTINUE
    2317. 

  2317.       IF(NFCNS.GT.3) RETURN
    2318. 

  2318.       ALPHA(NFCNS+1)=0.0
    2319. 

  2319.       ALPHA(NFCNS+2)=0.0
    2320. 

  2320.    END SUBROUTINE remez
    2321. 

  2321. 

  2322.    DOUBLE PRECISION FUNCTION D(K,N,M,X)
    2323. 

  2323. !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
    2324. 

  2324.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
    2325. 

  2325.       DIMENSION DES(1045),GRID(1045),WT(1045)
    2326. 

  2326.       DOUBLE PRECISION AD,DEV,X,Y
    2327. 

  2327.       DOUBLE PRECISION Q
    2328. 

  2328.       DOUBLE PRECISION PI2
    2329. 

  2329.       D = 1.0
    2330. 

  2330.       Q = X(K)
    2331. 

  2331.       DO 3 L = 1,M
    2332. 

  2332.       DO 2 J = L,N,M
    2333. 

  2333.             IF(J-K) 1,2,1
    2334. 

  2334. 1                  D = 2.0*D*(Q-X(J))
    2335. 

  2335. 2      CONTINUE
    2336. 

  2336. 3      CONTINUE
    2337. 

  2337.       D = 1.0/D
    2338. 

  2338.    END FUNCTION D
    2339. 

  2339. 

  2340. 

  2341.    DOUBLE PRECISION FUNCTION GEE(K,N,GRID,PI2,X,Y,AD)
    2342. 

  2342. !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
    2343. 

  2343.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
    2344. 

  2344.       DIMENSION DES(1045),GRID(1045),WT(1045)
    2345. 

  2345.       DOUBLE PRECISION AD,DEV,X,Y
    2346. 

  2346.       DOUBLE PRECISION P,C,D,XF
    2347. 

  2347.       DOUBLE PRECISION PI2
    2348. 

  2348.       P = 0.0
    2349. 

  2349.       XF = GRID(K)
    2350. 

  2350.       XF = DCOS(PI2*XF)
    2351. 

  2351.       D = 0.0
    2352. 

  2352.       DO 1 J =1,N
    2353. 

  2353.             C = XF-X(J)
    2354. 

  2354.             C = AD(J)/C
    2355. 

  2355.             D = D+C
    2356. 

  2356.             P = P+C*Y(J)
    2357. 

  2357. 1      CONTINUE
    2358. 

  2358.       GEE = P/D
    2359. 

  2359.    END FUNCTION GEE
    2360. 

  2360. 

  2361.       REAL FUNCTION RSLF(P,T)
    2362. 

  2362. 

  2363.       IMPLICIT NONE
    2364. 

  2364.       REAL, INTENT(IN):: P, T
    2365. 

  2365.       REAL:: ESL,X
    2366. 

  2366.       REAL, PARAMETER:: C0= .611583699E03
    2367. 

  2367.       REAL, PARAMETER:: C1= .444606896E02
    2368. 

  2368.       REAL, PARAMETER:: C2= .143177157E01
    2369. 

  2369.       REAL, PARAMETER:: C3= .264224321E-1
    2370. 

  2370.       REAL, PARAMETER:: C4= .299291081E-3
    2371. 

  2371.       REAL, PARAMETER:: C5= .203154182E-5
    2372. 

  2372.       REAL, PARAMETER:: C6= .702620698E-8
    2373. 

  2373.       REAL, PARAMETER:: C7= .379534310E-11
    2374. 

  2374.       REAL, PARAMETER:: C8=-.321582393E-13
    2375. 

  2375. 

  2376.       X=MAX(-80.,T-273.16)
    2377. 

  2377. 

  2378. !      ESL=612.2*EXP(17.67*X/(T-29.65))
    2379. 

  2379.       ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
    2380. 

  2380.       RSLF=.622*ESL/(P-ESL)
    2381. 

  2381. 

  2382.       END FUNCTION RSLF
    2383. 

  2383. 

  2384. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2385. 

  2385. ! DFI startfwd group of functions
    2386. 

  2386. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2387. 

  2387. 

  2388.    SUBROUTINE wrf_dfi_startfwd_init ( )
    2389. 

  2389. 

  2390.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
    2391. 

  2391.      USE module_utility
    2392. 

  2392.      
    2393. 

  2393.      IMPLICIT NONE
    2394. 

  2394. 

  2395.      INTERFACE
    2396. 

  2396.         SUBROUTINE dfi_startfwd_init_recurse(grid)
    2397. 

  2397.           USE module_domain, ONLY : domain
    2398. 

  2398.           TYPE (domain), POINTER :: grid
    2399. 

  2399.         END SUBROUTINE dfi_startfwd_init_recurse
    2400. 

  2400.      END INTERFACE
    2401. 

  2401. 

  2402.      ! Now, setup all nests
    2403. 

  2403. 

  2404.      CALL dfi_startfwd_init_recurse(head_grid)
    2405. 

  2405.      
    2406. 

  2406.      CALL set_current_grid_ptr( head_grid )
    2407. 

  2407. 

  2408.    END SUBROUTINE wrf_dfi_startfwd_init
    2409. 

  2409. 

  2410. 

  2411.    RECURSIVE SUBROUTINE dfi_startfwd_init_recurse(grid)
    2412. 

  2412. 

  2413.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
    2414. 

  2414. 

  2415.      IMPLICIT NONE
    2416. 

  2416. 

  2417.       INTERFACE
    2418. 

  2418.          SUBROUTINE dfi_startfwd_init(grid)
    2419. 

  2419.             USE module_domain, ONLY : domain
    2420. 

  2420.             TYPE (domain), POINTER :: grid
    2421. 

  2421.          END SUBROUTINE dfi_startfwd_init
    2422. 

  2422.       END INTERFACE
    2423. 

  2423. 

  2424.      INTEGER :: kid
    2425. 

  2425.      TYPE (domain), POINTER :: grid
    2426. 

  2426.      TYPE (domain), POINTER :: grid_ptr
    2427. 

  2427.     
    2428. 

  2428.      grid_ptr => grid
    2429. 

  2429. 

  2430.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2431. 

  2431.         !
    2432. 

  2432.         ! Assure that time-step is set back to positive
    2433. 

  2433.         !   for this forward step.
    2434. 

  2434.         !
    2435. 

  2435.         grid_ptr%dt = abs(grid_ptr%dt)
    2436. 

  2436.         grid_ptr%time_step = abs(grid_ptr%time_step)
    2437. 

  2437.         CALL set_current_grid_ptr( grid_ptr )
    2438. 

  2438.         CALL dfi_startfwd_init( grid_ptr )
    2439. 

  2439.         DO kid = 1, max_nests
    2440. 

  2440.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2441. 

  2441.               CALL dfi_startfwd_init_recurse(grid_ptr%nests(kid)%ptr)
    2442. 

  2442.            ENDIF
    2443. 

  2443.         END DO
    2444. 

  2444.         grid_ptr => grid_ptr%sibling
    2445. 

  2445.      END DO
    2446. 

  2446.      
    2447. 

  2447.    END SUBROUTINE dfi_startfwd_init_recurse
    2448. 

  2448. 

  2449. 

  2450.    SUBROUTINE dfi_startfwd_init ( grid )
    2451. 

  2451. 

  2452.       USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
    2453. 

  2453.       USE module_utility
    2454. 

  2454.       USE module_state_description
    2455. 

  2455. 

  2456.       IMPLICIT NONE
    2457. 

  2457. 

  2458.       TYPE (domain) , POINTER                 :: grid
    2459. 

  2459.       INTEGER rc
    2460. 

  2460. 

  2461.       INTERFACE
    2462. 

  2462.          SUBROUTINE Setup_Timekeeping(grid)
    2463. 

  2463.             USE module_domain, ONLY : domain
    2464. 

  2464.             TYPE (domain), POINTER :: grid
    2465. 

  2465.          END SUBROUTINE Setup_Timekeeping
    2466. 

  2466. 

  2467.       END INTERFACE
    2468. 

  2468. 

  2469.       grid%dfi_stage = DFI_STARTFWD
    2470. 

  2470. 

  2471. #if (EM_CORE == 1)
    2472. 

  2472.       ! No need for adaptive time-step
    2473. 

  2473.       CALL nl_set_use_adaptive_time_step( grid%id, .false. )
    2474. 

  2474. #endif
    2475. 

  2475. 

  2476.       CALL Setup_Timekeeping (grid)
    2477. 

  2477.       grid%start_subtime = domain_get_start_time ( head_grid )
    2478. 

  2478.       grid%stop_subtime = domain_get_stop_time ( head_grid )
    2479. 

  2479. 

  2480.       CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
    2481. 

  2481. 

  2482.    END SUBROUTINE dfi_startfwd_init
    2483. 

  2483. 

  2484. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2485. 

  2485. ! DFI startbck group of functions
    2486. 

  2486. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2487. 

  2487. 

  2488.    SUBROUTINE wrf_dfi_startbck_init ( )
    2489. 

  2489. 

  2490.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
    2491. 

  2491.      USE module_utility
    2492. 

  2492.      
    2493. 

  2493.      IMPLICIT NONE
    2494. 

  2494. 

  2495.      INTERFACE
    2496. 

  2496.         SUBROUTINE dfi_startbck_init_recurse(grid)
    2497. 

  2497.           USE module_domain, ONLY : domain
    2498. 

  2498.           TYPE (domain), POINTER :: grid
    2499. 

  2499.         END SUBROUTINE dfi_startbck_init_recurse
    2500. 

  2500.      END INTERFACE
    2501. 

  2501. 

  2502.      ! Now, setup all nests
    2503. 

  2503. 

  2504.      CALL dfi_startbck_init_recurse(head_grid)
    2505. 

  2505.      
    2506. 

  2506.      CALL set_current_grid_ptr( head_grid )
    2507. 

  2507. 

  2508.    END SUBROUTINE wrf_dfi_startbck_init
    2509. 

  2509. 

  2510. 

  2511.    RECURSIVE SUBROUTINE dfi_startbck_init_recurse(grid)
    2512. 

  2512. 

  2513.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
    2514. 

  2514. 

  2515.      IMPLICIT NONE
    2516. 

  2516. 

  2517.       INTERFACE
    2518. 

  2518.          SUBROUTINE dfi_startbck_init(grid)
    2519. 

  2519.             USE module_domain, ONLY : domain
    2520. 

  2520.             TYPE (domain), POINTER :: grid
    2521. 

  2521.          END SUBROUTINE dfi_startbck_init
    2522. 

  2522.       END INTERFACE
    2523. 

  2523. 

  2524.      INTEGER :: kid
    2525. 

  2525.      TYPE (domain), POINTER :: grid
    2526. 

  2526.      TYPE (domain), POINTER :: grid_ptr
    2527. 

  2527.     
    2528. 

  2528.      grid_ptr => grid
    2529. 

  2529. 

  2530.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2531. 

  2531.         !
    2532. 

  2532.         ! Assure that time-step is set back to positive
    2533. 

  2533.         !   for this forward step.
    2534. 

  2534.         !
    2535. 

  2535.         grid_ptr%dt = abs(grid_ptr%dt)
    2536. 

  2536.         grid_ptr%time_step = abs(grid_ptr%time_step)
    2537. 

  2537.         CALL set_current_grid_ptr( grid_ptr )
    2538. 

  2538.         CALL dfi_startbck_init( grid_ptr )
    2539. 

  2539.         DO kid = 1, max_nests
    2540. 

  2540.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2541. 

  2541.               CALL dfi_startbck_init_recurse(grid_ptr%nests(kid)%ptr)
    2542. 

  2542.            ENDIF
    2543. 

  2543.         END DO
    2544. 

  2544.         grid_ptr => grid_ptr%sibling
    2545. 

  2545.      END DO
    2546. 

  2546.      
    2547. 

  2547.    END SUBROUTINE dfi_startbck_init_recurse
    2548. 

  2548. 

  2549. 

  2550.    SUBROUTINE dfi_startbck_init ( grid )
    2551. 

  2551. 

  2552.       USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
    2553. 

  2553.       USE module_utility
    2554. 

  2554.       USE module_state_description
    2555. 

  2555. 

  2556.       IMPLICIT NONE
    2557. 

  2557. 

  2558.       TYPE (domain) , POINTER                 :: grid
    2559. 

  2559.       INTEGER rc
    2560. 

  2560. 

  2561.       INTERFACE
    2562. 

  2562.          SUBROUTINE Setup_Timekeeping(grid)
    2563. 

  2563.             USE module_domain, ONLY : domain
    2564. 

  2564.             TYPE (domain), POINTER :: grid
    2565. 

  2565.          END SUBROUTINE Setup_Timekeeping
    2566. 

  2566. 

  2567.       END INTERFACE
    2568. 

  2568. 

  2569.       grid%dfi_stage = DFI_STARTBCK
    2570. 

  2570. 

  2571.       ! set physics options to zero
    2572. 

  2572.       CALL nl_set_mp_physics( grid%id, 0 )
    2573. 

  2573.       CALL nl_set_ra_lw_physics( grid%id, 0 )
    2574. 

  2574.       CALL nl_set_ra_sw_physics( grid%id, 0 )
    2575. 

  2575.       CALL nl_set_sf_surface_physics( grid%id, 0 )
    2576. 

  2576.       CALL nl_set_sf_sfclay_physics( grid%id, 0 )
    2577. 

  2577.       CALL nl_set_sf_urban_physics( grid%id, 0 )
    2578. 

  2578.       CALL nl_set_bl_pbl_physics( grid%id, 0 )
    2579. 

  2579.       CALL nl_set_cu_physics( grid%id, 0 )
    2580. 

  2580.       CALL nl_set_damp_opt( grid%id, 0 )
    2581. 

  2581.       CALL nl_set_sst_update( grid%id, 0 )
    2582. 

  2582.       CALL nl_set_gwd_opt( grid%id, 0 )
    2583. 

  2583. #if (EM_CORE == 1)
    2584. 

  2584.       CALL nl_set_diff_6th_opt( grid%id, 0 )
    2585. 

  2585.       CALL nl_set_use_adaptive_time_step( grid%id, .false. )
    2586. 

  2586. #endif
    2587. 

  2587.       CALL nl_set_feedback( grid%id, 0 )
    2588. 

  2588. #if (EM_CORE == 1)
    2589. 

  2589.       ! set bc
    2590. 

  2590.       CALL nl_set_constant_bc( grid%id, head_grid%constant_bc)
    2591. 

  2591. #endif
    2592. 

  2592. 

  2593. #ifdef WRF_CHEM
    2594. 

  2594.       ! set chemistry option to zero
    2595. 

  2595.       CALL nl_set_chem_opt (grid%id, 0)
    2596. 

  2596.       CALL nl_set_aer_ra_feedback (grid%id, 0)
    2597. 

  2597.       CALL nl_set_io_form_auxinput5 (grid%id, 0)
    2598. 

  2598.       CALL nl_set_io_form_auxinput7 (grid%id, 0)
    2599. 

  2599.       CALL nl_set_io_form_auxinput8 (grid%id, 0)
    2600. 

  2600. #endif
    2601. 

  2601. 

  2602. #if (EM_CORE == 1)
    2603. 

  2603.       ! set diffusion to zero for backward integration
    2604. 

  2604.       CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
    2605. 

  2605.       CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
    2606. 

  2606.       IF ( grid%moist_adv_opt == 2 ) THEN
    2607. 

  2607.          CALL nl_set_moist_adv_opt( grid%id, 0)
    2608. 

  2608.       ENDIF
    2609. 

  2609. #endif
    2610. 

  2610. 

  2611.       !tgs need to call start_domain here to reset bc initialization for 
    2612. 

  2612.       !      negative dt, but only for outer domain.
    2613. 

  2613.       if (grid%id == 1) then
    2614. 

  2614.         CALL start_domain ( grid , .TRUE. )
    2615. 

  2615.       endif
    2616. 

  2616. 

  2617.       ! Call wrf_run to advance forward 1 step
    2618. 

  2618. 

  2619.       ! Negate time step
    2620. 

  2620.       CALL nl_set_time_step ( grid%id, -grid%time_step )
    2621. 

  2621. 

  2622.       CALL Setup_Timekeeping (grid) 
    2623. 

  2623. 

  2624.       grid%start_subtime = domain_get_start_time ( grid )
    2625. 

  2625.       grid%stop_subtime = domain_get_stop_time ( grid )
    2626. 

  2626. 

  2627.       CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
    2628. 

  2628. 

  2629.    END SUBROUTINE dfi_startbck_init
    2630. 

  2630. 

  2631. 

  2632.    SUBROUTINE wrf_dfi_bck_init ( )
    2633. 

  2633. 

  2634.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
    2635. 

  2635.      USE module_utility
    2636. 

  2636.      USE module_state_description
    2637. 

  2637.      
    2638. 

  2638.      IMPLICIT NONE
    2639. 

  2639. 

  2640.      INTERFACE
    2641. 

  2641.         SUBROUTINE dfi_bck_init_recurse(grid)
    2642. 

  2642.           USE module_domain, ONLY : domain
    2643. 

  2643.           TYPE (domain), POINTER :: grid
    2644. 

  2644.         END SUBROUTINE dfi_bck_init_recurse
    2645. 

  2645.      END INTERFACE
    2646. 

  2646. 

  2647.      ! We can only call dfi_bck_init for the head_grid
    2648. 

  2648.      !    since nests have not been instantiated at this point,
    2649. 

  2649.      !    so, dfi_bck_init will need to be called for each
    2650. 

  2650.      !    nest from integrate.
    2651. 

  2651.      CALL dfi_bck_init_recurse(head_grid)
    2652. 

  2652. 

  2653.    END SUBROUTINE wrf_dfi_bck_init
    2654. 

  2654. 

  2655.    RECURSIVE SUBROUTINE dfi_bck_init_recurse(grid)
    2656. 

  2656. 

  2657.      USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
    2658. 

  2658. 

  2659.      IMPLICIT NONE
    2660. 

  2660. 

  2661.       INTERFACE
    2662. 

  2662.          SUBROUTINE dfi_bck_init(grid)
    2663. 

  2663.             USE module_domain, ONLY : domain
    2664. 

  2664.             TYPE (domain), POINTER :: grid
    2665. 

  2665.          END SUBROUTINE dfi_bck_init
    2666. 

  2666.       END INTERFACE
    2667. 

  2667. 

  2668.      INTEGER :: kid
    2669. 

  2669.      TYPE (domain), POINTER :: grid
    2670. 

  2670.      TYPE (domain), POINTER :: grid_ptr
    2671. 

  2671.     
    2672. 

  2672.      grid_ptr => grid
    2673. 

  2673. 

  2674.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2675. 

  2675.         !
    2676. 

  2676.         ! Assure that time-step is set back to positive
    2677. 

  2677.         !   for this forward step.
    2678. 

  2678.         !
    2679. 

  2679.         grid_ptr%dt = abs(grid_ptr%dt)
    2680. 

  2680.         grid_ptr%time_step = abs(grid_ptr%time_step)
    2681. 

  2681.         CALL set_current_grid_ptr( grid_ptr )
    2682. 

  2682.         CALL dfi_bck_init( grid_ptr )
    2683. 

  2683.         DO kid = 1, max_nests
    2684. 

  2684.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2685. 

  2685.               CALL dfi_bck_init_recurse(grid_ptr%nests(kid)%ptr)
    2686. 

  2686.            ENDIF
    2687. 

  2687.         END DO
    2688. 

  2688.         grid_ptr => grid_ptr%sibling
    2689. 

  2689.      END DO
    2690. 

  2690.      
    2691. 

  2691.    END SUBROUTINE dfi_bck_init_recurse
    2692. 

  2692. 

  2693. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2694. 

  2694. ! DFI forward initialization group of functions
    2695. 

  2695. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2696. 

  2696. 

  2697.    SUBROUTINE wrf_dfi_fwd_init ( )
    2698. 

  2698. 

  2699.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
    2700. 

  2700.      USE module_utility
    2701. 

  2701.      
    2702. 

  2702.      IMPLICIT NONE
    2703. 

  2703. 

  2704.      INTERFACE
    2705. 

  2705.         SUBROUTINE dfi_fwd_init_recurse(grid)
    2706. 

  2706.           USE module_domain, ONLY : domain
    2707. 

  2707.           TYPE (domain), POINTER :: grid
    2708. 

  2708.         END SUBROUTINE dfi_fwd_init_recurse
    2709. 

  2709.      END INTERFACE
    2710. 

  2710. 

  2711.      ! Now, setup all nests
    2712. 

  2712. 

  2713.      CALL dfi_fwd_init_recurse(head_grid)
    2714. 

  2714.      
    2715. 

  2715.      CALL set_current_grid_ptr( head_grid )
    2716. 

  2716. 

  2717.    END SUBROUTINE wrf_dfi_fwd_init
    2718. 

  2718. 

  2719.    RECURSIVE SUBROUTINE dfi_fwd_init_recurse(grid)
    2720. 

  2720. 

  2721.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
    2722. 

  2722. 

  2723.      IMPLICIT NONE
    2724. 

  2724. 

  2725.       INTERFACE
    2726. 

  2726.          SUBROUTINE dfi_fwd_init(grid)
    2727. 

  2727.             USE module_domain, ONLY : domain
    2728. 

  2728.             TYPE (domain), POINTER :: grid
    2729. 

  2729.          END SUBROUTINE dfi_fwd_init
    2730. 

  2730.       END INTERFACE
    2731. 

  2731. 

  2732.      INTEGER :: kid
    2733. 

  2733.      TYPE (domain), POINTER :: grid
    2734. 

  2734.      TYPE (domain), POINTER :: grid_ptr
    2735. 

  2735.     
    2736. 

  2736.      grid_ptr => grid
    2737. 

  2737. 

  2738.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2739. 

  2739.         !
    2740. 

  2740.         ! Assure that time-step is set back to positive
    2741. 

  2741.         !   for this forward step.
    2742. 

  2742.         !
    2743. 

  2743.         grid_ptr%dt = abs(grid_ptr%dt)
    2744. 

  2744.         grid_ptr%time_step = abs(grid_ptr%time_step)
    2745. 

  2745.         CALL set_current_grid_ptr( grid_ptr )
    2746. 

  2746.         CALL dfi_fwd_init( grid_ptr )
    2747. 

  2747.         DO kid = 1, max_nests
    2748. 

  2748.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2749. 

  2749.               CALL dfi_fwd_init_recurse(grid_ptr%nests(kid)%ptr)
    2750. 

  2750.            ENDIF
    2751. 

  2751.         END DO
    2752. 

  2752.         grid_ptr => grid_ptr%sibling
    2753. 

  2753.      END DO
    2754. 

  2754.      
    2755. 

  2755.    END SUBROUTINE dfi_fwd_init_recurse
    2756. 

  2756. 

  2757. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2758. 

  2758. ! DFI forecast initialization group of functions
    2759. 

  2759. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2760. 

  2760. 

  2761.    SUBROUTINE wrf_dfi_fst_init ( )
    2762. 

  2762. 

  2763.      USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
    2764. 

  2764.      USE module_utility
    2765. 

  2765.      
    2766. 

  2766.      IMPLICIT NONE
    2767. 

  2767. 

  2768.      INTERFACE
    2769. 

  2769.         SUBROUTINE dfi_fst_init_recurse(grid)
    2770. 

  2770.           USE module_domain, ONLY : domain
    2771. 

  2771.           TYPE (domain), POINTER :: grid
    2772. 

  2772.         END SUBROUTINE dfi_fst_init_recurse
    2773. 

  2773.      END INTERFACE
    2774. 

  2774. 

  2775.      ! Now, setup all nests
    2776. 

  2776. 

  2777.      CALL dfi_fst_init_recurse(head_grid)
    2778. 

  2778.      
    2779. 

  2779.      CALL set_current_grid_ptr( head_grid )
    2780. 

  2780. 

  2781.    END SUBROUTINE wrf_dfi_fst_init
    2782. 

  2782. 

  2783.    RECURSIVE SUBROUTINE dfi_fst_init_recurse ( grid )
    2784. 

  2784. 

  2785.      USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
    2786. 

  2786. 

  2787.      IMPLICIT NONE
    2788. 

  2788. 

  2789.       INTERFACE
    2790. 

  2790.          SUBROUTINE dfi_fst_init(grid)
    2791. 

  2791.             USE module_domain, ONLY : domain
    2792. 

  2792.             TYPE (domain), POINTER :: grid
    2793. 

  2793.          END SUBROUTINE dfi_fst_init
    2794. 

  2794.       END INTERFACE
    2795. 

  2795. 

  2796.      INTEGER :: kid
    2797. 

  2797.      TYPE (domain), POINTER :: grid
    2798. 

  2798.      TYPE (domain), POINTER :: grid_ptr
    2799. 

  2799.     
    2800. 

  2800.      grid_ptr => grid
    2801. 

  2801. 

  2802.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2803. 

  2803.         !
    2804. 

  2804.         ! Assure that time-step is set back to positive
    2805. 

  2805.         !   for this forward step.
    2806. 

  2806.         !
    2807. 

  2807.         grid_ptr%dt = abs(grid_ptr%dt)
    2808. 

  2808.         grid_ptr%time_step = abs(grid_ptr%time_step)
    2809. 

  2809.         CALL set_current_grid_ptr( grid_ptr )
    2810. 

  2810.         CALL dfi_fst_init( grid_ptr )
    2811. 

  2811.         DO kid = 1, max_nests
    2812. 

  2812.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2813. 

  2813.               CALL dfi_fst_init_recurse(grid_ptr%nests(kid)%ptr)
    2814. 

  2814.            ENDIF
    2815. 

  2815.         END DO
    2816. 

  2816.         grid_ptr => grid_ptr%sibling
    2817. 

  2817.      END DO
    2818. 

  2818.      
    2819. 

  2819.    END SUBROUTINE dfi_fst_init_recurse
    2820. 

  2820. 

  2821. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2822. 

  2822. ! DFI write initialization group of functions
    2823. 

  2823. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2824. 

  2824. 

  2825.    SUBROUTINE wrf_dfi_write_initialized_state( )
    2826. 

  2826. 

  2827.      USE module_domain, ONLY : domain, head_grid
    2828. 

  2828. 

  2829.      INTERFACE
    2830. 

  2830.         SUBROUTINE dfi_write_initialized_state_recurse(grid)
    2831. 

  2831.           USE module_domain, ONLY : domain
    2832. 

  2832.           TYPE (domain), POINTER :: grid
    2833. 

  2833.         END SUBROUTINE dfi_write_initialized_state_recurse
    2834. 

  2834.      END INTERFACE
    2835. 

  2835. 

  2836.      ! Now, setup all nests
    2837. 

  2837. 

  2838.      CALL dfi_write_initialized_state_recurse(head_grid)
    2839. 

  2839.      
    2840. 

  2840.    END SUBROUTINE wrf_dfi_write_initialized_state
    2841. 

  2841. 

  2842.    RECURSIVE SUBROUTINE dfi_write_initialized_state_recurse( grid )
    2843. 

  2843. 

  2844.      USE module_domain, ONLY : domain, max_nests
    2845. 

  2845.      
    2846. 

  2846.      IMPLICIT NONE
    2847. 

  2847.      
    2848. 

  2848.      INTERFACE
    2849. 

  2849.         SUBROUTINE dfi_write_initialized_state( grid )
    2850. 

  2850.           USE module_domain, ONLY : domain
    2851. 

  2851.           TYPE (domain), POINTER :: grid
    2852. 

  2852.         END SUBROUTINE dfi_write_initialized_state
    2853. 

  2853.      END INTERFACE
    2854. 

  2854.      
    2855. 

  2855.      INTEGER :: kid
    2856. 

  2856.      TYPE (domain), POINTER :: grid
    2857. 

  2857.      TYPE (domain), POINTER :: grid_ptr
    2858. 

  2858.      
    2859. 

  2859.      grid_ptr => grid
    2860. 

  2860.      
    2861. 

  2861.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2862. 

  2862.         !
    2863. 

  2863.         ! Assure that time-step is set back to positive
    2864. 

  2864.         !   for this forward step.
    2865. 

  2865.         !
    2866. 

  2866.         CALL dfi_write_initialized_state( grid_ptr )
    2867. 

  2867.         DO kid = 1, max_nests
    2868. 

  2868.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2869. 

  2869.               CALL dfi_write_initialized_state_recurse(grid_ptr%nests(kid)%ptr)
    2870. 

  2870.            ENDIF
    2871. 

  2871.         END DO
    2872. 

  2872.         grid_ptr => grid_ptr%sibling
    2873. 

  2873.      END DO
    2874. 

  2874.      
    2875. 

  2875.    END SUBROUTINE dfi_write_initialized_state_recurse
    2876. 

  2876. 

  2877. 

  2878.    RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid)
    2879. 

  2879. 

  2880.      USE module_domain, ONLY : domain, max_nests, set_current_grid_ptr
    2881. 

  2881. 

  2882.      IMPLICIT NONE
    2883. 

  2883. 

  2884.       INTERFACE
    2885. 

  2885.          SUBROUTINE dfi_array_reset(grid)
    2886. 

  2886.             USE module_domain, ONLY : domain
    2887. 

  2887.             TYPE (domain), POINTER :: grid
    2888. 

  2888.          END SUBROUTINE dfi_array_reset
    2889. 

  2889.       END INTERFACE
    2890. 

  2890. 

  2891.      INTEGER :: kid
    2892. 

  2892.      TYPE (domain), POINTER :: grid
    2893. 

  2893.      TYPE (domain), POINTER :: grid_ptr
    2894. 

  2894.     
    2895. 

  2895.      grid_ptr => grid
    2896. 

  2896. 

  2897.      DO WHILE ( ASSOCIATED( grid_ptr ) )
    2898. 

  2898.         CALL set_current_grid_ptr( grid_ptr )
    2899. 

  2899.         CALL dfi_array_reset( grid_ptr )
    2900. 

  2900.         DO kid = 1, max_nests
    2901. 

  2901.            IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
    2902. 

  2902.               CALL dfi_array_reset_recurse(grid_ptr%nests(kid)%ptr)
    2903. 

  2903.            ENDIF
    2904. 

  2904.         END DO
    2905. 

  2905.         grid_ptr => grid_ptr%sibling
    2906. 

  2906.      END DO
    2907. 

  2907.      
    2908. 

  2908.    END SUBROUTINE dfi_array_reset_recurse
    2909. 

  2909. 

  2910.