PageRenderTime 58ms CodeModel.GetById 19ms RepoModel.GetById 0ms app.codeStats 1ms

/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

Large files files are truncated, but you can click here to view the full file

    

  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.

Large files files are truncated, but you can click here to view the full file