wrfv2_fire/share/dfi.F

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