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