SUBROUTINE dfi_accumulate( grid ) 2,5
USE module_domain
, ONLY : domain
! USE module_configure
USE module_driver_constants
USE module_machine
! USE module_dm
USE module_model_constants
USE module_state_description
IMPLICIT NONE
REAL hn
CHARACTER*80 mess
! Input data.
TYPE(domain) , POINTER :: grid
IF ( grid%dfi_opt .EQ. DFI_NODFI .OR. grid%dfi_stage .EQ. DFI_FST ) RETURN
#if (EM_CORE == 1)
hn = grid%hcoeff(grid%itimestep+1)
! accumulate dynamic variables
grid%dfi_mu(:,:) = grid%dfi_mu(:,:) + grid%mu_2(:,:) * hn
grid%dfi_u(:,:,:) = grid%dfi_u(:,:,:) + grid%u_2(:,:,:) * hn
grid%dfi_v(:,:,:) = grid%dfi_v(:,:,:) + grid%v_2(:,:,:) * hn
grid%dfi_w(:,:,:) = grid%dfi_w(:,:,:) + grid%w_2(:,:,:) * hn
grid%dfi_ww(:,:,:) = grid%dfi_ww(:,:,:) + grid%ww(:,:,:) * hn
grid%dfi_t(:,:,:) = grid%dfi_t(:,:,:) + grid%t_2(:,:,:) * hn
grid%dfi_phb(:,:,:) = grid%dfi_phb(:,:,:) + grid%phb(:,:,:) * hn
grid%dfi_ph0(:,:,:) = grid%dfi_ph0(:,:,:) + grid%ph0(:,:,:) * hn
grid%dfi_php(:,:,:) = grid%dfi_php(:,:,:) + grid%php(:,:,:) * hn
grid%dfi_p(:,:,:) = grid%dfi_p(:,:,:) + grid%p(:,:,:) * hn
grid%dfi_ph(:,:,:) = grid%dfi_ph(:,:,:) + grid%ph_2(:,:,:) * hn
grid%dfi_tke(:,:,:) = grid%dfi_tke(:,:,:) + grid%tke_2(:,:,:) * hn
grid%dfi_al(:,:,:) = grid%dfi_al(:,:,:) + grid%al(:,:,:) * hn
grid%dfi_alt(:,:,:) = grid%dfi_alt(:,:,:) + grid%alt(:,:,:) * hn
grid%dfi_pb(:,:,:) = grid%dfi_pb(:,:,:) + grid%pb(:,:,:) * hn
! neg. check on hydrometeor and scalar variables
grid%moist(:,:,:,:) = max(0.,grid%moist(:,:,:,:))
grid%dfi_scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:))
! dfi_savehydmeteors is a namelist parameter, default =0 which means hydrometeor
! and scalar fields will be spinning up in DFI; if dfi_savehydmeteors=1 then
! hydrometeor fields will stay unchanged in DFI, but water vapor mixing ratio
! QV will be modified.
IF ( grid%dfi_savehydmeteors .EQ. 0 ) then
grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
ELSE
grid%dfi_moist(:,:,:,P_QV) = grid%dfi_moist(:,:,:,P_QV) + grid%moist(:,:,:,P_QV) * hn
ENDIF
! IF ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
! grid%dfi_QVG(:,:) = grid%dfi_QVG(:,:) + grid%QVG(:,:) * hn
! ENDIF
! accumulate DFI coefficient
grid%hcoeff_tot = grid%hcoeff_tot + hn
#endif
#if (NMM_CORE == 1)
hn = grid%hcoeff(grid%ntsd+1)
grid%dfi_pd(:,:) = grid%dfi_pd(:,:) + grid%pd(:,:) * hn
grid%dfi_pint(:,:,:) = grid%dfi_pint(:,:,:) + grid%pint(:,:,:) * hn
grid%dfi_dwdt(:,:,:) = grid%dfi_dwdt(:,:,:) + grid%dwdt(:,:,:) * hn
grid%dfi_t(:,:,:) = grid%dfi_t(:,:,:) + grid%t(:,:,:) * hn
grid%dfi_q(:,:,:) = grid%dfi_q(:,:,:) + grid%q(:,:,:) * hn
grid%dfi_q2(:,:,:) = grid%dfi_q2(:,:,:) + grid%q2(:,:,:) * hn
grid%dfi_cwm(:,:,:) = grid%dfi_cwm(:,:,:) + grid%cwm(:,:,:) * hn
grid%dfi_u(:,:,:) = grid%dfi_u(:,:,:) + grid%u(:,:,:) * hn
grid%dfi_v(:,:,:) = grid%dfi_v(:,:,:) + grid%v(:,:,:) * hn
grid%dfi_moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) + grid%moist(:,:,:,:) * hn
grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
! accumulate DFI coefficient
grid%hcoeff_tot = grid%hcoeff_tot + hn
write(mess,*) 'grid%hcoeff_tot: ', grid%hcoeff_tot
call wrf_message(mess)
#endif
END SUBROUTINE dfi_accumulate
SUBROUTINE dfi_bck_init ( grid ) 1,15
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
USE module_utility
USE module_state_description
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER rc
CHARACTER*80 mess
INTERFACE
SUBROUTINE Setup_Timekeeping(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE Setup_Timekeeping
SUBROUTINE dfi_save_arrays(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_save_arrays
SUBROUTINE dfi_clear_accumulation(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_clear_accumulation
SUBROUTINE optfil_driver(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE optfil_driver
SUBROUTINE start_domain(grid, allowed_to_read)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
LOGICAL, INTENT(IN) :: allowed_to_read
END SUBROUTINE start_domain
#if (EM_CORE == 1)
SUBROUTINE rebalance_driver_dfi (grid)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
END SUBROUTINE rebalance_driver_dfi
#endif
END INTERFACE
grid%dfi_stage = DFI_BCK
! Negate time step
IF ( grid%time_step_dfi .gt. 0 ) THEN
CALL nl_set_time_step ( 1, -grid%time_step_dfi )
ELSE
CALL nl_set_time_step ( 1, -grid%time_step )
CALL nl_set_time_step_fract_num ( 1, -grid%time_step_fract_num )
ENDIF
CALL Setup_Timekeeping (grid)
!tgs 7apr11 - need to call start_domain here to reset bc initialization for negative dt
CALL start_domain ( grid , .TRUE. )
!tgs 7apr11 - save arrays should be done after start_domain to get correct grid%p field
CALL dfi_save_arrays ( grid )
! set physics options to zero
CALL nl_set_mp_physics( grid%id, 0 )
CALL nl_set_ra_lw_physics( grid%id, 0 )
CALL nl_set_ra_sw_physics( grid%id, 0 )
CALL nl_set_sf_surface_physics( grid%id, 0 )
CALL nl_set_sf_sfclay_physics( grid%id, 0 )
CALL nl_set_sf_urban_physics( grid%id, 0 )
CALL nl_set_bl_pbl_physics( grid%id, 0 )
CALL nl_set_cu_physics( grid%id, 0 )
CALL nl_set_damp_opt( grid%id, 0 )
CALL nl_set_sst_update( grid%id, 0 )
CALL nl_set_fractional_seaice( grid%id, 0 )
CALL nl_set_gwd_opt( grid%id, 0 )
CALL nl_set_feedback( grid%id, 0 )
! set bc
#if (EM_CORE == 1)
CALL nl_set_diff_6th_opt( grid%id, 0 )
CALL nl_set_constant_bc(1, grid%constant_bc)
CALL nl_set_use_adaptive_time_step( grid%id, .false. )
#endif
#ifdef WRF_CHEM
! set chemistry option to zero
CALL nl_set_chem_opt (grid%id, 0)
CALL nl_set_aer_ra_feedback (grid%id, 0)
CALL nl_set_io_form_auxinput5 (grid%id, 0)
CALL nl_set_io_form_auxinput7 (grid%id, 0)
CALL nl_set_io_form_auxinput8 (grid%id, 0)
#endif
! set diffusion to zero for backward integration
#if (EM_CORE == 1)
CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
IF ( grid%moist_adv_opt == 2 ) THEN
CALL nl_set_moist_adv_opt( grid%id, 0)
ENDIF
#endif
! If a request to do pressure level diags, then shut it off
! until we get to the real forecast part of this integration.
#if (EM_CORE == 1)
CALL nl_set_p_lev_diags( grid%id, grid%p_lev_diags_dfi)
#endif
#if (NMM_CORE == 1)
!
! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
!
CALL nl_get_time_step( grid%id, grid%time_step )
if (grid%time_step .lt. 0) then
! DT =-DT
write(mess,*) 'changing signs for backward integration'
call wrf_message(mess)
grid%CPGFV = -grid%CPGFV
grid%EN = -grid%EN
grid%ENT = -grid%ENT
grid%F4D = -grid%F4D
grid%F4Q = -grid%F4Q
grid%EF4T = -grid%EF4T
grid%EM(:) = -grid%EM(:)
grid%EMT(:) = -grid%EMT(:)
grid%F4Q2(:) = -grid%F4Q2(:)
grid%WPDAR(:,:)= -grid%WPDAR(:,:)
grid%CPGFU(:,:)= -grid%CPGFU(:,:)
grid%CURV(:,:)= -grid%CURV(:,:)
grid%FCP(:,:)= -grid%FCP(:,:)
grid%FAD(:,:)= -grid%FAD(:,:)
grid%F(:,:)= -grid%F(:,:)
endif
#endif
grid%start_subtime = domain_get_start_time ( grid )
grid%stop_subtime = domain_get_stop_time ( grid )
CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
!tgs 7apr11 - this call has been moved up
! CALL dfi_save_arrays ( grid )
CALL dfi_clear_accumulation( grid )
CALL optfil_driver(grid)
!tgs need to call start_domain here to reset bc initialization for negative dt
CALL start_domain ( grid , .TRUE. )
!tgs need to call rebalance here to remove imbalances in initial fields
! when config_flags%cycling=.true.
#if (EM_CORE == 1)
if(grid%cycling) then
! print *,' Rebalancing is on '
CALL rebalance_driver_dfi ( grid )
endif
#endif
END SUBROUTINE dfi_bck_init
SUBROUTINE dfi_fwd_init ( grid ) 1,13
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
USE module_utility, ONLY : WRFU_ClockSet
USE module_state_description
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER rc
CHARACTER*80 mess
INTERFACE
SUBROUTINE Setup_Timekeeping(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE Setup_Timekeeping
SUBROUTINE dfi_save_arrays(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_save_arrays
SUBROUTINE dfi_clear_accumulation(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_clear_accumulation
SUBROUTINE optfil_driver(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE optfil_driver
#if (EM_CORE == 1)
SUBROUTINE rebalance_driver_dfi (grid)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
END SUBROUTINE rebalance_driver_dfi
#endif
SUBROUTINE start_domain(grid, allowed_to_read)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
LOGICAL, INTENT(IN) :: allowed_to_read
END SUBROUTINE start_domain
END INTERFACE
grid%dfi_stage = DFI_FWD
! for Setup_Timekeeping to use when setting the clock
IF ( grid%time_step_dfi .gt. 0 ) THEN
CALL nl_set_time_step ( grid%id, grid%time_step_dfi )
ELSE
CALL nl_get_time_step( grid%id, grid%time_step )
CALL nl_get_time_step_fract_num( grid%id, grid%time_step_fract_num )
grid%time_step = abs(grid%time_step)
CALL nl_set_time_step( grid%id, grid%time_step )
grid%time_step_fract_num = abs(grid%time_step_fract_num)
CALL nl_set_time_step_fract_num( grid%id, grid%time_step_fract_num )
ENDIF
grid%itimestep=0
grid%xtime=0.
! reset physics options to normal
CALL nl_set_mp_physics( grid%id, grid%mp_physics)
CALL nl_set_ra_lw_physics( grid%id, grid%ra_lw_physics)
CALL nl_set_ra_sw_physics( grid%id, grid%ra_sw_physics)
CALL nl_set_sf_surface_physics( grid%id, grid%sf_surface_physics)
CALL nl_set_sf_sfclay_physics( grid%id, grid%sf_sfclay_physics)
CALL nl_set_sf_urban_physics( grid%id, grid%sf_urban_physics)
CALL nl_set_bl_pbl_physics( grid%id, grid%bl_pbl_physics)
CALL nl_set_cu_physics( grid%id, grid%cu_physics)
CALL nl_set_damp_opt( grid%id, grid%damp_opt)
CALL nl_set_sst_update( grid%id, 0)
CALL nl_set_fractional_seaice( grid%id, grid%fractional_seaice)
CALL nl_set_gwd_opt( grid%id, grid%gwd_opt)
CALL nl_set_feedback( grid%id, 0 )
#if (EM_CORE == 1)
CALL nl_set_diff_6th_opt( grid%id, grid%diff_6th_opt)
CALL nl_set_use_adaptive_time_step( grid%id, .false. )
! set bc
CALL nl_set_constant_bc(grid%id, grid%constant_bc)
#endif
#if (NMM_CORE == 1)
!
! CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
!
!mptest CALL nl_get_time_step( grid%id, grid%time_step )
!!! here need to key off something else being the "wrong" sign
if (grid%cpgfv .gt. 0) then
write(mess,*) 'changing signs for forward integration'
call wrf_message(mess)
grid%CPGFV = -grid%CPGFV
grid%EN = -grid%EN
grid%ENT = -grid%ENT
grid%F4D = -grid%F4D
grid%F4Q = -grid%F4Q
grid%EF4T = -grid%EF4T
grid%EM(:) = -grid%EM(:)
grid%EMT(:) = -grid%EMT(:)
grid%F4Q2(:) = -grid%F4Q2(:)
grid%WPDAR(:,:)= -grid%WPDAR(:,:)
grid%CPGFU(:,:)= -grid%CPGFU(:,:)
grid%CURV(:,:)= -grid%CURV(:,:)
grid%FCP(:,:)= -grid%FCP(:,:)
grid%FAD(:,:)= -grid%FAD(:,:)
grid%F(:,:)= -grid%F(:,:)
endif
#endif
!#ifdef WRF_CHEM
! ! reset chem option to normal
! CALL nl_set_chem_opt( grid%id, grid%chem_opt)
! CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
!#endif
#if (EM_CORE == 1)
! reset km_opt to normal
CALL nl_set_km_opt( grid%id, grid%km_opt)
CALL nl_set_moist_adv_opt( grid%id, grid%moist_adv_opt)
#endif
CALL Setup_Timekeeping (grid)
grid%start_subtime = domain_get_start_time ( head_grid )
grid%stop_subtime = domain_get_stop_time ( head_grid )
CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
IF ( grid%dfi_opt .EQ. DFI_DFL ) THEN
CALL dfi_save_arrays ( grid )
END IF
CALL dfi_clear_accumulation( grid )
CALL optfil_driver(grid)
!tgs need to call it here to reset bc initialization for positive time_step
CALL start_domain ( grid , .TRUE. )
END SUBROUTINE dfi_fwd_init
SUBROUTINE dfi_fst_init ( grid ) 1,9
USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time
USE module_state_description
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
CHARACTER (LEN=80) :: wrf_error_message
INTERFACE
SUBROUTINE Setup_Timekeeping(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE Setup_Timekeeping
SUBROUTINE dfi_save_arrays(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_save_arrays
SUBROUTINE dfi_clear_accumulation(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_clear_accumulation
SUBROUTINE optfil_driver(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE optfil_driver
#if (EM_CORE == 1)
SUBROUTINE rebalance_driver_dfi (grid)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
END SUBROUTINE rebalance_driver_dfi
#endif
SUBROUTINE start_domain(grid, allowed_to_read)
USE module_domain, ONLY : domain
TYPE (domain) :: grid
LOGICAL, INTENT(IN) :: allowed_to_read
END SUBROUTINE start_domain
END INTERFACE
grid%dfi_stage = DFI_FST
! reset time_step to normal and adaptive_time_step
CALL nl_set_time_step( grid%id, grid%time_step )
grid%itimestep=0
grid%xtime=0. ! BUG: This will probably not work for all DFI options
! only use adaptive time stepping for forecast
#if (EM_CORE == 1)
if (grid%id == 1) then
CALL nl_set_use_adaptive_time_step( 1, grid%use_adaptive_time_step )
endif
CALL nl_set_sst_update( grid%id, grid%sst_update)
! reset to normal bc
CALL nl_set_constant_bc(grid%id, .false.)
#endif
CALL nl_set_feedback( grid%id, grid%feedback )
#ifdef WRF_CHEM
! reset chem option to normal
CALL nl_set_chem_opt( grid%id, grid%chem_opt)
CALL nl_set_aer_ra_feedback (grid%id, grid%aer_ra_feedback)
CALL nl_set_io_form_auxinput5 (grid%id, grid%io_form_auxinput5)
CALL nl_set_io_form_auxinput7 (grid%id, grid%io_form_auxinput7)
CALL nl_set_io_form_auxinput8 (grid%id, grid%io_form_auxinput8)
#endif
#if (EM_CORE == 1)
! reset pressure level diags to normal
CALL nl_set_p_lev_diags ( grid%id, grid%p_lev_diags)
#endif
CALL Setup_Timekeeping (grid)
grid%start_subtime = domain_get_start_time ( grid )
grid%stop_subtime = domain_get_stop_time ( grid )
CALL start_domain ( grid , .TRUE. )
END SUBROUTINE dfi_fst_init
SUBROUTINE dfi_write_initialized_state( grid ) 1,8
! Driver layer
USE module_domain, ONLY : domain, head_grid
USE module_io_domain
USE module_configure, ONLY : grid_config_rec_type, model_config_rec
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER :: fid, ierr
CHARACTER (LEN=80) :: wrf_error_message
CHARACTER (LEN=80) :: rstname
CHARACTER (LEN=132) :: message
TYPE (grid_config_rec_type) :: config_flags
CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
WRITE (wrf_err_message,'(A,I4)') 'Writing out initialized model state'
CALL wrf_message(TRIM(wrf_err_message))
WRITE (rstname,'(A,I2.2)')'wrfinput_initialized_d',grid%id
CALL open_w_dataset ( fid, TRIM(rstname), grid, config_flags, output_input, "DATASET=INPUT", ierr )
IF ( ierr .NE. 0 ) THEN
WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
CALL WRF_ERROR_FATAL ( message )
END IF
CALL output_input ( fid, grid, config_flags, ierr )
CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
END SUBROUTINE dfi_write_initialized_state
SUBROUTINE tdfi_write_analyzed_state ( grid ) 1,7
! Driver layer
USE module_domain, ONLY : domain, head_grid
USE module_io_domain
USE module_configure, ONLY : grid_config_rec_type, model_config_rec
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER :: fid, ierr
CHARACTER (LEN=80) :: wrf_error_message
CHARACTER (LEN=80) :: rstname
CHARACTER (LEN=132) :: message
TYPE (grid_config_rec_type) :: config_flags
CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
rstname = 'wrfout_analyzed_d01'
CALL open_w_dataset ( fid, TRIM(rstname),grid, config_flags, output_input, "DATASET=INPUT", ierr )
IF ( ierr .NE. 0 ) THEN
WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
CALL WRF_ERROR_FATAL ( message )
END IF
CALL output_input ( fid, grid, config_flags, ierr )
CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
END SUBROUTINE tdfi_write_analyzed_state
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI array reset group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_array_reset ( ) 4,4
USE module_domain, ONLY : domain, head_grid, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_array_reset_recurse
END INTERFACE
! Copy filtered arrays back into state arrays in grid structure, and
! restore original land surface fields
CALL dfi_array_reset_recurse(head_grid)
CALL set_current_grid_ptr( head_grid )
END SUBROUTINE wrf_dfi_array_reset
SUBROUTINE dfi_array_reset( grid ) 1,6
USE module_domain, ONLY : domain
! USE module_configure
! USE module_driver_constants
! USE module_machine
! USE module_dm
USE module_model_constants
USE module_state_description
IMPLICIT NONE
INTEGER :: its, ite, jts, jte, kts, kte, &
i, j, k
! Input data.
TYPE(domain) , POINTER :: grid
! local
! real p1000mb,eps,svp1,svp2,svp3,svpt0
real eps
! parameter (p1000mb = 1.e+05, eps=0.622,svp1=0.6112,svp3=29.65,svpt0=273.15)
parameter (eps=0.622)
integer nsoil , ncount
REAL es,qs,pol,tx,temp,pres,rslf,rsif
CHARACTER*80 mess
ncount = 0
nsoil = grid%num_soil_layers
IF ( grid%dfi_opt .EQ. DFI_NODFI ) RETURN
! Set dynamic variables
! divide by total DFI coefficient
#if (EM_CORE == 1)
grid%mu_2(:,:) = grid%dfi_mu(:,:) / grid%hcoeff_tot
grid%u_2(:,:,:) = grid%dfi_u(:,:,:) / grid%hcoeff_tot
grid%v_2(:,:,:) = grid%dfi_v(:,:,:) / grid%hcoeff_tot
grid%w_2(:,:,:) = grid%dfi_w(:,:,:) / grid%hcoeff_tot
grid%ww(:,:,:) = grid%dfi_ww(:,:,:) / grid%hcoeff_tot
grid%t_2(:,:,:) = grid%dfi_t(:,:,:) / grid%hcoeff_tot
grid%phb(:,:,:) = grid%dfi_phb(:,:,:) / grid%hcoeff_tot
grid%ph0(:,:,:) = grid%dfi_ph0(:,:,:) / grid%hcoeff_tot
grid%php(:,:,:) = grid%dfi_php(:,:,:) / grid%hcoeff_tot
grid%p(:,:,:) = grid%dfi_p(:,:,:) / grid%hcoeff_tot
grid%ph_2(:,:,:) = grid%dfi_ph(:,:,:) / grid%hcoeff_tot
grid%tke_2(:,:,:) = grid%dfi_tke(:,:,:) / grid%hcoeff_tot
grid%al(:,:,:) = grid%dfi_al(:,:,:) / grid%hcoeff_tot
grid%alt(:,:,:) = grid%dfi_alt(:,:,:) / grid%hcoeff_tot
grid%pb(:,:,:) = grid%dfi_pb(:,:,:) / grid%hcoeff_tot
IF ( grid%dfi_savehydmeteors .EQ. 0 ) then
! print *,'Normal DFI'
grid%moist(:,:,:,:) = max(0.,grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot)
grid%scalar(:,:,:,:) = max(0.,grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot)
ELSE
! print *,'In dfi_array_reset, QV comp, dfi_save_hydrometeors=1'
grid%moist(:,:,:,P_QV) = max(0.,grid%dfi_moist(:,:,:,P_QV) / grid%hcoeff_tot)
ENDIF
if ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
! grid%qvg(:,:) = grid%dfi_qvg(:,:) / grid%hcoeff_tot
endif
IF ( grid%dfi_savehydmeteors .EQ. 1 ) then
! print *,'restore initial surface fields'
! restore initial fields
grid%SNOW (:,:) = grid%dfi_SNOW (:,:)
grid%SNOWH (:,:) = grid%dfi_SNOWH (:,:)
grid%SNOWC (:,:) = grid%dfi_SNOWC (:,:)
grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
grid%TSK (:,:) = grid%dfi_TSK (:,:)
grid%TSLB (:,:,:) = grid%dfi_TSLB (:,:,:)
grid%SMOIS (:,:,:) = grid%dfi_SMOIS (:,:,:)
IF ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
grid%QVG (:,:) = grid%dfi_QVG (:,:)
grid%TSNAV (:,:) = grid%dfi_TSNAV (:,:)
grid%SOILT1(:,:) = grid%dfi_SOILT1(:,:)
grid%SMFR3D(:,:,:) = grid%dfi_SMFR3D (:,:,:)
grid%KEEPFR3DFLAG(:,:,:) = grid%dfi_KEEPFR3DFLAG(:,:,:)
ENDIF
ELSE
! print *, 'take into account DFI snowfall'
! restore initial fields
its = grid%sp31 ; ite = grid%ep31 ;
jts = grid%sp33 ; jte = grid%ep33 ;
DO j=jts,jte
DO i=its,ite
grid%SNOW (i,j) = max(grid%SNOW (i,j),grid%dfi_SNOW (i,j))
grid%SNOWH (i,j) = max(grid%SNOWH (i,j),grid%dfi_SNOWH (i,j))
grid%SNOWC (i,j) = max(grid%SNOWC (i,j),grid%dfi_SNOWC (i,j))
grid%CANWAT(i,j) = max(grid%CANWAT(i,j),grid%dfi_CANWAT(i,j))
IF(grid%SNOWC (i,j) .eq.0.) then
grid%TSK (i,j) = grid%dfi_TSK (i,j)
do k=1,nsoil
grid%TSLB (i,k,j) = grid%dfi_TSLB (i,k,j)
grid%SMOIS (i,k,j) = grid%dfi_SMOIS (i,k,j)
IF ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
grid%KEEPFR3DFLAG(i,k,j) = grid%dfi_KEEPFR3DFLAG(i,k,j)
grid%SMFR3D(i,k,j) = grid%dfi_SMFR3D (i,k,j)
ENDIF
enddo
IF ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
grid%QVG (i,j) = grid%dfi_QVG (i,j)
grid%TSNAV (i,j) = grid%dfi_TSNAV (i,j)
grid%SOILT1(i,j) = grid%dfi_SOILT1(i,j)
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
! restore analized hydrometeor and scalar fileds
IF ( grid%dfi_savehydmeteors .EQ. 1 ) then
! print *,'In dfi_array_reset - restore initial hydrometeors'
! grid%moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) !tgs
grid%moist(:,:,:,P_QC) = grid%dfi_moist(:,:,:,P_QC) !tgs
grid%moist(:,:,:,P_QR) = grid%dfi_moist(:,:,:,P_QR) !tgs
grid%moist(:,:,:,P_QI) = grid%dfi_moist(:,:,:,P_QI) !tgs
grid%moist(:,:,:,P_QS) = grid%dfi_moist(:,:,:,P_QS) !tgs
grid%moist(:,:,:,P_QG) = grid%dfi_moist(:,:,:,P_QG) !tgs
grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:)
if(grid%dfi_stage .EQ. DFI_FWD) then
!tgs change QV to restore initial RH field after the diabatic DFI
its = grid%sp31 ; ite = grid%ep31 ;
kts = grid%sp32 ; kte = grid%ep32 ;
jts = grid%sp33 ; jte = grid%ep33 ;
DO j=jts,jte
DO i=its,ite
do k = kts , kte
temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
** (r_d / Cp)
pres = grid%p(i,k,j)+grid%pb(i,k,j)
!tgs rslf - function to compute qs from Thompson microphysics
qs = rslf(pres, temp)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-thomp',temp,pres,qs
! endif
!! es=svp1*10.*EXP(svp2*(temp-svpt0)/(temp-svp3))
!! qs=0.622*es/(pres/100.-es)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-wrf',temp,pres,qs
! endif
Tx = temp-273.15
POL = 0.99999683 + TX*(-0.90826951E-02 + &
TX*(0.78736169E-04 + TX*(-0.61117958E-06 + &
TX*(0.43884187E-08 + TX*(-0.29883885E-10 + &
TX*(0.21874425E-12 + TX*(-0.17892321E-14 + &
TX*(0.11112018E-16 + TX*(-0.30994571E-19)))))))))
es = 6.1078/POL**8
!! qs = 0.62197*es / (pres/100. - es*0.37803)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-ruc',temp,pres,qs
! endif
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'before grid%moist (i,k,j,P_QV),grid%dfi_rh(i,k,j)', &
! grid%moist(i,k,j,P_QV),grid%dfi_rh(i,k,j)
! endif
! Ensure that saturated points going into DFI remain saturated after DFI
IF(grid%dfi_rh(i,k,j) == 1. ) then
grid%moist (i,k,j,P_QV) = MAX(grid%moist (i,k,j,P_QV),grid%dfi_rh(i,k,j)*QS)
ncount=ncount+1
ENDIF
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'after grid%moist (i,k,j,P_QV)', &
! grid%moist(i,k,j,P_QV)
! endif
enddo
ENDDO
ENDDO
! print *,'QV reset to saturation at ncount= ',ncount
endif
!tgs need to call rebalance here again after hydrometeor and QV reset
! print *,' Rebalancing is on '
CALL rebalance_driver_dfi ( grid )
ENDIF
#endif
#if (NMM_CORE == 1)
write(mess,*) ' divide by grid%hcoeff_tot: ', grid%hcoeff_tot
call wrf_message(mess)
if (grid%hcoeff_tot .lt. 0.0001) then
call wrf_error_fatal("bad grid%hcoeff_tot")
endif
grid%pd(:,:) = grid%dfi_pd(:,:) / grid%hcoeff_tot
grid%pint(:,:,:) = grid%dfi_pint(:,:,:) / grid%hcoeff_tot
! grid%dwdt(:,:,:) = grid%dfi_dwdt(:,:,:) / grid%hcoeff_tot
grid%t(:,:,:) = grid%dfi_t(:,:,:) / grid%hcoeff_tot
grid%q(:,:,:) = grid%dfi_q(:,:,:) / grid%hcoeff_tot
grid%q2(:,:,:) = grid%dfi_q2(:,:,:) / grid%hcoeff_tot
grid%cwm(:,:,:) = grid%dfi_cwm(:,:,:) / grid%hcoeff_tot
grid%u(:,:,:) = grid%dfi_u(:,:,:) / grid%hcoeff_tot
grid%v(:,:,:) = grid%dfi_v(:,:,:) / grid%hcoeff_tot
grid%moist(:,:,:,:) = grid%dfi_moist(:,:,:,:) / grid%hcoeff_tot
grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot
! restore initial fields
grid%SNOW(:,:) = grid%dfi_SNOW(:,:)
grid%SNOWH(:,:) = grid%dfi_SNOWH(:,:)
! grid%SNOWC(:,:) = grid%dfi_SNOWC(:,:)
grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
grid%NMM_TSK(:,:) = grid%dfi_NMM_TSK(:,:)
! save soil fields
grid%STC(:,:,:) = grid%dfi_STC(:,:,:)
grid%SMC(:,:,:) = grid%dfi_SMC(:,:,:)
grid%SH2O(:,:,:) = grid%dfi_SH2O(:,:,:)
#endif
END SUBROUTINE dfi_array_reset
SUBROUTINE optfil_driver( grid ) 2,4
USE module_domain, ONLY : domain
USE module_utility
! USE module_wrf_error
! USE module_timing
! USE module_date_time
! USE module_configure
USE module_state_description
USE module_model_constants
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
! Local variables
integer :: nstep2, nstepmax, rundfi, i, rc,hr,min,sec
integer :: yr,jday
real :: timestep, tauc
TYPE(WRFU_TimeInterval) :: run_interval
CHARACTER*80 mess
timestep=abs(grid%dt)
run_interval = grid%stop_subtime - grid%start_subtime
CALL WRFU_TimeGet(grid%start_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
CALL WRFU_TimeGet(grid%stop_subtime, YY=yr, dayofYear=jday, H=hr, M=min, S=sec, rc=rc)
CALL WRFU_TimeIntervalGet( run_interval, S=rundfi, rc=rc )
rundfi = abs(rundfi)
nstep2= ceiling((1.0 + real(rundfi)/timestep) / 2.0)
! nstep2 is equal to a half of timesteps per initialization period,
! should not exceed nstepmax
tauc = real(grid%dfi_cutoff_seconds)
! Get DFI coefficient
grid%hcoeff(:) = 0.0
IF ( grid%dfi_nfilter < 0 .OR. grid%dfi_nfilter > 8 ) THEN
write(mess,*) 'Invalid filter specified in namelist.'
call wrf_message(mess)
ELSE
call dfcoef(nstep2-1, grid%dt, tauc, grid%dfi_nfilter, grid%hcoeff)
END IF
IF ( MOD(int(1.0 + real(rundfi)/timestep),2) /= 0 ) THEN
DO i=1,nstep2-1
grid%hcoeff(2*nstep2-i) = grid%hcoeff(i)
END DO
ELSE
DO i=1,nstep2
grid%hcoeff(2*nstep2-i+1) = grid%hcoeff(i)
END DO
END IF
END SUBROUTINE optfil_driver
SUBROUTINE dfi_clear_accumulation( grid ) 2,1
USE module_domain, ONLY : domain
! USE module_configure
! USE module_driver_constants
! USE module_machine
! USE module_dm
! USE module_model_constants
USE module_state_description
IMPLICIT NONE
! Input data.
TYPE(domain) , POINTER :: grid
#if (EM_CORE == 1)
grid%dfi_mu(:,:) = 0.
grid%dfi_u(:,:,:) = 0.
grid%dfi_v(:,:,:) = 0.
grid%dfi_w(:,:,:) = 0.
grid%dfi_ww(:,:,:) = 0.
grid%dfi_t(:,:,:) = 0.
grid%dfi_phb(:,:,:) = 0.
grid%dfi_ph0(:,:,:) = 0.
grid%dfi_php(:,:,:) = 0.
grid%dfi_p(:,:,:) = 0.
grid%dfi_ph(:,:,:) = 0.
grid%dfi_tke(:,:,:) = 0.
grid%dfi_al(:,:,:) = 0.
grid%dfi_alt(:,:,:) = 0.
grid%dfi_pb(:,:,:) = 0.
IF ( grid%dfi_savehydmeteors .EQ. 0 ) then
! print *,'normal DFI'
grid%dfi_moist(:,:,:,:) = 0.
grid%dfi_scalar(:,:,:,:) = 0.
ELSE
! print *,'In dfi_clear_accumulation, clear dfi_QV - dfi_savehydmeteors=1'
grid%dfi_moist(:,:,:,P_QV) = 0.
ENDIF
if ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
! grid%dfi_qvg(:,:) = 0.
endif
#endif
#if (NMM_CORE == 1)
grid%dfi_pd(:,:) = 0.
grid%dfi_pint(:,:,:) = 0.
grid%dfi_dwdt(:,:,:) = 0.
grid%dfi_t(:,:,:) = 0.
grid%dfi_q(:,:,:) = 0.
grid%dfi_q2(:,:,:) = 0.
grid%dfi_cwm(:,:,:) = 0.
grid%dfi_u(:,:,:) = 0.
grid%dfi_v(:,:,:) = 0.
grid%dfi_moist(:,:,:,:) = 0.
grid%dfi_scalar(:,:,:,:) = 0.
#endif
grid%hcoeff_tot = 0.0
END SUBROUTINE dfi_clear_accumulation
SUBROUTINE dfi_save_arrays( grid ) 2,3
USE module_domain, ONLY : domain
! USE module_configure
! USE module_driver_constants
! USE module_machine
! USE module_dm
USE module_model_constants
USE module_state_description
IMPLICIT NONE
INTEGER :: its, ite, jts, jte, kts, kte, &
i, j, k
! Input data.
TYPE(domain) , POINTER :: grid
! local
REAL es,qs,pol,tx,temp,pres,rslf
#if (EM_CORE == 1)
! save surface 2-D fields
grid%dfi_SNOW(:,:) = grid%SNOW(:,:)
grid%dfi_SNOWH(:,:) = grid%SNOWH(:,:)
grid%dfi_SNOWC(:,:) = grid%SNOWC(:,:)
grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
grid%dfi_TSK(:,:) = grid%TSK(:,:)
! save soil fields
grid%dfi_TSLB(:,:,:) = grid%TSLB(:,:,:)
grid%dfi_SMOIS(:,:,:) = grid%SMOIS(:,:,:)
! RUC LSM only, need conditional
IF ( grid%sf_surface_physics .EQ. RUCLSMSCHEME ) then
grid%dfi_QVG(:,:) = grid%QVG(:,:)
grid%dfi_SOILT1(:,:) = grid%SOILT1(:,:)
! use this array to save initial mu_2
! grid%dfi_TSNAV(:,:) = grid%mu_2(:,:)
grid%dfi_TSNAV(:,:) = grid%TSNAV(:,:)
grid%dfi_SMFR3D(:,:,:) = grid%SMFR3D(:,:,:)
grid%dfi_KEEPFR3DFLAG(:,:,:) = grid%KEEPFR3DFLAG(:,:,:)
ENDIF
#endif
#if (NMM_CORE == 1)
! save surface 2-D fields
grid%dfi_SNOW(:,:) = grid%SNOW(:,:)
grid%dfi_SNOWH(:,:) = grid%SNOWH(:,:)
! grid%dfi_SNOWC(:,:) = grid%SNOWC(:,:)
grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
grid%dfi_NMM_TSK(:,:) = grid%NMM_TSK(:,:)
! save soil fields
grid%dfi_STC(:,:,:) = grid%STC(:,:,:)
grid%dfi_SMC(:,:,:) = grid%SMC(:,:,:)
grid%dfi_SH2O(:,:,:) = grid%SH2O(:,:,:)
#endif
#if (EM_CORE == 1)
! save hydrometeor and scalar fields
IF ( grid%dfi_savehydmeteors .EQ. 1 ) then !tgs
! print *,'In dfi_save_arrays - save initial hydrometeors'
!! grid%dfi_moist(:,:,:,:) = grid%moist(:,:,:,:)
grid%dfi_moist(:,:,:,P_QC) = grid%moist(:,:,:,P_QC)
grid%dfi_moist(:,:,:,P_QR) = grid%moist(:,:,:,P_QR)
grid%dfi_moist(:,:,:,P_QI) = grid%moist(:,:,:,P_QI)
grid%dfi_moist(:,:,:,P_QS) = grid%moist(:,:,:,P_QS)
grid%dfi_moist(:,:,:,P_QG) = grid%moist(:,:,:,P_QG)
grid%dfi_scalar(:,:,:,:) = grid%scalar(:,:,:,:)
ENDIF
if(grid%dfi_stage .EQ. DFI_BCK) then
! compute initial RH field to be reintroduced after the diabatic DFI
its = grid%sp31 ; ite = grid%ep31 ;
kts = grid%sp32 ; kte = grid%ep32 ;
jts = grid%sp33 ; jte = grid%ep33 ;
DO j=jts,jte
DO i=its,ite
do k = kts , kte
temp = (grid%t_2(i,k,j)+t0)*( (grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb )&
** (r_d / Cp)
pres = grid%p(i,k,j)+grid%pb(i,k,j)
!tgs rslf - function to compute qs from Thompson microphysics
qs = rslf(pres, temp)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-thomp',temp,pres,qs
! endif
!! es=svp1*10.*EXP(svp2*(temp-svpt0)/(temp-svp3))
!! qs=0.622*es/(pres/100.-es)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-wrf',temp,pres,qs
! endif
Tx = temp-svpt0
POL = 0.99999683 + TX*(-0.90826951E-02 + &
TX*(0.78736169E-04 + TX*(-0.61117958E-06 + &
TX*(0.43884187E-08 + TX*(-0.29883885E-10 + &
TX*(0.21874425E-12 + TX*(-0.17892321E-14 + &
TX*(0.11112018E-16 + TX*(-0.30994571E-19)))))))))
es = 6.1078/POL**8
!! qs = 0.62197*es / (pres/100. - es*0.37803)
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs-ruc',temp,pres,qs
! endif
grid%dfi_rh (i,k,j) = MIN(1.,MAX(0.,grid%moist(i,k,j,P_QV)/qs))
! if(i.eq. 464 .and. j.eq. 212 .and. k.eq.6) then
! print *,'temp,pres,qs,grid%moist (i,k,j,P_QV),grid%dfi_rh (i,k,j)',temp,pres,qs, &
! grid%moist(i,k,j,P_QV),grid%dfi_rh (i,k,j)
! endif
!tgs saturation check for points with water or ice clouds
IF ((grid%moist (i,k,j,P_QC) .GT. 1.e-6 .or. &
grid%moist (i,k,j,P_QI) .GT. 1.e-6) .and. &
grid%dfi_rh (i,k,j) .lt. 1.) THEN
grid%dfi_rh (i,k,j)=1.
ENDIF
end do
END DO
ENDDO
endif
#endif
END SUBROUTINE dfi_save_arrays
SUBROUTINE dfcoef (NSTEPS,DT,TAUC,NFILT,H) 1,11
!
! calculate filter weights with selected window.
!
! peter lynch and xiang-yu huang
!
! ref: see hamming, r.w., 1989: digital filters,
! prentice-hall international. 3rd edition.
!
! input: nsteps - number of timesteps
! forward or backward.
! dt - time step in seconds.
! tauc - cut-off period in seconds.
! nfilt - indicator for selected filter.
!
! output: h - array(0:nsteps) with the
! required filter weights
!
!------------------------------------------------------------
implicit none
integer, intent(in) :: nsteps, nfilt
real , intent(in) :: dt, tauc
real, intent(out) :: h(1:nsteps+1)
! Local data
integer :: n
real :: pi, omegac, x, window, deltat
real :: hh(0:nsteps)
pi=4*ATAN(1.)
deltat=ABS(dt)
! windows are defined by a call and held in hh.
if ( nfilt .eq. -1) call debug (nsteps,h)
IF ( NFILT .EQ. 0 ) CALL UNIFORM (NSTEPS,HH)
IF ( NFILT .EQ. 1 ) CALL LANCZOS (NSTEPS,HH)
IF ( NFILT .EQ. 2 ) CALL HAMMING (NSTEPS,HH)
IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
IF ( NFILT .EQ. 4 ) CALL KAISER (NSTEPS,HH)
IF ( NFILT .EQ. 5 ) CALL POTTER2 (NSTEPS,HH)
IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
IF ( NFILT .LE. 6 ) THEN ! sinc-windowed filters
! calculate the cutoff frequency
OMEGAC = 2.*PI/TAUC
DO N=0,NSTEPS
WINDOW = HH(N)
IF ( N .EQ. 0 ) THEN
X = (OMEGAC*DELTAT/PI)
ELSE
X = SIN(N*OMEGAC*DELTAT)/(N*PI)
END IF
HH(N) = X*WINDOW
END DO
! normalize the sums to be unity
CALL NORMLZ(HH,NSTEPS)
DO N=0,NSTEPS
H(N+1) = HH(NSTEPS-N)
END DO
ELSE IF ( NFILT .EQ. 7 ) THEN ! dolph filter
CALL DOLPH(DT,TAUC,NSTEPS,H)
ELSE IF ( NFILT .EQ. 8 ) THEN ! 2nd order, 2nd type quick start filter (RHO)
CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS)
END IF
RETURN
END SUBROUTINE dfcoef
SUBROUTINE NORMLZ(HH,NMAX) 1
! normalize the sum of hh to be unity
implicit none
integer, intent(in) :: nmax
real , dimension(0:nmax), intent(out) :: hh
! local data
real :: sumhh
integer :: n
SUMHH = HH(0)
DO N=1,NMAX
SUMHH = SUMHH + 2*HH(N)
ENDDO
DO N=0,NMAX
HH(N) = HH(N)/SUMHH
ENDDO
RETURN
END subroutine normlz
subroutine debug(nsteps, ww) 1
implicit none
integer, intent(in) :: nsteps
real , dimension(0:nsteps), intent(out) :: ww
integer :: n
do n=0,nsteps
ww(n)=0
end do
ww(int(nsteps/2))=1
return
end subroutine debug
SUBROUTINE UNIFORM(NSTEPS,WW) 1
! define uniform or rectangular window function.
implicit none
integer, intent(in) :: nsteps
real , dimension(0:nsteps), intent(out) :: ww
integer :: n
DO N=0,NSTEPS
WW(N) = 1.
ENDDO
RETURN
END subroutine uniform
SUBROUTINE LANCZOS(NSTEPS,WW) 1
! define (genaralised) lanczos window function.
implicit none
integer, parameter :: nmax = 1000
integer, intent(in) :: nsteps
real , dimension(0:nmax), intent(out) :: ww
integer :: n
real :: power, pi, w
! (for the usual lanczos window, power = 1 )
POWER = 1
PI=4*ATAN(1.)
DO N=0,NSTEPS
IF ( N .EQ. 0 ) THEN
W = 1.0
ELSE
W = SIN(N*PI/(NSTEPS+1)) / ( N*PI/(NSTEPS+1))
ENDIF
WW(N) = W**POWER
ENDDO
RETURN
END SUBROUTINE lanczos
SUBROUTINE HAMMING(NSTEPS,WW) 1
! define (genaralised) hamming window function.
implicit none
integer, intent(in) :: nsteps
real, dimension(0:nsteps) :: ww
integer :: n
real :: alpha, pi, w
! (for the usual hamming window, alpha=0.54,
! for the hann window, alpha=0.50).
ALPHA=0.54
PI=4*ATAN(1.)
DO N=0,NSTEPS
IF ( N .EQ. 0 ) THEN
W = 1.0
ELSE
W = ALPHA + (1-ALPHA)*COS(N*PI/(NSTEPS))
ENDIF
WW(N) = W
ENDDO
RETURN
END SUBROUTINE hamming
SUBROUTINE BLACKMAN(NSTEPS,WW) 1
! define blackman window function.
implicit none
integer, intent(in) :: nsteps
real, dimension(0:nsteps) :: ww
integer :: n
real :: pi, w
PI=4*ATAN(1.)
DO N=0,NSTEPS
IF ( N .EQ. 0 ) THEN
W = 1.0
ELSE
W = 0.42 + 0.50*COS( N*PI/(NSTEPS)) &
+ 0.08*COS(2*N*PI/(NSTEPS))
ENDIF
WW(N) = W
ENDDO
RETURN
END SUBROUTINE blackman
SUBROUTINE KAISER(NSTEPS,WW) 1,2
! define kaiser window function.
implicit none
real, external :: bessi0
integer, intent(in) :: nsteps
real, dimension(0:nsteps) :: ww
integer :: n
real :: alpha, xi0a, xn, as
ALPHA = 1
XI0A = BESSI0(ALPHA)
DO N=0,NSTEPS
XN = N
AS = ALPHA*SQRT(1.-(XN/NSTEPS)**2)
WW(N) = BESSI0(AS) / XI0A
ENDDO
RETURN
END SUBROUTINE kaiser
REAL FUNCTION BESSI0(X) 2
! from numerical recipes (press, et al.)
implicit none
real(8) :: Y
real(8) :: P1 = 1.0d0
real(8) :: P2 = 3.5156230D0
real(8) :: P3 = 3.0899424D0
real(8) :: P4 = 1.2067492D0
real(8) :: P5 = 0.2659732D0
real(8) :: P6 = 0.360768D-1
real(8) :: P7 = 0.45813D-2
real*8 :: Q1 = 0.39894228D0
real*8 :: Q2 = 0.1328592D-1
real*8 :: Q3 = 0.225319D-2
real*8 :: Q4 = -0.157565D-2
real*8 :: Q5 = 0.916281D-2
real*8 :: Q6 = -0.2057706D-1
real*8 :: Q7 = 0.2635537D-1
real*8 :: Q8 = -0.1647633D-1
real*8 :: Q9 = 0.392377D-2
real :: x, ax
IF (ABS(X).LT.3.75) THEN
Y=(X/3.75)**2
BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
ELSE
AX=ABS(X)
Y=3.75/AX
BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 &
+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
ENDIF
RETURN
END FUNCTION bessi0
SUBROUTINE POTTER2(NSTEPS,WW) 1,1
! define potter window function.
! modified to fall off over twice the range.
implicit none
integer, intent(in) :: nsteps
real, dimension(0:nsteps),intent(out) :: ww
integer :: n
real :: ck, sum, arg
! local data
real, dimension(0:3) :: d
real :: pi
integer :: ip
d(0) = 0.35577019
d(1) = 0.2436983
d(2) = 0.07211497
d(3) = 0.00630165
PI=4*ATAN(1.)
CK = 1.0
DO N=0,NSTEPS
IF (N.EQ.NSTEPS) CK = 0.5
ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
!min--- modification in next statement
ARG = ARG/2.
!min--- end of modification
SUM = D(0)
DO IP=1,3
SUM = SUM + 2.*D(IP)*COS(ARG*FLOAT(IP))
END DO
WW(N) = CK*SUM
END DO
RETURN
END SUBROUTINE potter2
SUBROUTINE dolphwin(m, window) 1,1
! calculation of dolph-chebyshev window or, for short,
! dolph window, using the expression in the reference:
!
! antoniou, andreas, 1993: digital filters: analysis,
! design and applications. mcgraw-hill, inc., 689pp.
!
! the dolph window is optimal in the following sense:
! for a given main-lobe width, the stop-band attenuation
! is minimal; for a given stop-band level, the main-lobe
! width is minimal.
!
! it is possible to specify either the ripple-ratio r
! or the stop-band edge thetas.
IMPLICIT NONE
! Arguments
INTEGER, INTENT(IN) :: m
REAL, DIMENSION(0:M), INTENT(OUT) :: window
! local data
REAL, DIMENSION(0:2*M) :: t
REAL, DIMENSION(0:M) :: w, time
REAL :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
INTEGER :: n, nm1, nt, i
PI = 4*ATAN(1.D0)
THETAS = 2*PI/M
N = 2*M+1
NM1 = N-1
X0 = 1/COS(THETAS/2)
TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
RR = 0.5*(TERM1+TERM2)
R = 1/RR
DB = 20*LOG10(R)
WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
DO NT=0,M
SUM = RR
DO I=1,M
ARG = X0*COS(I*PI/N)
CALL CHEBY(T,NM1,ARG)
TERM1 = T(NM1)
TERM2 = COS(2*NT*PI*I/N)
SUM = SUM + 2*TERM1*TERM2
ENDDO
W(NT) = SUM/N
TIME(NT) = NT
ENDDO
! fill up the array for return
DO NT=0,M
WINDOW(NT) = W(NT)
ENDDO
RETURN
END SUBROUTINE dolphwin
SUBROUTINE dolph(deltat, taus, m, window) 1,7
! calculation of dolph-chebyshev window or, for short,
! dolph window, using the expression in the reference:
!
! antoniou, andreas, 1993: digital filters: analysis,
! design and applications. mcgraw-hill, inc., 689pp.
!
! the dolph window is optimal in the following sense:
! for a given main-lobe width, the stop-band attenuation
! is minimal; for a given stop-band level, the main-lobe
! width is minimal.
IMPLICIT NONE
! Arguments
INTEGER, INTENT(IN) :: m
REAL, DIMENSION(0:M), INTENT(OUT) :: window
REAL, INTENT(IN) :: deltat, taus
! local data
integer, PARAMETER :: NMAX = 5000
REAL, dimension(0:NMAX) :: t, w, time
real, dimension(0:2*nmax) :: w2
INTEGER :: NPRPE=0 ! no of pe
CHARACTER*80 :: MES
real :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
integer :: n, nm1, i, nt
PI = 4*ATAN(1.D0)
WRITE (mes,'(A,F8.2,A,F10.2)') 'In dolph, deltat = ',deltat,' taus = ',taus
CALL wrf_message(TRIM(mes))
N = 2*M+1
NM1 = N-1
THETAS = 2*PI*ABS(DELTAT/TAUS)
X0 = 1/COS(THETAS/2)
TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
RR = 0.5*(TERM1+TERM2)
R = 1/RR
DB = 20*LOG10(R)
WRITE (mes,'(A,2I8)') 'In dolph: M,N = ', M,N
CALL wrf_message(TRIM(mes))
WRITE (mes,'(A,F10.3)') 'In dolph: THETAS (STOP-BAND EDGE) = ', thetas
CALL wrf_message(TRIM(mes))
WRITE (mes,'(A,2F10.3)') 'In dolph: R,DB = ', R,DB
CALL wrf_message(TRIM(mes))
DO NT=0,M
SUM = 1
DO I=1,M
ARG = X0*COS(I*PI/N)
CALL CHEBY(T,NM1,ARG)
TERM1 = T(NM1)
TERM2 = COS(2*NT*PI*I/N)
SUM = SUM + R*2*TERM1*TERM2
ENDDO
W(NT) = SUM/N
TIME(NT) = NT
WRITE (mes,'(A,F10.6,2x,E17.7)') 'In dolph: TIME, W = ', TIME(NT), W(NT)
CALL wrf_message(TRIM(mes))
ENDDO
! fill in the negative-time values by symmetry.
DO NT=0,M
W2(M+NT) = W(NT)
W2(M-NT) = W(NT)
ENDDO
! fill up the array for return
SUMW = 0.
DO NT=0,2*M
SUMW = SUMW + W2(NT)
ENDDO
WRITE (mes,'(A,F10.4)') 'In dolph: SUM OF WEIGHTS W2 = ', sumw
CALL wrf_message(TRIM(mes))
DO NT=0,M
WINDOW(NT) = W2(NT)
ENDDO
RETURN
END SUBROUTINE dolph
SUBROUTINE cheby(t, n, x) 2
! calculate all chebyshev polynomials up to order n
! for the argument value x.
! reference: numerical recipes, page 184, recurrence
! t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) , n>=2.
IMPLICIT NONE
! Arguments
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: x
REAL, DIMENSION(0:N) :: t
integer :: nn
T(0) = 1
T(1) = X
IF(N.LT.2) RETURN
DO NN=2,N
T(NN) = 2*X*T(NN-1) - T(NN-2)
ENDDO
RETURN
END SUBROUTINE cheby
SUBROUTINE rhofil(dt, tauc, norder, nstep, ictype, frow, nosize) 1,2
! RHO = recurssive high order.
!
! This routine calculates and returns the
! Last Row, FROW, of the FILTER matrix.
!
! Input Parameters:
! DT : Time Step in seconds
! TAUC : Cut-off period (hours)
! NORDER : Order of QS Filter
! NSTEP : Number of step/Size of row.
! ICTYPE : Initial Conditions
! NOSIZE : Max. side of FROW.
!
! Working Fields:
! ACOEF : X-coefficients of filter
! BCOEF : Y-coefficients of filter
! FILTER : Filter Matrix.
!
! Output Parameters:
! FROW : Last Row of Filter Matrix.
!
! Note: Two types of initial conditions are permitted.
! ICTYPE = 1 : Order increasing each row to NORDER.
! ICTYPE = 2 : Order fixed at NORDER throughout.
!
! DOUBLE PRECISION USED THROUGHOUT.
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION MUC
! N.B. Single Precision for List Parameters.
REAL, intent(in) :: DT,TAUC
! Space for the last row of FILTER.
integer, intent(in) :: norder, nstep, ictype, nosize
REAL , dimension(0:nosize), intent(out):: FROW
! Arrays for rho filter.
integer, PARAMETER :: NOMAX=100
real , dimension(0:NOMAX) :: acoef, bcoef
real , dimension(0:NOMAX,0:NOMAX) :: filter
! Working space.
real , dimension(0:NOMAX) :: alpha, beta
real :: DTT
DTT = ABS(DT)
PI = 2*DASIN(1.D0)
IOTA = CMPLX(0.,1.)
! Filtering Parameters (derived).
THETAC = 2*PI*DTT/(TAUC)
MUC = tan(THETAC/2)
FC = THETAC/(2*PI)
! Clear the arrays.
DO NC=0,NOMAX
ACOEF(NC) = 0.
BCOEF(NC) = 0.
ALPHA(NC) = 0.
BETA (NC) = 0.
FROW (NC) = 0.
DO NR=0,NOMAX
FILTER(NR,NC) = 0.
ENDDO
ENDDO
! Fill up the Filter Matrix.
FILTER(0,0) = 1.
! Get the coefficients of the Filter.
IF ( ICTYPE.eq.2 ) THEN
CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
ENDIF
DO 100 NROW=1,NSTEP
IF ( ICTYPE.eq.1 ) THEN
NORD = MIN(NROW,NORDER)
IF ( NORD.le.NORDER) THEN
CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
ENDIF
ENDIF
DO K=0,NROW
ALPHA(K) = ACOEF(NROW-K)
IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
ENDDO
! Correction for terms of negative index.
IF ( ICTYPE.eq.2 ) THEN
IF ( NROW.lt.NORDER ) THEN
CN = 0.
DO NN=NROW+1,NORDER
CN = CN + (ACOEF(NN)+BCOEF(NN))
ENDDO
ALPHA(0) = ALPHA(0) + CN
ENDIF
ENDIF
! Check sum of ALPHAs and BETAs = 1
SUMAB = 0.
DO NN=0,NROW
SUMAB = SUMAB + ALPHA(NN)
IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
ENDDO
DO KK=0,NROW-1
SUMBF = 0.
DO LL=0,NROW-1
SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
ENDDO
FILTER(NROW,KK) = ALPHA(KK)+SUMBF
ENDDO
FILTER(NROW,NROW) = ALPHA(NROW)
! Check sum of row elements = 1
SUMROW = 0.
DO NN=0,NROW
SUMROW = SUMROW + FILTER(NROW,NN)
ENDDO
100 CONTINUE
DO NC=0,NSTEP
FROW(NC) = FILTER(NSTEP,NC)
ENDDO
RETURN
END SUBROUTINE rhofil
SUBROUTINE rhocof(nord, nomax, muc, ca, cb) 2,1
! Get the coefficients of the RHO Filter
! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT NONE
! Arguments
integer, intent(in) :: nord, nomax
real, dimension(0:nomax) :: ca, cb
! Functions
double precision, external :: cnr
! Local variables
INTEGER :: nn
COMPLEX :: IOTA
DOUBLE PRECISION :: MUC, ZN
DOUBLE PRECISION :: pi, root2, rn, sigma, gain, sumcof
PI = 2*ASIN(1.)
ROOT2 = SQRT(2.)
IOTA = (0.,1.)
RN = 1./FLOAT(NORD)
SIGMA = 1./( SQRT(2.**RN-1.) )
GAIN = (MUC*SIGMA/(1+MUC*SIGMA))**NORD
ZN = (1-MUC*SIGMA)/(1+MUC*SIGMA)
DO NN=0,NORD
CA(NN) = CNR(NORD,NN)*GAIN
IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)*(-ZN)**NN
ENDDO
! Check sum of coefficients = 1
SUMCOF = 0.
DO NN=0,NORD
SUMCOF = SUMCOF + CA(NN)
IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
ENDDO
RETURN
END SUBROUTINE RHOCOF
DOUBLE PRECISION FUNCTION cnr(n,r) 1
! Binomial Coefficient (n,r).
! IMPLICIT DOUBLE PRECISION(C,X)
IMPLICIT NONE
! Arguments
INTEGER , intent(in) :: n, R
! Local variables
INTEGER :: k
DOUBLE PRECISION :: coeff, xn, xr, xk
IF ( R.eq.0 ) THEN
CNR = 1.0
RETURN
ENDIF
Coeff = 1.0
XN = DFLOAT(N)
XR = DFLOAT(R)
DO K=1,R
XK = DFLOAT(K)
COEFF = COEFF * ( (XN-XR+XK)/XK )
ENDDO
CNR = COEFF
RETURN
END FUNCTION cnr
SUBROUTINE optfil (grid,NH,DELTAT,NHMAX),4
!----------------------------------------------------------------------
! SUBROUTINE optfil (NH,DELTAT,TAUP,TAUS,LPRINT, &
! H,NHMAX)
!
! - Huang and Lynch optimal filter
! Monthly Weather Review, Feb 1993
!----------------------------------------------------------
! Input Parameters in List:
! NH : Half-length of the Filter
! DELTAT : Time-step (in seconds).
! TAUP : Period of pass-band edge (hours).
! TAUS : Period of stop-band edge (hours).
! LPRINT : Logical switch for messages.
! NHMAX : Maximum permitted Half-length.
!
! Output Parameters in List:
! H : Impulse Response.
! DP : Deviation in pass-band (db)
! DS : Deviation in stop-band (db)
!----------------------------------------------------------
!
USE module_domain, ONLY : domain
TYPE(domain) , POINTER :: grid
REAL,DIMENSION( 20) :: EDGE
REAL,DIMENSION( 10) :: FX, WTX, DEVIAT
REAL,DIMENSION(2*NHMAX+1) :: H
logical LPRINT
REAL, INTENT (IN) :: DELTAT
INTEGER, INTENT (IN) :: NH, NHMAX
!
TAUP = 3.
TAUS = 1.5
LPRINT = .true.
!initialize H array
NL=2*NHMAX+1
do 101 n=1,NL
H(n)=0.
101 continue
NFILT = 2*NH+1
print *,' start optfil, NFILT=', nfilt
!
! 930325 PL & XYH : the upper limit is changed from 64 to 128.
IF(NFILT.LE.0 .OR. NFILT.GT.128 ) THEN
WRITE(6,*) 'NH=',NH
CALL wrf_error_fatal (' Sorry, error 1 in call to OPTFIL ')
ENDIF
!
! The following four should always be the same.
JTYPE = 1
NBANDS = 2
!CC JPRINT = 0
LGRID = 16
!
! calculate transition frequencies.
DT = ABS(DELTAT)
FS = DT/(TAUS*3600.)
FP = DT/(TAUP*3600.)
IF(FS.GT.0.5) then
! print *,' FS too large in OPTFIL '
CALL wrf_error_fatal (' FS too large in OPTFIL ')
! return
end if
IF(FP.LT.0.0) then
! print *, ' FP too small in OPTFIL '
CALL wrf_error_fatal (' FP too small in OPTFIL ')
! return
end if
!
! Relative Weights in pass- and stop-bands.
WTP = 1.0
WTS = 1.0
!
!CC NOTE: (FP,FC,FS) is an arithmetic progression, so
!CC (1/FS,1/FC,1/FP) is a harmonic one.
!CC TAUP = 1/( (1/TAUC)-(1/DTAU) )
!CC TAUS = 1/( (1/TAUC)+(1/DTAU) )
!CC TAUC : Cut-off Period (hours).
!CC DTAU : Transition half-width (hours).
!CC FC = 1/TAUC ; DF = 1/DTAU
!CC FP = FC - DF : FS = FC + DF
!
IF ( LPRINT ) THEN
TAUC = 2./((1/TAUS)+(1/TAUP))
DTAU = 2./((1/TAUS)-(1/TAUP))
FC = DT/(TAUC*3600.)
DF = DT/(DTAU*3600.)
WRITE(6,*) ' DT ' , dt
WRITE(6,*) ' TAUS, TAUP ' , TAUS,TAUP
WRITE(6,*) ' TAUC, DTAU ' , TAUC,DTAU
WRITE(6,*) ' FP, FS ' , FP, FS
WRITE(6,*) ' FC, DF ' , FC, DF
WRITE(6,*) ' WTS, WTP ' , WTS, WTP
ENDIF
!
! Fill the control vectors for MCCPAR
EDGE(1) = 0.0
EDGE(2) = FP
EDGE(3) = FS
EDGE(4) = 0.5
FX(1) = 1.0
FX(2) = 0.0
WTX(1) = WTP
WTX(2) = WTS
CALL MCCPAR(NFILT,JTYPE,NBANDS,LPRINT,LGRID, &
EDGE,FX,WTX,DEVIAT, h )
!
! Save the deviations in the pass- and stop-bands.
DP = DEVIAT(1)
DS = DEVIAT(2)
!
! Fill out the array H (only first half filled in MCCPAR).
IF(MOD(NFILT,2).EQ.0) THEN
NHALF = ( NFILT )/2
ELSE
NHALF = (NFILT+1)/2
ENDIF
DO 100 nn=1,NHALF
H(NFILT+1-nn) = h(nn)
100 CONTINUE
! normalize the sums to be unity
sumh = 0
do 150 n=1,NFILT
sumh = sumh + H(n)
150 continue
print *,'SUMH =', sumh
do 200 n=1,NFILT
H(n) = H(n)/sumh
200 continue
do 201 n=1,NFILT
grid%hcoeff(n)=H(n)
201 continue
! print *,'HCOEFF(n) ', grid%hcoeff
!
END SUBROUTINE optfil
SUBROUTINE MCCPAR (NFILT,JTYPE,NBANDS,LPRINT,LGRID, & 1,7
EDGE,FX,WTX,DEVIAT,h )
! PROGRAM FOR THE DESIGN OF LINEAR PHASE FINITE IMPULSE
! REPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM
!
!************************************************************
!* Reference: McClellan, J.H., T.W. Parks and L.R.Rabiner, *
!* 1973: A computer program for designing *
!* optimum FIR linear phase digital filters. *
!* IEEE Trans. on Audio and Electroacoustics, *
!* Vol AU-21, No. 6, 506-526. *
!************************************************************
!
! THREE TYPES OF FILTERS ARE INCLUDED -- BANDPASS FILTERS
! DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS
!
!---------------------------------------------------------------
!
! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION H(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
DOUBLE PRECISION PI2,PI
DOUBLE PRECISION AD,DEV,X,Y
LOGICAL LPRINT
PI = 3.141592653589793
PI2 = 6.283185307179586
! ......
NFMAX = 128
100 CONTINUE
! PROGRAM INPUT SECTION
!CC READ(5,*) NFILT,JTYPE,NBANDS,JPRINT,LGRID
IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) THEN
CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
END IF
IF(NBANDS.LE.0) NBANDS = 1
! ....
IF(LGRID.LE.0) LGRID = 16
JB = 2*NBANDS
!cc READ(5,*) (EDGE(J),J=1,JB)
!cc READ(5,*) (FX(J),J=1,NBANDS)
!cc READ(5,*) (WTX(J),J=1,NBANDS)
IF(JTYPE.EQ.0) THEN
CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
END IF
NEG = 1
IF(JTYPE.EQ.1) NEG = 0
NODD = NFILT/2
NODD = NFILT-2*NODD
NFCNS = NFILT/2
IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS = NFCNS+1
! ...
GRID(1) = EDGE(1)
DELF = LGRID*NFCNS
DELF = 0.5/DELF
IF(NEG.EQ.0) GOTO 135
IF(EDGE(1).LT.DELF) GRID(1) = DELF
135 CONTINUE
J = 1
L = 1
LBAND = 1
140 FUP = EDGE(L+1)
145 TEMP = GRID(J)
! ....
DES(J) = EFF(TEMP,FX,WTX,LBAND,JTYPE)
WT(J) = WATE(TEMP,FX,WTX,LBAND,JTYPE)
J = J+1
GRID(J) = TEMP+DELF
IF(GRID(J).GT.FUP) GOTO 150
GOTO 145
150 GRID(J-1) = FUP
DES(J-1) = EFF(FUP,FX,WTX,LBAND,JTYPE)
WT(J-1) = WATE(FUP,FX,WTX,LBAND,JTYPE)
LBAND = LBAND+1
L = L+2
IF(LBAND.GT.NBANDS) GOTO 160
GRID(J) = EDGE(L)
GOTO 140
160 NGRID = J-1
IF(NEG.NE.NODD) GOTO 165
IF(GRID(NGRID).GT.(0.5-DELF)) NGRID = NGRID-1
165 CONTINUE
! ......
IF(NEG) 170,170,180
170 IF(NODD.EQ.1) GOTO 200
DO 175 J=1,NGRID
CHANGE = DCOS(PI*GRID(J))
DES(J) = DES(J)/CHANGE
WT(J) = WT(J)*CHANGE
175 CONTINUE
GOTO 200
180 IF(NODD.EQ.1) GOTO 190
DO 185 J = 1,NGRID
CHANGE = DSIN(PI*GRID(J))
DES(J) = DES(J)/CHANGE
WT(J) = WT(J)*CHANGE
185 CONTINUE
GOTO 200
190 DO 195 J =1,NGRID
CHANGE = DSIN(PI2*GRID(J))
DES(J) = DES(J)/CHANGE
WT(J) = WT(J)*CHANGE
195 CONTINUE
! ......
200 TEMP = FLOAT(NGRID-1)/FLOAT(NFCNS)
DO 210 J = 1,NFCNS
IEXT(J) = (J-1)*TEMP+1
210 CONTINUE
IEXT(NFCNS+1) = NGRID
NM1 = NFCNS-1
NZ = NFCNS+1
! CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM
CALL REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
! CALCULATE THE IMPULSE RESPONSE
IF(NEG) 300,300,320
300 IF(NODD.EQ.0) GOTO 310
DO 305 J=1,NM1
H(J) = 0.5*ALPHA(NZ-J)
305 CONTINUE
H(NFCNS)=ALPHA(1)
GOTO 350
310 H(1) = 0.25*ALPHA(NFCNS)
DO 315 J = 2,NM1
H(J) = 0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
315 CONTINUE
H(NFCNS) = 0.5*ALPHA(1)+0.25*ALPHA(2)
GOTO 350
320 IF(NODD.EQ.0) GOTO 330
H(1) = 0.25*ALPHA(NFCNS)
H(2) = 0.25*ALPHA(NM1)
DO 325 J = 3,NM1
H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
325 CONTINUE
H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(3)
H(NZ) = 0.0
GOTO 350
330 H(1) = 0.25*ALPHA(NFCNS)
DO 335 J =2,NM1
H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
335 CONTINUE
H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(2)
! PROGRAM OUTPUT SECTION
350 CONTINUE
!
IF(LPRINT) THEN
print *, '****************************************************'
print *, 'FINITE IMPULSE RESPONSE (FIR)'
print *, 'LINEAR PHASE DIGITAL FILTER DESIGN'
print *, 'REMEZ EXCHANGE ALGORITHM'
IF(JTYPE.EQ.1) WRITE(6,365)
365 FORMAT(25X,'BANDPASS FILTER'/)
IF(JTYPE.EQ.2) WRITE(6,370)
370 FORMAT(25X,'DIFFERENTIATOR '/)
IF(JTYPE.EQ.3) WRITE(6,375)
375 FORMAT(25X,'HILBERT TRANSFORMER '/)
WRITE(6,378) NFILT
378 FORMAT(15X,'FILTER LENGTH =',I3/)
WRITE(6,380)
380 FORMAT(15X,'***** IMPULSE RESPONSE *****')
DO 381 J = 1,NFCNS
K = NFILT+1-J
IF(NEG.EQ.0) WRITE(6,382) J,H(J),K
IF(NEG.EQ.1) WRITE(6,383) J,H(J),K
381 CONTINUE
382 FORMAT(20X,'H(',I3,') = ',E15.8,' = H(',I4,')')
383 FORMAT(20X,'H(',I3,') = ',E15.8,' = -H(',I4,')')
IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ
384 FORMAT(20X,'H(',I3,') = 0.0')
DO 450 K=1,NBANDS,4
KUP = K+3
IF(KUP.GT.NBANDS) KUP = NBANDS
print *
WRITE(6,385) (J,J=K,KUP)
385 FORMAT(24X,4('BAND',I3,8X))
WRITE(6,390) (EDGE(2*J-1),J=K,KUP)
390 FORMAT(2X,'LOWER BAND EDGE',5F15.8)
WRITE(6,395) (EDGE(2*J),J=K,KUP)
395 FORMAT(2X,'UPPER BAND EDGE',5F15.8)
IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
400 FORMAT(2X,'DESIRED VALUE',2X,5F15.8)
IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
405 FORMAT(2X,'DESIRED SLOPE',2X,5F15.8)
WRITE(6,410) (WTX(J),J=K,KUP)
410 FORMAT(2X,'WEIGHTING',6X,5F15.8)
DO 420 J = K,KUP
DEVIAT(J) = DEV/WTX(J)
420 CONTINUE
WRITE(6,425) (DEVIAT(J),J=K,KUP)
425 FORMAT(2X,'DEVIATION',6X,5F15.8)
IF(JTYPE.NE.1) GOTO 450
DO 430 J = K,KUP
DEVIAT(J) = 20.0*ALOG10(DEVIAT(J))
430 CONTINUE
WRITE(6,435) (DEVIAT(J),J=K,KUP)
435 FORMAT(2X,'DEVIATION IN DB',5F15.8)
450 CONTINUE
print *, 'EXTREMAL FREQUENCIES'
WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
455 FORMAT((2X,5F15.7))
WRITE(6,460)
460 FORMAT(1X,70(1H*))
!
ENDIF
!
!CC IF(NFILT.NE.0) GOTO 100 ! removal of re-run loop.
!
END SUBROUTINE mccpar
FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE) 2
DIMENSION FX(5),WTX(5)
IF(JTYPE.EQ.2) GOTO 1
EFF = FX(LBAND)
RETURN
1 EFF = FX(LBAND)*TEMP
END FUNCTION eff
FUNCTION WATE(TEMP,FX,WTX,LBAND,JTYPE) 2
DIMENSION FX(5),WTX(5)
IF(JTYPE.EQ.2) GOTO 1
WATE = WTX(LBAND)
RETURN
1 IF(FX(LBAND).LT.0.0001) GOTO 2
WATE = WTX(LBAND)/TEMP
RETURN
2 WATE = WTX(LBAND)
END FUNCTION wate
! SUBROUTINE ERROR
!! WRITE(6,*)' **** ERROR IN INPUT DATA ****'
! CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
! END SUBROUTINE error
SUBROUTINE REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID) 1,9
! THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
! FOR THE WEIGHTED CHEBCHEV APPROXIMATION OF A CONTINUOUS
! FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE
! ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
! DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
! GRID, THE NUMBER OF COSINES, AND THE INITIAL GUESS OF THE
! EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV
! ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
! FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
! THE COEFFICIENTS OF THE BEST APPROXIMATION.
!
! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
DIMENSION EDGE(20)
DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DIMENSION A(66),P(65),Q(65)
DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION, EXTERNAL :: D, GEE
!
! THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
!
ITRMAX=25
DEVL=-1.0
NZ=NFCNS+1
NZZ=NFCNS+2
NITER=0
100 CONTINUE
IEXT(NZZ)=NGRID+1
NITER=NITER+1
IF(NITER.GT.ITRMAX) GO TO 400
DO 110 J=1,NZ
DTEMP=GRID(IEXT(J))
DTEMP=DCOS(DTEMP*PI2)
110 X(J)=DTEMP
JET=(NFCNS-1)/15+1
DO 120 J=1,NZ
120 AD(J)=D(J,NZ,JET,X)
DNUM=0.0
DDEN=0.0
K=1
DO 130 J=1,NZ
L=IEXT(J)
DTEMP=AD(J)*DES(L)
DNUM=DNUM+DTEMP
DTEMP=K*AD(J)/WT(L)
DDEN=DDEN+DTEMP
130 K=-K
DEV=DNUM/DDEN
NU=1
IF(DEV.GT.0.0) NU=-1
DEV=-NU*DEV
K=NU
DO 140 J=1,NZ
L=IEXT(J)
DTEMP=K*DEV/WT(L)
Y(J)=DES(L)+DTEMP
140 K=-K
IF(DEV.GE.DEVL) GO TO 150
WRITE(6,*) ' ******** FAILURE TO CONVERGE *********** '
WRITE(6,*) ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR '
WRITE(6,*) ' THE IMPULSE RESPONSE MAY BE CORRECT '
WRITE(6,*) ' CHECK WITH A FREQUENCY RESPONSE '
WRITE(6,*) ' **************************************** '
GO TO 400
150 DEVL=DEV
JCHNGE=0
K1=IEXT(1)
KNZ=IEXT(NZ)
KLOW=0
NUT=-NU
J=1
!
! SEARCH FOR THE EXTERMAL FREQUENCIES OF THE BEST
! APPROXIMATION.
200 IF(J.EQ.NZZ) YNZ=COMP
IF(J.GE.NZZ) GO TO 300
KUP=IEXT(J+1)
L=IEXT(J)+1
NUT=-NUT
IF(J.EQ.2) Y1=COMP
COMP=DEV
IF(L.GE.KUP) GO TO 220
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 220
COMP=NUT*ERR
210 L=L+1
IF(L.GE.KUP) GO TO 215
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 215
COMP=NUT*ERR
GO TO 210
215 IEXT(J)=L-1
J=J+1
KLOW=L-1
JCHNGE=JCHNGE+1
GO TO 200
220 L=L-1
225 L=L-1
IF(L.LE.KLOW) GO TO 250
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.GT.0.0) GO TO 230
IF(JCHNGE.LE.0) GO TO 225
GO TO 260
230 COMP=NUT*ERR
235 L=L-1
IF(L.LE.KLOW) GO TO 240
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 240
COMP=NUT*ERR
GO TO 235
240 KLOW=IEXT(J)
IEXT(J)=L+1
J=J+1
JCHNGE=JCHNGE+1
GO TO 200
250 L=IEXT(J)+1
IF(JCHNGE.GT.0) GO TO 215
255 L=L+1
IF(L.GE.KUP) GO TO 260
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 255
COMP=NUT*ERR
GO TO 210
260 KLOW=IEXT(J)
J=J+1
GO TO 200
300 IF(J.GT.NZZ) GO TO 320
IF(K1.GT.IEXT(1)) K1=IEXT(1)
IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
NUT1=NUT
NUT=-NU
L=0
KUP=K1
COMP=YNZ*(1.00001)
LUCK=1
310 L=L+1
IF(L.GE.KUP) GO TO 315
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 310
COMP=NUT*ERR
J=NZZ
GO TO 210
315 LUCK=6
GO TO 325
320 IF(LUCK.GT.9) GO TO 350
IF(COMP.GT.Y1) Y1=COMP
K1=IEXT(NZZ)
325 L=NGRID+1
KLOW=KNZ
NUT=-NUT1
COMP=Y1*(1.00001)
330 L=L-1
IF(L.LE.KLOW) GO TO 340
ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
ERR=(ERR-DES(L))*WT(L)
DTEMP=NUT*ERR-COMP
IF(DTEMP.LE.0.0) GO TO 330
J=NZZ
COMP=NUT*ERR
LUCK=LUCK+10
GO TO 235
340 IF(LUCK.EQ.6) GO TO 370
DO 345 J=1,NFCNS
345 IEXT(NZZ-J)=IEXT(NZ-J)
IEXT(1)=K1
GO TO 100
350 KN=IEXT(NZZ)
DO 360 J=1,NFCNS
360 IEXT(J)=IEXT(J+1)
IEXT(NZ)=KN
GO TO 100
370 IF(JCHNGE.GT.0) GO TO 100
!
! CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
! USING THE INVERSE DISCRETE FOURIER TRANSFORM.
!
400 CONTINUE
NM1=NFCNS-1
FSH=1.0E-06
GTEMP=GRID(1)
X(NZZ)=-2.0
CN=2*NFCNS-1
DELF=1.0/CN
L=1
KKK=0
IF(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1
IF(NFCNS.LE.3) KKK=1
IF(KKK.EQ.1) GO TO 405
DTEMP=DCOS(PI2*GRID(1))
DNUM=DCOS(PI2*GRID(NGRID))
AA=2.0/(DTEMP-DNUM)
BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
405 CONTINUE
DO 430 J=1,NFCNS
FT=(J-1)*DELF
XT=DCOS(PI2*FT)
IF(KKK.EQ.1) GO TO 410
XT=(XT-BB)/AA
! original : FT=ARCOS(XT)/PI2
FT=ACOS(XT)/PI2
410 XE=X(L)
IF(XT.GT.XE) GO TO 420
IF((XE-XT).LT.FSH) GO TO 415
L=L+1
GO TO 410
415 A(J)=Y(L)
GO TO 425
420 IF((XT-XE).LT.FSH) GO TO 415
GRID(1)=FT
A(J)=GEE(1,NZ,GRID,PI2,X,Y,AD)
425 CONTINUE
IF(L.GT.1) L=L-1
430 CONTINUE
GRID(1)=GTEMP
DDEN=PI2/CN
DO 510 J=1,NFCNS
DTEMP=0.0
DNUM=(J-1)*DDEN
IF(NM1.LT.1) GO TO 505
DO 500 K=1,NM1
500 DTEMP=DTEMP+A(K+1)*DCOS(DNUM*K)
505 DTEMP=2.0*DTEMP+A(1)
510 ALPHA(J)=DTEMP
DO 550 J=2,NFCNS
550 ALPHA(J)=2*ALPHA(J)/CN
ALPHA(1)=ALPHA(1)/CN
IF(KKK.EQ.1) GO TO 545
P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
P(2)=2.0*AA*ALPHA(NFCNS)
Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
DO 540 J=2,NM1
IF(J.LT.NM1) GO TO 515
AA=0.5*AA
BB=0.5*BB
515 CONTINUE
P(J+1)=0.0
DO 520 K=1,J
A(K)=P(K)
520 P(K)=2.0*BB*A(K)
P(2)=P(2)+A(1)*2.0*AA
JM1=J-1
DO 525 K=1,JM1
525 P(K)=P(K)+Q(K)+AA*A(K+1)
JP1=J+1
DO 530 K=3,JP1
530 P(K)=P(K)+AA*A(K-1)
IF(J.EQ.NM1) GO TO 540
DO 535 K=1,J
535 Q(K)=-A(K)
Q(1)=Q(1)+ALPHA(NFCNS-1-J)
540 CONTINUE
DO 543 J=1,NFCNS
543 ALPHA(J)=P(J)
545 CONTINUE
IF(NFCNS.GT.3) RETURN
ALPHA(NFCNS+1)=0.0
ALPHA(NFCNS+2)=0.0
END SUBROUTINE remez
DOUBLE PRECISION FUNCTION D(K,N,M,X) 23
! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION Q
DOUBLE PRECISION PI2
D = 1.0
Q = X(K)
DO 3 L = 1,M
DO 2 J = L,N,M
IF(J-K) 1,2,1
1 D = 2.0*D*(Q-X(J))
2 CONTINUE
3 CONTINUE
D = 1.0/D
END FUNCTION D
DOUBLE PRECISION FUNCTION GEE(K,N,GRID,PI2,X,Y,AD) 8
! COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION P,C,D,XF
DOUBLE PRECISION PI2
P = 0.0
XF = GRID(K)
XF = DCOS(PI2*XF)
D = 0.0
DO 1 J =1,N
C = XF-X(J)
C = AD(J)/C
D = D+C
P = P+C*Y(J)
1 CONTINUE
GEE = P/D
END FUNCTION GEE
!+---+-----------------------------------------------------------------+
! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
! A FUNCTION OF TEMPERATURE AND PRESSURE in Thompson microphysics
!
REAL FUNCTION RSLF(P,T) 7
IMPLICIT NONE
REAL, INTENT(IN):: P, T
REAL:: ESL,X
REAL, PARAMETER:: C0= .611583699E03
REAL, PARAMETER:: C1= .444606896E02
REAL, PARAMETER:: C2= .143177157E01
REAL, PARAMETER:: C3= .264224321E-1
REAL, PARAMETER:: C4= .299291081E-3
REAL, PARAMETER:: C5= .203154182E-5
REAL, PARAMETER:: C6= .702620698E-8
REAL, PARAMETER:: C7= .379534310E-11
REAL, PARAMETER:: C8=-.321582393E-13
X=MAX(-80.,T-273.16)
! ESL=612.2*EXP(17.67*X/(T-29.65))
ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
RSLF=.622*ESL/(P-ESL)
END FUNCTION RSLF
!+---+-----------------------------------------------------------------+
! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
! FUNCTION OF TEMPERATURE AND PRESSURE
!
REAL FUNCTION RSIF(P,T) 3
IMPLICIT NONE
REAL, INTENT(IN):: P, T
REAL:: ESI,X
REAL, PARAMETER:: C0= .609868993E03
REAL, PARAMETER:: C1= .499320233E02
REAL, PARAMETER:: C2= .184672631E01
REAL, PARAMETER:: C3= .402737184E-1
REAL, PARAMETER:: C4= .565392987E-3
REAL, PARAMETER:: C5= .521693933E-5
REAL, PARAMETER:: C6= .307839583E-7
REAL, PARAMETER:: C7= .105785160E-9
REAL, PARAMETER:: C8= .161444444E-12
X=MAX(-80.,T-273.16)
ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
RSIF=.622*ESI/(P-ESI)
! ALTERNATIVE
! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
! supercooled water for atmospheric applications, Q. J. R.
! Meteorol. Soc (2005), 131, pp. 1539-1565.
! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
END FUNCTION RSIF
!+---+-----------------------------------------------------------------+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI startfwd group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_startfwd_init ( ) 1,4
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
USE module_utility
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_startfwd_init_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_startfwd_init_recurse
END INTERFACE
! Now, setup all nests
CALL dfi_startfwd_init_recurse(head_grid)
CALL set_current_grid_ptr( head_grid )
END SUBROUTINE wrf_dfi_startfwd_init
RECURSIVE SUBROUTINE dfi_startfwd_init_recurse(grid) 2,5
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_startfwd_init(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_startfwd_init
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
grid_ptr%dt = abs(grid_ptr%dt)
grid_ptr%time_step = abs(grid_ptr%time_step)
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_startfwd_init( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_startfwd_init_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_startfwd_init_recurse
SUBROUTINE dfi_startfwd_init ( grid ) 1,3
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
USE module_utility
USE module_state_description
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER rc
INTERFACE
SUBROUTINE Setup_Timekeeping(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE Setup_Timekeeping
END INTERFACE
grid%dfi_stage = DFI_STARTFWD
#if (EM_CORE == 1)
! No need for adaptive time-step
CALL nl_set_use_adaptive_time_step( grid%id, .false. )
#endif
CALL Setup_Timekeeping (grid)
grid%start_subtime = domain_get_start_time ( head_grid )
grid%stop_subtime = domain_get_stop_time ( head_grid )
CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
END SUBROUTINE dfi_startfwd_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI startbck group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_startbck_init ( ),4
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
USE module_utility
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_startbck_init_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_startbck_init_recurse
END INTERFACE
! Now, setup all nests
CALL dfi_startbck_init_recurse(head_grid)
CALL set_current_grid_ptr( head_grid )
END SUBROUTINE wrf_dfi_startbck_init
RECURSIVE SUBROUTINE dfi_startbck_init_recurse(grid) 2,5
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_startbck_init(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_startbck_init
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
grid_ptr%dt = abs(grid_ptr%dt)
grid_ptr%time_step = abs(grid_ptr%time_step)
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_startbck_init( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_startbck_init_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_startbck_init_recurse
SUBROUTINE dfi_startbck_init ( grid ) 1,4
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
USE module_utility
USE module_state_description
IMPLICIT NONE
TYPE (domain) , POINTER :: grid
INTEGER rc
INTERFACE
SUBROUTINE Setup_Timekeeping(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE Setup_Timekeeping
END INTERFACE
grid%dfi_stage = DFI_STARTBCK
! set physics options to zero
CALL nl_set_mp_physics( grid%id, 0 )
CALL nl_set_ra_lw_physics( grid%id, 0 )
CALL nl_set_ra_sw_physics( grid%id, 0 )
CALL nl_set_sf_surface_physics( grid%id, 0 )
CALL nl_set_sf_sfclay_physics( grid%id, 0 )
CALL nl_set_sf_urban_physics( grid%id, 0 )
CALL nl_set_bl_pbl_physics( grid%id, 0 )
CALL nl_set_cu_physics( grid%id, 0 )
CALL nl_set_damp_opt( grid%id, 0 )
CALL nl_set_sst_update( grid%id, 0 )
CALL nl_set_gwd_opt( grid%id, 0 )
#if (EM_CORE == 1)
CALL nl_set_diff_6th_opt( grid%id, 0 )
CALL nl_set_use_adaptive_time_step( grid%id, .false. )
#endif
CALL nl_set_feedback( grid%id, 0 )
#if (EM_CORE == 1)
! set bc
CALL nl_set_constant_bc( grid%id, head_grid%constant_bc)
#endif
#ifdef WRF_CHEM
! set chemistry option to zero
CALL nl_set_chem_opt (grid%id, 0)
CALL nl_set_aer_ra_feedback (grid%id, 0)
CALL nl_set_io_form_auxinput5 (grid%id, 0)
CALL nl_set_io_form_auxinput7 (grid%id, 0)
CALL nl_set_io_form_auxinput8 (grid%id, 0)
#endif
#if (EM_CORE == 1)
! set diffusion to zero for backward integration
CALL nl_set_km_opt( grid%id, grid%km_opt_dfi)
CALL nl_set_moist_adv_dfi_opt( grid%id, grid%moist_adv_dfi_opt)
IF ( grid%moist_adv_opt == 2 ) THEN
CALL nl_set_moist_adv_opt( grid%id, 0)
ENDIF
#endif
! If a request to do pressure level diags, then shut it off
! until we get to the real forecast part of this integration.
#if (EM_CORE == 1)
CALL nl_set_p_lev_diags( grid%id, grid%p_lev_diags_dfi)
#endif
!tgs need to call start_domain here to reset bc initialization for
! negative dt, but only for outer domain.
if (grid%id == 1) then
CALL start_domain ( grid , .TRUE. )
endif
! Call wrf_run to advance forward 1 step
! Negate time step
CALL nl_set_time_step ( grid%id, -grid%time_step )
CALL Setup_Timekeeping (grid)
grid%start_subtime = domain_get_start_time ( grid )
grid%stop_subtime = domain_get_stop_time ( grid )
CALL WRFU_ClockSet(grid%domain_clock, currTime=grid%start_subtime, rc=rc)
END SUBROUTINE dfi_startbck_init
SUBROUTINE wrf_dfi_bck_init ( ) 2,3
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time
USE module_utility
USE module_state_description
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_bck_init_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_bck_init_recurse
END INTERFACE
! We can only call dfi_bck_init for the head_grid
! since nests have not been instantiated at this point,
! so, dfi_bck_init will need to be called for each
! nest from integrate.
CALL dfi_bck_init_recurse(head_grid)
END SUBROUTINE wrf_dfi_bck_init
RECURSIVE SUBROUTINE dfi_bck_init_recurse(grid) 2,5
USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_bck_init(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_bck_init
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
grid_ptr%dt = abs(grid_ptr%dt)
grid_ptr%time_step = abs(grid_ptr%time_step)
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_bck_init( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_bck_init_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_bck_init_recurse
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI forward initialization group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_fwd_init ( ) 3,4
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
USE module_utility
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_fwd_init_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_fwd_init_recurse
END INTERFACE
! Now, setup all nests
CALL dfi_fwd_init_recurse(head_grid)
CALL set_current_grid_ptr( head_grid )
END SUBROUTINE wrf_dfi_fwd_init
RECURSIVE SUBROUTINE dfi_fwd_init_recurse(grid) 2,5
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_fwd_init(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_fwd_init
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
grid_ptr%dt = abs(grid_ptr%dt)
grid_ptr%time_step = abs(grid_ptr%time_step)
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_fwd_init( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_fwd_init_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_fwd_init_recurse
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI forecast initialization group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_fst_init ( ) 3,4
USE module_domain, ONLY : domain, head_grid, domain_get_stop_time, domain_get_start_time, set_current_grid_ptr
USE module_utility
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_fst_init_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_fst_init_recurse
END INTERFACE
! Now, setup all nests
CALL dfi_fst_init_recurse(head_grid)
CALL set_current_grid_ptr( head_grid )
END SUBROUTINE wrf_dfi_fst_init
RECURSIVE SUBROUTINE dfi_fst_init_recurse ( grid ) 2,5
USE module_domain, ONLY : domain, domain_get_stop_time, domain_get_start_time, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_fst_init(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_fst_init
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
grid_ptr%dt = abs(grid_ptr%dt)
grid_ptr%time_step = abs(grid_ptr%time_step)
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_fst_init( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_fst_init_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_fst_init_recurse
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DFI write initialization group of functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE wrf_dfi_write_initialized_state( ) 3,3
USE module_domain, ONLY : domain, head_grid
INTERFACE
SUBROUTINE dfi_write_initialized_state_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_write_initialized_state_recurse
END INTERFACE
! Now, setup all nests
CALL dfi_write_initialized_state_recurse(head_grid)
END SUBROUTINE wrf_dfi_write_initialized_state
RECURSIVE SUBROUTINE dfi_write_initialized_state_recurse( grid ) 2,4
USE module_domain, ONLY : domain, max_nests
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_write_initialized_state( grid )
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_write_initialized_state
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
CALL dfi_write_initialized_state( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_write_initialized_state_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_write_initialized_state_recurse
SUBROUTINE wrf_tdfi_write_analyzed_state( ),3
USE module_domain, ONLY : domain, head_grid
INTERFACE
SUBROUTINE tdfi_write_analyzed_state_recurse(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE tdfi_write_analyzed_state_recurse
END INTERFACE
! Now, setup all nests
CALL tdfi_write_analyzed_state_recurse(head_grid)
END SUBROUTINE wrf_tdfi_write_analyzed_state
RECURSIVE SUBROUTINE tdfi_write_analyzed_state_recurse( grid ) 2,4
USE module_domain, ONLY : domain, max_nests
IMPLICIT NONE
INTERFACE
SUBROUTINE tdfi_write_analyzed_state( grid )
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE tdfi_write_analyzed_state
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
!
! Assure that time-step is set back to positive
! for this forward step.
!
CALL tdfi_write_analyzed_state( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL tdfi_write_analyzed_state_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE tdfi_write_analyzed_state_recurse
RECURSIVE SUBROUTINE dfi_array_reset_recurse(grid) 2,5
USE module_domain, ONLY : domain, max_nests, set_current_grid_ptr
IMPLICIT NONE
INTERFACE
SUBROUTINE dfi_array_reset(grid)
USE module_domain, ONLY : domain
TYPE (domain), POINTER :: grid
END SUBROUTINE dfi_array_reset
END INTERFACE
INTEGER :: kid
TYPE (domain), POINTER :: grid
TYPE (domain), POINTER :: grid_ptr
grid_ptr => grid
DO WHILE ( ASSOCIATED( grid_ptr ) )
CALL set_current_grid_ptr( grid_ptr )
CALL dfi_array_reset( grid_ptr )
DO kid = 1, max_nests
IF ( ASSOCIATED( grid_ptr%nests(kid)%ptr ) ) THEN
CALL dfi_array_reset_recurse(grid_ptr%nests(kid)%ptr)
ENDIF
END DO
grid_ptr => grid_ptr%sibling
END DO
END SUBROUTINE dfi_array_reset_recurse
!-------------------------------------------------------------------
#if (EM_CORE == 1)
SUBROUTINE rebalance_driver_dfi ( grid ) 2,2
USE module_domain, ONLY : domain
IMPLICIT NONE
TYPE (domain) :: grid
CALL rebalance_dfi( grid &
!
#include "actual_new_args.inc"
!
)
END SUBROUTINE rebalance_driver_dfi
!---------------------------------------------------------------------
SUBROUTINE rebalance_dfi ( grid & 1,5
!
#include "dummy_new_args.inc"
!
)
USE module_domain, ONLY : domain
USE module_utility
USE module_state_description
USE module_model_constants
USE module_driver_constants
#ifdef DM_PARALLEL
USE module_dm
USE module_comm_dm, ONLY : &
HALO_EM_INIT_1_sub &
,HALO_EM_INIT_2_sub &
,HALO_EM_INIT_3_sub &
,HALO_EM_INIT_4_sub &
,HALO_EM_INIT_5_sub
#endif
IMPLICIT NONE
TYPE (domain) :: grid
#include "dummy_new_decl.inc"
REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold
REAL :: qvf , qvf1 , qvf2, qtot
REAL :: pfu, pfd, phm
! Local domain indices and counters.
INTEGER :: n_moist
INTEGER :: &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
ips, ipe, jps, jpe, kps, kpe, &
i, j, k, ispe, ktf
SELECT CASE ( model_data_order )
CASE ( DATA_ORDER_ZXY )
kds = grid%sd31 ; kde = grid%ed31 ;
ids = grid%sd32 ; ide = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
kms = grid%sm31 ; kme = grid%em31 ;
ims = grid%sm32 ; ime = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XYZ )
ids = grid%sd31 ; ide = grid%ed31 ;
jds = grid%sd32 ; jde = grid%ed32 ;
kds = grid%sd33 ; kde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
jms = grid%sm32 ; jme = grid%em32 ;
kms = grid%sm33 ; kme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XZY )
ids = grid%sd31 ; ide = grid%ed31 ;
kds = grid%sd32 ; kde = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
kms = grid%sm32 ; kme = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
END SELECT
ktf=MIN(kte,kde-1)
DO j=jts,jte
DO i=its,ite
grid%ph_2(i,1,j) = 0.
END DO
END DO
! Base state potential temperature and inverse density (alpha = 1/rho) from
! the half eta levels and the base-profile surface pressure. Compute 1/rho
! from equation of state. The potential temperature is a perturbation from t0.
n_moist = num_moist-1
! print *,'n_moist,PARAM_FIRST_SCALAR',n_moist,PARAM_FIRST_SCALAR
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
! equation) down from the top to get the pressure perturbation. First get the pressure
! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
k = kte-1
qtot = 0.
DO ispe=PARAM_FIRST_SCALAR,n_moist
qtot = qtot + 0.5*(moist(i,k,j,ispe)+moist(i,k,j,ispe))
ENDDO
qvf2 = 1./(1.+qtot)
qvf1 = qtot*qvf2
grid%p(i,k,j) = - 0.5*(grid%mu_2(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
qvf = 1.+rvovrd*moist(i,k,j,P_QV)
grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
! inverse density fields (total and perturbation).
DO k=kte-2,1,-1
qtot = 0.
DO ispe=PARAM_FIRST_SCALAR,n_moist
qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
ENDDO
qvf2 = 1./(1.+qtot)
qvf1 = qtot*qvf2
grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_2(i,j) + &
qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
qvf = 1. + rvovrd*moist(i,k,j,P_QV)
grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
ENDDO
! This is the hydrostatic equation used in the model after the small timesteps. In
! the model, grid%al (inverse density) is computed from the geopotential.
IF (grid%hypsometric_opt == 1) THEN
DO k = 2,kte
grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
+ grid%mu_2(i,j)*grid%alb(i,k-1,j) )
grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
END DO
ELSE IF (grid%hypsometric_opt == 2) THEN
! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is
! dry pressure.
! Note that al*p approximates Rd*T and dLOG(p) does z.
! Here T varies mostly linear with z, the first-order
! integration produces better result.
grid%ph_2(i,1,j) = grid%phb(i,1,j)
DO k = 2,kte
pfu = grid%mu0(i,j)*grid%znw(k) + grid%p_top
pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
END DO
DO k = 1,kte
grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
END DO
END IF
ENDDO !i
ENDDO !j
ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
#ifdef DM_PARALLEL
# include "HALO_EM_INIT_1.inc"
# include "HALO_EM_INIT_2.inc"
# include "HALO_EM_INIT_3.inc"
# include "HALO_EM_INIT_4.inc"
# include "HALO_EM_INIT_5.inc"
#endif
END SUBROUTINE rebalance_dfi
#endif
!---------------------------------------------------------------------