subroutine da_transfer_wrftoxb(xbx, grid, config_flags)

   !---------------------------------------------------------------------------
   ! Purpose: Transfers fields from WRF to first guess structure.
   !    Updated for Analysis on Arakawa-C grid
   !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   !---------------------------------------------------------------------------

   implicit none
   
   type (xbx_type), intent(inout)     :: xbx        ! Header & non-gridded vars.

   type(domain), intent(inout)        :: grid
   type(grid_config_rec_type), intent(in) :: config_flags

   integer :: i, j, k, ij

   real    :: theta, tmpvar, alt

   real, dimension(ims:ime,jms:jme) :: rgh_fac

   character(len=19) :: current_date

   real :: loc_psac_mean

   real, dimension(jds:jde) :: loc_latc_mean

   integer :: size2d

   real, dimension(kms:kme) :: DDT

   real   :: qvf1, cvpm, cpovcv, ppb, ttb, albn, aln, height, temp
   real, allocatable :: arrayglobal(:,:)

   logical:: no_ppb

   real, dimension(ims:ime,kms:kme) :: pf, pp

#ifdef A2C
   real   :: uu, vv
#endif


   ! If grid%pb does not existed in FG (YRG, 08/26/2010):
     ppb = sum(grid%pb*grid%pb)
     no_ppb = ppb == 0.0

   !---------------------------------------------------------------------------
   ! Set xb array range indices for processor subdomain.
   !---------------------------------------------------------------------------

   if (trace_use) call da_trace_entry("da_transfer_wrftoxb")

   grid%xb % map  = grid%map_proj
   grid%xb % ds   = grid%dx

   grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
   grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
   grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1

   ! WHY?
   !rizvi   xbx%big_header%bhi(5,5) = grid%start_year
   !rizvi   xbx%big_header%bhi(6,5) = grid%start_month
   !rizvi   xbx%big_header%bhi(7,5) = grid%start_day
   !rizvi   xbx%big_header%bhi(8,5) = grid%start_hour

   !---------------------------------------------------------------------------
   ! WRF-specific fitelds:
   !---------------------------------------------------------------------------

   ptop = grid%p_top

   grid%xb%sigmaf(kte+1) = grid%znw(kte+1)

   grid%xb%znw(kte+1) = grid%znw(kte+1)
   grid%xb%znu(kte+1) = 0.0
 
   do k=kts,kte
      grid%xb%sigmah(k) = grid%znu(k)
      grid%xb%sigmaf(k) = grid%znw(k)

      grid%xb%znu(k) = grid%znu(k)
      grid%xb%znw(k) = grid%znw(k)
      grid%xb%dn(k)  = grid%dn(k)
      grid%xb%dnw(k) = grid%dnw(k)
   end do

   grid%xb % ptop = ptop
      
   !---------------------------------------------------------------------------
   ! Convert WRF fields to xb:
   !---------------------------------------------------------------------------

   if (print_detail_xb) then
      write(unit=stdout, fmt='(3a, i8)') &
         'file:', __FILE__, ', line:', __LINE__

      write(unit=stdout, fmt=*) 'its,ite=', its,ite
      write(unit=stdout, fmt=*) 'jts,jte=', jts,jte
      write(unit=stdout, fmt=*) 'kts,kte=', kts,kte

      write(unit=stdout, fmt='(/5a/)') &
         'lvl         dnw                dn             rdnw       rdn'

      do k=kts,kte+1
         write(unit=stdout, fmt='(i3,8f16.8)') k, &
            grid%dnw(k), grid%dn(k), grid%rdnw(k), grid%rdn(k)
      end do

      write(unit=stdout, fmt='(/5a/)') &
         'lvl         znu                 znw           rdnw       rdn'

      do k=kts,kte+1
         write(unit=stdout, fmt='(i3,8f16.8)') k, &
            grid%xb%sigmah(k), grid%xb%sigmaf(k), grid%rdnw(k), &
            grid%rdn(k)
      end do

      write(unit=stdout, fmt='(/5a/)') &
         'lvl         phb                 ph_2'

      do k=kts,kte
         write(unit=stdout, fmt='(i3,8e20.12)') k, &
               grid%phb(its,jts,k), grid%ph_2(its,jts,k)
      end do

      write(unit=stdout, fmt=*) 'simple variables:'

      if (jte == jde) then
         write(unit=stdout, fmt=*) ' '

         do k=kts+5,kte,10
            do i=its,ite,10
               write(unit=stdout, fmt=*) &
                    '  grid%v_2(', i, ',', jde+1, ',', k, ')=', &
                       grid%v_2(i, jde+1,k)
            end do
            write(unit=stdout, fmt=*) ' '
         end do
      end if

      if (ite == ide) then
         write(unit=stdout, fmt=*) ' '

         do k=kts+5,kte,10
            do j=jts,jte,10
               write(unit=stdout, fmt=*) &
                  '  grid%u_2(', ide+1, ',', j, ',', k, ')=', &
                  grid%u_2(ide+1,j,k)
            end do
            write(unit=stdout, fmt=*) ' '
         end do
      end if

      write(unit=stdout, fmt=*) 'simple variables:'

      write(unit=stdout,fmt=*) &
         '  grid%u_1(its,jts,kts)=', grid%u_1(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%v_1(its,jts,kts)=', grid%v_1(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%w_1(its,jts,kts)=', grid%w_1(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%t_1(its,jts,kts)=', grid%t_1(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%ph_1(its,jts,kts)=',grid%ph_1(its,jts,kts)


      write(unit=stdout,fmt=*) &
         '  grid%u_2(its,jte,kts)=', grid%u_2(its,jte,kts)
      write(unit=stdout,fmt=*) &
         '  grid%v_2(ite,jts,kts)=', grid%v_2(ite,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%w_2(its,jts,kts)=', grid%w_2(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%t_2(its,jts,kts)=', grid%t_2(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%ph_2(its,jts,kts)=',grid%ph_2(its,jts,kts)
      write(unit=stdout,fmt=*) &
         '  grid%phb(its,jts,kts)=', grid%phb(its,jts,kts)

      write(unit=stdout, fmt=*) &
         '  grid%sm31,grid%em31,grid%sm32,grid%em32, grid%sm33,grid%em33=', &
         grid%sm31,grid%em31,grid%sm32,grid%em32,grid%sm33,grid%em33

      write(unit=stdout, fmt=*) '  grid%p_top=', grid%p_top
      write(unit=stdout, fmt=*) '  grid%znu(kts)=', grid%znu(kts)
      write(unit=stdout, fmt=*) '  grid%mub(its,jts)=', grid%mub(its,jts)
      write(unit=stdout, fmt=*) '  grid%mu_2(its,jts)=', &
         grid%mu_2(its,jts)

      write(unit=stdout, fmt=*) '  hbot(its,jts)=', grid%hbot(its,jts)
      write(unit=stdout, fmt=*) '  htop(its,jts)=', grid%htop(its,jts)

      write(unit=stdout, fmt=*) '  grid%p_top=', grid%p_top
      write(unit=stdout, fmt=*) '  num_moist=', num_moist
      write(unit=stdout, fmt=*) '  P_QV=', P_QV

      write(unit=stdout, fmt=*) '  moist(its,jts,kts,2)=', &
         grid%moist(its,jts,kts,2)
      write(unit=stdout, fmt=*) ' '
   end if

   !---------------------------------------------------------------
   ! Need this to exchange values in the halo region.
   ! grid%xa%u and grid%xa%v are used as temporary arrays and so
   ! it is easy to use the existing exchange scheme.
   !
   ! Note, this is needed as u_2 and v_2 has no guarantee
   ! the most east column, and the most north row are
   ! properly initailized for each tile.
   !---------------------------------------------------------------
#ifdef A2C

   grid%xa%u(its:ite+1,jts:jte,kts:kte) = grid%u_2(its:ite+1,jts:jte,kts:kte)
   grid%xa%v(its:ite,jts:jte+1,kts:kte) = grid%v_2(its:ite,jts:jte+1,kts:kte)
   grid%xb%u(its:ite+1,jts:jte,kts:kte) = grid%u_2(its:ite+1,jts:jte,kts:kte)
   grid%xb%v(its:ite,jts:jte+1,kts:kte) = grid%v_2(its:ite,jts:jte+1,kts:kte)

!rizvi's update
  if ((fg_format==fg_format_wrf_arw_regional  .or. &
       fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
     ipe = ipe + 1
     ide = ide + 1
  end if

  if ((fg_format==fg_format_wrf_arw_regional  .or. &
       fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
     jpe = jpe + 1
     jde = jde + 1
  end if
!rizvi's update
#ifdef DM_PARALLEL
#include "HALO_XB_UV.inc"
#endif

!rizvi's update
  if ((fg_format==fg_format_wrf_arw_regional  .or. &
       fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
     ipe = ipe - 1
     ide = ide - 1
  end if

  if ((fg_format==fg_format_wrf_arw_regional  .or. &
       fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
     jpe = jpe - 1
     jde = jde - 1
  end if

#else

   ! Fill the halo region for u and v.

#ifdef DM_PARALLEL
#include "HALO_EM_C.inc"
#endif

#endif

   if (print_detail_xb) then
      write(unit=stdout, fmt=*) &
         ' ids,ide,jds,jde,kds,kde=', ids,ide,jds,jde,kds,kde
      write(unit=stdout, fmt=*) &
         ' its,ite,jts,jte,kts,kte=', its,ite,jts,jte,kts,kte
      write(unit=stdout, fmt=*) &
          ' ims,ime,jms,jme,kms,kme=', ims,ime,jms,jme,kms,kme
         
      write(unit=stdout, fmt=*) &
         ' lbound(grid%xb%u)=',   lbound(grid%xb%u)
      write(unit=stdout, fmt=*) &
         ' lbound(grid%xb%v)=',   lbound(grid%xb%v)
      write(unit=stdout, fmt=*) &
         ' lbound(grid%u_2)=', lbound(grid%u_2)
      write(unit=stdout, fmt=*) &
         ' lbound(grid%v_2)=', lbound(grid%v_2)
      write(unit=stdout, fmt=*) &
         ' ubound(grid%xb%u)=',   ubound(grid%xb%u)
      write(unit=stdout, fmt=*) &
         ' ubound(grid%xb%v)=',   ubound(grid%xb%v)
      write(unit=stdout, fmt=*) &
         ' ubound(grid%u_2)=', ubound(grid%u_2)
      write(unit=stdout, fmt=*) &
         ' ubound(grid%v_2)=', ubound(grid%v_2)
   end if

   !$OMP PARALLEL DO &
   !$OMP PRIVATE ( ij, i, j, k, cvpm, cpovcv, ppb, temp, ttb ) &
   !$OMP PRIVATE ( albn, qvf1, aln, theta )
   do ij = 1 , grid%num_tiles

   do j=grid%j_start(ij), grid%j_end(ij)
      k = kte+1

      do i=its,ite
         grid%p(i,j,k) = 0.0
         grid%xb%map_factor(i,j) = grid%msft(i,j)
         grid%xb%cori(i,j) = grid%f(i,j)
         grid%xb%tgrn(i,j) = grid%sst(i,j)
         if (grid%xb%tgrn(i,j) < 100.0) &
            grid%xb%tgrn(i,j) = grid%tmn(i,j)
         grid%xb%lat(i,j) = grid%xlat(i,j)
         grid%xb%lon(i,j) = grid%xlong(i,j)
         grid%xb%terr(i,j) = grid%ht(i,j)
         grid%xb%snow(i,j) = grid%snowc(i,j)
         grid%xb%lanu(i,j) = grid%lu_index(i,j)
         grid%xb%landmask(i,j) = grid%landmask(i,j)
         grid%xb%xland(i,j) = grid%xland(i,j)
         ! Z. Liu below are variables used by RTTOV
         grid%xb%tsk(i,j) = grid%tsk(i,j)
         grid%xb%smois(i,j) = grid%smois(i,j,1)
         grid%xb%tslb(i,j) = grid%tslb(i,j,1)
         grid%xb%xice(i,j) = grid%xice(i,j)
         grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
         grid%xb%isltyp(i,j) = grid%isltyp(i,j)
         grid%xb%vegfra(i,j) = grid%vegfra(i,j)
         grid%xb%snowh(i,j) = grid%snowh(i,j)*1000.0 ! meter to mm    
      end do

      ! WHY?
      ! Adapted the code from "real.init.code" by Y.-R. Guo 05/13/2004:

      ! do i=its,ite
      !    k = kte
      !    qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV))
      !    qvf2 = 1.0/(1.0+qvf1)
      !    qvf1 = qvf1*qvf2
      !    grid%xb%p(i,j,k) = -0.5*(grid%mu_2(i,j)+qvf1* &
      !       grid%mub(i,j))/grid%rdnw(k)/qvf2

      !    do k = kte-1,1,-1
      !       qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV))
      !       qvf2 = 1.0/(1.0+qvf1)
      !       qvf1 = qvf1*qvf2
      !       grid%p(i,j,k) = grid%p(i,j,k+1) - &
      !          (grid%mu_2(i,j)+qvf1*grid%mub(i,j))/qvf2/rdn(k+1)
      !    end do
      ! end do

      ! Adapted the code from WRF module_big_step_utilitites_em.F ----
      !         subroutine calc_p_rho_phi      Y.-R. Guo (10/20/2004)

      ! NOTE: as of V3.1, P (pressure perturbation) and PB (base state pressure)
      ! are included in the wrfinput file. However, P and PB are still
      ! re-calculated here.

      ! NOTE: as of 7/01/2010, P and PB are directly from the first guess
      ! and no longer re-calculated here.

      cvpm =  - (1.0 - gas_constant/cp)
      cpovcv = cp / (cp - gas_constant)

      ! In case of var4d, pb etc. will be recalculated in start_em with realsize=8, 
      ! However, the originals are computed with realsize=4.
      if ( no_ppb ) then
         do k=kts,kte
            do i=its,ite
               ! The base specific volume (from real.init.code)
               ppb  = grid%znu(k) * grid%mub(i,j) + ptop
               grid%pb(i,j,k) = ppb
               temp = MAX ( iso_temp, base_temp + base_lapse*log(ppb/base_pres) )
               ttb  = temp * (base_pres/ppb)**kappa
               ! ttb  = (base_temp + base_lapse*log(ppb/base_pres)) * &
               !   (base_pres/ppb)**kappa
               albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm

               qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
               aln  = -1.0 / (grid%mub(i,j)+grid%mu_2(i,j)) * &
                      (albn*grid%mu_2(i,j) + grid%rdnw(k) * &
                      (grid%ph_2(i,j,k+1) - grid%ph_2(i,j,k)))
               ! total pressure:
               grid%xb%p(i,j,k) = base_pres * &
                                 ((gas_constant*(t0+grid%t_2(i,j,k))*qvf1) / &
                                 (base_pres*(aln+albn)))**cpovcv
               ! total density
               grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
               ! pressure purtubation:
               grid%p(i,j,k) = grid%xb%p(i,j,k) - ppb
            end do
         end do
      else
         do k=kts,kte
            do i=its,ite
              if (grid%hypsometric_opt == 1) then
                ppb = grid%pb(i,j,k)
                temp = MAX ( iso_temp, base_temp + base_lapse*log(ppb/base_pres) )
                ttb  = temp * (base_pres/ppb)**kappa
                ! ttb  = (base_temp + base_lapse*log(ppb/base_pres)) * &
                !   (base_pres/ppb)**kappa
                albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm

                qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
                aln  = -1.0 / (grid%mub(i,j)+grid%mu_2(i,j)) * &
                       (albn*grid%mu_2(i,j) + grid%rdnw(k) * &
                       (grid%ph_2(i,j,k+1) - grid%ph_2(i,j,k)))
                grid%xb%p(i,j,k) = grid%pb(i,j,k) + grid%p(i,j,k)
                ! total density
                grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
              elseif (grid%hypsometric_opt == 2) then
                grid%xb%p(i,j,k) = grid%pb(i,j,k) + grid%p(i,j,k)
                theta = t0 + grid%t_2(i,j,k)
                qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
                ! inverse density via the equation of state
                alt = gas_constant*theta*qvf1/base_pres*(grid%xb%p(i,j,k)/base_pres)**cvpm         
                ! total density
                grid%xb%rho(i,j,k)= 1.0/alt
              endif
            end do
         end do
      endif

      do k=kts,kte+1
         do i=its,ite
            grid%xb%hf(i,j,k) = (grid%phb(i,j,k)+grid%ph_2(i,j,k))/gravity
            grid%xb%w (i,j,k) = grid%w_2(i,j,k)
         end do
      end do

      if (grid%hypsometric_opt == 2) then
        ! Compute full-level pressure
        ! The full and half level dry hydrostatic pressure is easily computed (YRG, 12/20/2011):
        do i=its,ite
           do k=kts, kte+1
              pf(i,k) = (grid%mub(i,j)+grid%mu_2(i,j)) * grid%znw(k) + ptop
              pp(i,k) = (grid%mub(i,j)+grid%mu_2(i,j)) * grid%znu(k) + ptop
           enddo
        enddo
      endif

      do k=kts,kte
         do i=its,ite
            grid%xb%u(i,j,k) = 0.5*(grid%u_2(i,j,k)+grid%u_2(i+1,j,k))
            grid%xb%v(i,j,k) = 0.5*(grid%v_2(i,j,k)+grid%v_2(i,j+1,k))
            grid%xb%wh(i,j,k)= 0.5*(grid%xb%w(i,j,k)+grid%xb%w(i,j,k+1))
            if (grid%hypsometric_opt == 1) then
               grid%xb%h(i,j,k) = 0.5*(grid%xb%hf(i,j,k)+grid%xb%hf(i,j,k+1))
            elseif (grid%hypsometric_opt == 2) then
            !  This is almost correct for pressure, but not for height. 
            !  Arithmetic mean of full-level heights is always higher than the actual,
            !  leading to a biased model for height-based observations (e.g., GPS RO)
            !  and surface variables (2-meter TQ and 10-meter wind).
            !  wee 11/22/2011

               grid%xb%h(i,j,k) = grid%xb%hf(i,j,k)+(grid%xb%hf(i,j,k+1)-grid%xb%hf(i,j,k)) &
                    *LOG(pf(i,k)/pp(i,k))/LOG(pf(i,k)/pf(i,k+1))
            endif

            if ( num_pseudo == 0 ) then
               grid%moist(i,j,k,P_QV) = max(grid%moist(i,j,k,P_QV), qlimit)
            end if
            grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)

            theta = t0 + grid%t_2(i,j,k)
            grid%xb%t(i,j,k) = theta*(grid%xb%p(i,j,k)/base_pres)**kappa

            ! Convert to specific humidity from mixing ratio of water vapor:
            grid%xb%q(i,j,k)=grid%xb%q(i,j,k)/(1.0+grid%xb%q(i,j,k))
   
            ! Background qrn needed for radar radial velocity assmiilation:

            if (size(grid%moist,dim=4) >= 4) then
               grid%xb%qcw(i,j,k) = max(grid%moist(i,j,k,p_qc), 0.0)
               grid%xb%qrn(i,j,k) = max(grid%moist(i,j,k,p_qr), 0.0)
!rizvi For doing single obs test with radiance one need to ensure non-zero hydrometeor
!  For doing so uncomment following two cards below
!               grid%xb%qcw(i,j,k) = 0.0001                             
!               grid%xb%qrn(i,j,k) = 0.0001                             

               grid%xb%qt (i,j,k) = grid%xb%q(i,j,k) + grid%xb%qcw(i,j,k) + &
                  grid%xb%qrn(i,j,k)
            end if

            if (size(grid%moist,dim=4) >= 6) then
               grid%xb%qci(i,j,k) = max(grid%moist(i,j,k,p_qi), 0.0)
               grid%xb%qsn(i,j,k) = max(grid%moist(i,j,k,p_qs), 0.0)
!rizvi For doing single obs test with radiance one need to ensure non-zero hydrometeor
!  For doing so uncomment following two cards below
!               grid%xb%qci(i,j,k) = 0.0001                             
!               grid%xb%qsn(i,j,k) = 0.0001                             
            end if

            if (size(grid%moist,dim=4) >= 7) then
               grid%xb%qgr(i,j,k) = max(grid%moist(i,j,k,p_qg), 0.0)
            end if

            if ( config_flags%mp_physics == 3 ) then   ! WSM3-class scheme
               if ( grid%xb%t(i,j,k) <= t_kelvin ) then
                  grid%xb%qci(i,j,k) = grid%xb%qcw(i,j,k)
                  grid%xb%qcw(i,j,k) = 0.0
                  grid%xb%qsn(i,j,k) = grid%xb%qrn(i,j,k)
                  grid%xb%qrn(i,j,k) = 0.0
               end if
            end if

         end do
      end do

      do i=its,ite
         grid%xb%psac(i,j) = grid%mub(i,j)+grid%mu_2(i,j)
! To make the Psfc consistent with WRF (YRG, 04/06/2010):
         grid%xb%psfc(i,j) = grid%psfc(i,j)

         if (grid%xb%tgrn(i,j) < 100.0) &    
            grid%xb%tgrn(i,j) = grid%xb%t(i,j,kts)+ &
            0.0065*(grid%xb%h(i,j,kts)-grid%xb%hf(i,j,kts))
      end do
   end do

   end do
   !$OMP END PARALLEL DO

!   write (999,'("MABS=",e13.5)') sum(abs(grid%xb%psfc(:,:)-grid%psfc(:,:))) / &
!                                       float((ite-its+1)*(jte-jts+1))
      
   grid%xb%ztop = grid%xb%hf(its,jts,kte+1)

   if (print_detail_xb) then
      write(unit=stdout, fmt=*) ' '
      if (print_detail_xb) then
         write(unit=stdout, fmt='(/5a/)') &
            'lvl         h                 p                t'

         do k=kts,kte
            write(unit=stdout, fmt='(i3,8e20.12)') k, &
               grid%xb%h(its,jts,k), grid%xb%p(its,jts,k), grid%xb%t(its,jts,k)
         end do
      end if

      write(unit=stdout,fmt=*) ' '
      write(unit=stdout,fmt=*) 'grid%xb%u(its,jte,kte)=', grid%xb%u(its,jte,kte)
      write(unit=stdout,fmt=*) 'grid%xb%v(ite,jts,kte)=', grid%xb%v(ite,jts,kte)
      write(unit=stdout,fmt=*) 'grid%xb%w(its,jts,kte)=', grid%xb%w(its,jts,kte)
      write(unit=stdout,fmt=*) 'grid%xb%t(its,jts,kte)=', grid%xb%t(its,jts,kte)
      write(unit=stdout,fmt=*) 'grid%xb%p(its,jts,kte)=', grid%xb%p(its,jts,kte)
      write(unit=stdout,fmt=*) 'grid%xb%q(its,jts,kte)=', grid%xb%q(its,jts,kte)
      write(unit=stdout,fmt=*) 'grid%xb%h(its,jts,kte)=', grid%xb%h(its,jts,kte)
      write(unit=stdout,fmt=*) &
         'grid%xb%hf(its,jts,kte)=', grid%xb%hf(its,jts,kte)
      write(unit=stdout,fmt=*) &
         'grid%xb%map_factor(its,jts)=', grid%xb%map_factor(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%cori(its,jts)=', grid%xb%cori(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%tgrn(its,jts)=', grid%xb%tgrn(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%lat(its,jts)=', grid%xb%lat(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%lon(its,jts)=', grid%xb%lon(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%terr(its,jts)=', grid%xb%terr(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%snow(its,jts)=', grid%xb%snow(its,jts)
      write(unit=stdout,fmt=*) 'grid%xb%lanu(its,jts)=', grid%xb%lanu(its,jts)
      write(unit=stdout,fmt=*) &
         'grid%xb%landmask(its,jts)=', grid%xb%landmask(its,jts)
      write(unit=stdout,fmt=*) '(ite,jte)=', ite,jte                   
      write(unit=stdout,fmt=*) 'grid%xb%lat(ite,jte)=', grid%xb%lat(ite,jte)
      write(unit=stdout,fmt=*) 'grid%xb%lon(ite,jte)=', grid%xb%lon(ite,jte)
      write(unit=stdout,fmt=*) ' '
   end if

   !---------------------------------------------------------------------------
   ! [3.0] Calculate vertical inner product for use in vertical transform:
   !---------------------------------------------------------------------------
      
   if (vertical_ip == vertical_ip_sqrt_delta_p) then
      ! Vertical inner product is sqrt(Delta p):
      do k=kts,kte
         grid%xb % vertical_inner_product(its:ite,jts:jte,k) = &
            sqrt(grid%xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k))
      end do 
   else if (vertical_ip == vertical_ip_delta_p) then

      ! Vertical inner product is Delta p:
      do k=1,grid%xb%mkz
         grid % xb % vertical_inner_product(its:ite,jts:jte,k) = &
         grid % xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k)
      end do
   end if

   !---------------------------------------------------------------------------
   ! Roughness
   !---------------------------------------------------------------------------

   current_date = 'yyyy-mm-dd_hh:mm:ss'

   write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
      grid%start_year, '-', &
      grid%start_month, '-', &
      grid%start_day, '_', &
      grid%start_hour, ':', &
      grid%start_minute, ':', &
      grid%start_second

   !xbx % mminlu = 'USGS'
   xbx % mminlu = trim(grid%mminlu)

   call da_roughness_from_lanu(19, xbx % mminlu, current_date, &
      grid%xb % lanu, grid%xb % rough)

   !---------------------------------------------------------------------------
   ! Calculate 1/grid box areas:
   !---------------------------------------------------------------------------

   if (print_detail_xb) then
      write(unit=stdout, fmt='(/a, e24.16)') &
         'grid%xb % ds=', grid%xb % ds

      write(unit=stdout, fmt='(a, e24.16/)') &
           'grid%xb % map_factor(its,jts)=', grid%xb % map_factor(its,jts)
   end if

   !$OMP PARALLEL DO &
   !$OMP PRIVATE ( ij, i, j, tmpvar, height, message )
   do ij = 1 , grid%num_tiles

   do j=grid%j_start(ij),grid%j_end(ij)
      do i=its,ite
         if (grid%xb%ztop < grid%xb%hf(i,j,kte+1)) &
             grid%xb%ztop = grid%xb%hf(i,j,kte+1)

         tmpvar = grid%xb%ds / grid%xb%map_factor(i,j)

         grid%xb % grid_box_area(i,j) = tmpvar*tmpvar

         ! Calculate surface variable(wind, moisture, temperature)
         ! sfc variables: 10-m wind, and 2-m T, Q, at cross points

         height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
         if (height <= 0.0) then
            message(1) = "Negative height found"
            write (unit=message(2),FMT='(2I6,A,F10.2,A,F10.2)') &
               i,j,' ht = ',grid%xb%h(i,j,kts) ,' terr =  ',grid%xb%terr(i,j)
            call da_error(__FILE__,__LINE__, message(1:2))
         end if

#ifdef A2C
         uu = 0.5*(grid%xb%u(i,j,kts)+grid%xb%u(i+1,j,kts) )
         vv = 0.5*(grid%xb%v(i,j,kts)+grid%xb%v(i,j+1,kts) )
#endif
         call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
#ifdef A2C
            grid%xb%p(i,j,kts), grid%xb%t(i,j,kts), grid%xb%q(i,j,kts),uu,vv, &
#else
            grid%xb%p(i,j,kts), grid%xb%t(i,j,kts), grid%xb%q(i,j,kts), &
            grid%xb%u(i,j,kts), grid%xb%v(i,j,kts), &
#endif
            grid%xb%p(i,j,kts+1), grid%xb%t(i,j,kts+1), grid%xb%q(i,j,kts+1), &
            height,  grid%xb%rough(i,j),grid%xb%xland(i,j), grid%xb%ds, & ! Modified by Eric Chiang (NOV 2010) 
            grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), &
            grid%xb%q2(i,j), grid%xb%regime(i,j))
  
            grid%xb%u10(i,j) = grid%u10(i,j)         ! Add by Eric Chiang (NOV 2010)
            grid%xb%v10(i,j) = grid%v10(i,j)         ! Add by Eric Chiang (NOV 2010)
            grid%xb%t2(i,j) = grid%t2(i,j)           ! Add by Eric Chiang (NOV 2010)
            grid%xb%q2(i,j) = grid%q2(i,j)           ! Add by Eric Chiang (NOV 2010)

      end do
   end do

   end do
   !$OMP END PARALLEL DO

#ifdef VAR4D
   CALL set_physical_bc2d( grid%xb%grid_box_area, 't', config_flags, &
                           ids, ide, jds, jde,    &
                           ims, ime, jms, jme,    &
                           ips, ipe, jps, jpe,    &
                           its, ite, jts, jte    )
#endif

   !---------------------------------------------------------------------------
   ! Calculate saturation vapour pressure and relative humidity:
   !---------------------------------------------------------------------------

   !$OMP PARALLEL DO &
   !$OMP PRIVATE ( ij, k, j, i )
   do ij = 1 , grid%num_tiles
      do k=kts,kte
         do j=grid%j_start(ij),grid%j_end(ij)
            do i=its,ite
               call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
                  grid%xb % q(i,j,k), grid%xb %es(i,j,k), grid%xb %qs(i,j,k), &
                  grid%xb %rh(i,j,k))
            end do
         end do
      end do
   end do
   !$OMP END PARALLEL DO 

   !---------------------------------------------------------------------------
   ! Calculate dew point temperature:
   !---------------------------------------------------------------------------

   call da_trh_to_td (grid)

   if (print_detail_xb) then
      i=its; j=jts; k=kts

      write(unit=stdout, fmt=*) 'i,j,k=', i,j,k
      write(unit=stdout, fmt=*) 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
      write(unit=stdout, fmt=*) 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
      write(unit=stdout, fmt=*) 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
      write(unit=stdout, fmt=*) 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
      write(unit=stdout, fmt=*) ' '
   end if

   !---------------------------------------------------------------------------
   ! Sea level pressure and total precipitable water
   !---------------------------------------------------------------------------

   call da_wrf_tpq_2_slp (grid)

   ! WHY?
   ! do j = jts,jte
   !    do i = its,ite
   !       call da_tpq_to_slp(grid%xb%t(i,j,:), grid%xb%q(i,j,:), &
   !          grid%xb%p(i,j,:), grid%xb%terr(i,j), &
   !          grid%xb%psfc(i,j), grid%xb%slp(i,j), grid%xp)
   !    end do
   ! end do

   call da_integrat_dz(grid)

   !---------------------------------------------------------------------------
   ! Surface wind speed
   !---------------------------------------------------------------------------

   tmpvar = log(10.0/0.0001)

   !$OMP PARALLEL DO &
#ifndef A2C
   !$OMP PRIVATE (ij, i, j, height)
#else
   !$OMP PRIVATE (ij, i, j, height, uu, vv)
#endif
   do ij = 1, grid%num_tiles

   do j=grid%j_start(ij), grid%j_end(ij)
      do i=its,ite
         height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
         rgh_fac(i,j) = 1.0/log(height/0.0001)
#ifndef A2C
         grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,kts)*grid%xb%u(i,j,kts) &
                         + grid%xb%v(i,j,kts)*grid%xb%v(i,j,kts) + 1.0e-6) &
                    *tmpvar*rgh_fac(i,j)
#else
         uu = 0.5*(grid%xb%u(i,j,kts)+grid%xb%u(i+1,j,kts) )
         vv = 0.5*(grid%xb%v(i,j,kts)+grid%xb%v(i,j+1,kts) )
         grid%xb%speed(i,j) = sqrt(uu*uu + vv*vv + 1.0e-6)*tmpvar*rgh_fac(i,j)
#endif
      end do
   end do

   end do
   !$OMP END PARALLEL DO

   !---------------------------------------------------------------------------
   ! Brightness temperature SH Chen
   !---------------------------------------------------------------------------

   if (use_ssmitbobs)   &
      call da_transform_xtotb(grid)

   !---------------------------------------------------------------------------
   ! GPS Refractivity linked by Y.-R. Guo 05/28/2004
   !---------------------------------------------------------------------------

   call da_transform_xtogpsref(grid)

   !---------------------------------------------------------------------------
   ! Ground-based GPS ZTD must follow the GPS Refractivity calculation.
   !---------------------------------------------------------------------------

   ! WHY? For certain computation method, not current one.
   if (use_gpsztdobs) then
      call da_transform_xtoztd(grid)
      if (print_detail_xb) then
        i=its; j=jts
        write(unit=stdout, fmt=*) 'grid%xb % tpw(i,j)=', grid%xb % tpw(i,j)
        write(unit=stdout, fmt=*) 'grid%xb % ztd(i,j)=', grid%xb % ztd(i,j)
        write(unit=stdout, fmt=*) ' '
      end if
   end if

   !---------------------------------------------------------------------------
   ! Calculate means for later use in setting up background errors.
   !---------------------------------------------------------------------------

   ! WHY?
   ! if (.not. associated(xbx % latc_mean)) then
   allocate (xbx % latc_mean(jds:jde))
   if (trace_use) call da_trace("da_transfer_wrftoxb",&
      message="allocated xbx%latc_mean")
   ! end if

   size2d = (ide-ids+1)*(jde-jds+1)

   tmpvar = 1.0/real(size2d)

   ! Bitwitse-exact reduction preserves operation order of serial code for
   ! testing, at the cost of much-increased run-time.  Turn it off when not
   ! testing.  Thits will always be .false. for a serial or 1-process MPI run.
   if (test_dm_exact) then
      allocate(arrayglobal(ids:ide, jds:jde))
      call da_patch_to_global(grid,grid%xb%psac, arrayglobal)
      loc_psac_mean = tmpvar*sum(arrayglobal(ids:ide,jds:jde))
      deallocate(arrayglobal)
   else
      loc_latc_mean = 0.0
      loc_psac_mean = tmpvar*sum(grid%xb % psac(its:ite,jts:jte))
   end if

   tmpvar = 1.0/real(ide-ids+1)

   if (test_dm_exact) then
      allocate(arrayglobal(ids:ide, jds:jde))
      call da_patch_to_global(grid,grid%xb%lat, arrayglobal)
      do j=jds,jde
         loc_latc_mean(j) = tmpvar*sum(arrayglobal(ids:ide, j))
      end do
      deallocate(arrayglobal)
   else
      loc_latc_mean = 0.0
      do j=jts,jte
         loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(its:ite, j))
      end do
   end if

   if (test_dm_exact) then
      ! Broadcast result from monitor to other tasks.
      call wrf_dm_bcast_real(loc_psac_mean, 1)
      xbx % psac_mean = loc_psac_mean
      ! Broadcast result from monitor to other tasks.
      call wrf_dm_bcast_real(loc_latc_mean, (jde-jds+1))
      xbx % latc_mean = loc_latc_mean
   else
      xbx % psac_mean = wrf_dm_sum_real(loc_psac_mean)
      call wrf_dm_sum_reals(loc_latc_mean, xbx % latc_mean)
   end if

   if (print_detail_xb) then
      ! write(unit=stdout, fmt=*) 'loc_psac_mean  =', loc_psac_mean
      write(unit=stdout, fmt=*) 'xbx % psac_mean=', xbx % psac_mean

      ! write(unit=stdout, fmt=*) 'loc_latc_mean  =', loc_latc_mean(jts)
      write(unit=stdout, fmt=*) 'xbx % latc_mean=', xbx % latc_mean(jts)
   end if


   ! Fill the halo region for xb        

#ifdef DM_PARALLEL
#include "HALO_XB.inc"
#endif

   ! Calculate time step from one dimensional cloud model parameterization

   if (dt_cloud_model) then
      do j = jts, jte
         do i = its, ite
            call da_cloud_model (grid%xb%t(I,J,:),  grid%xb%p(I,J,:), &
               grid%xb%q(I,J,:), grid%xb%qcw(I,J,:), grid%xb%qrn(I,J,:), &
               grid%xb%h(I,J,:), grid%xb%hf(I,J,:), ddt, kts, kte)

            do k = kts, kte
               grid%xb%delt(i,j,k) = DDT(k)
            end do
         end do
      end do
   end if

   deallocate (xbx % latc_mean)

   if (trace_use) call da_trace_exit("da_transfer_wrftoxb")

end subroutine da_transfer_wrftoxb