subroutine da_sfc_wtq_lin(psfc, tg, ps, ts, qs, us, vs, regime,           &
   psfc_prime,tg_prime,ps_prime,ts_prime,qs_prime, &
   us_prime,vs_prime, hs,roughness,xland,dx,         & ! Modified by Eric Chiang (JULY 2010)
   u10_prime,v10_prime,t2_prime,q2_prime) 

   !---------------------------------------------------------------------------
   ! Purpose: Calculate the  10m wind, 2m temperature and moisture based on the
   ! similarity theory/
   !
   ! Reference:
   ! ---------
   !
   !  input Variables(basic state):
   ! 
   !   psfc, tg               : surface pressure and ground temperature
   !   ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
   !   regime                 : PBL regime
   !
   !  input Variables(pertubation):
   ! 
   !   psfc_prime, tg_prime   : Surface pressure and ground temperature
   !   ps_prime, ts_prime,    : Model variables at the lowest half sigma
   !   qs_prime, us_prime,    : level 
   !   vs_prime               : 
   !
   !  Constants:
   !
   !   hs                     : height at the lowest half sigma level
   !   roughness              : roughness
   !   xland                  : land-water-mask (=2 water, =1 land)
   !
   !  output Variables(pertubation):
   !
   !   u10_prime, v10_prime   : 10-m high observed wind components
   !   t2_prime , q2_prime    : 2-m high observed temperature and mixing ratio
   !
   !---------------------------------------------------------------------------
   !  
   !                      psim  : mechanical psi at lowlest sigma level
   !                      psim2 : mechanical psi at 2m 
   !                      psimz : mechanical psi at 10m 
   !
   !---------------------------------------------------------------------------

   implicit none

   real, intent (in)  :: regime
   real, intent (in)  :: ps , ts , qs , us, vs, psfc, tg
   real, intent (in)  :: ps_prime, ts_prime, qs_prime, us_prime, vs_prime, psfc_prime, tg_prime
   real, intent (in)  :: hs, roughness, xland
   real, intent (out) :: u10_prime, v10_prime, t2_prime, q2_prime

   ! Maximum number of iterations in computing psim, psih

   integer, parameter :: k_iteration = 10 
   !      integer, parameter :: k_iteration = 1

   ! h10 is the height of 10m where the wind observed
   ! h2  is the height of 2m where the temperature and 
   !        moisture observed.

   real, parameter :: h10 = 10.0, h2 = 2.0
   !
   ! Default roughness over the land

   real, parameter :: zint0 = 0.01 
   !
   ! Von Karman constant

   real, parameter :: k_kar = 0.4
   !
   ! Working variables

   real :: Vc2, Va2, V2 
   real :: rib, rcp, xx, yy, cc, Pi
   real :: psiw, psiz, mol, ust, hol, holz, hol2
   real :: psim, psimz, psim2, psih, psihz, psih2
   real :: psit, psit2, psiq, psiq2
   real :: gzsoz0, gz10oz0, gz2oz0
   real :: eg, qg, tvg, tvs
   real :: ths, thg, thvs, thvg, vsgd, vsgd2, dx
   real :: zq0, z0

   real :: Vc2_prime, Va2_prime, V2_prime
   real :: rib_prime, xx_prime, yy_prime
   real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
           hol_prime, holz_prime, hol2_prime
   real :: psim_prime, psimz_prime, psim2_prime, &
           psih_prime, psihz_prime, psih2_prime
   real :: psit_prime, psit2_prime, &
           psiq_prime, psiq2_prime
   real :: qg_prime, tvg_prime, tvs_prime
   real :: ths_prime, thg_prime, thvs_prime, thvg_prime 

   real, parameter :: ka = 2.4E-5

   integer :: iregime

   if (trace_use) call da_trace_entry("da_sfc_wtq_lin")

   rcp = gas_constant/cp

   ! 1 Compute the roughness length based upon season and land use 
   ! =====================================

   ! 1.1 Define the rouhness length
   !     -----------------

   z0 = roughness

   if (z0 < 0.0001) z0 = 0.0001

   ! 1.2 Define the rouhgness length for moisture
   !     -----------------

   if (xland .ge. 1.5) then
      zq0 = z0
    else
      zq0 =  zint0
   end if

   ! 1.3 Define the some constant variable for psi
   !     -----------------

   gzsoz0 = log(hs/z0)

   gz10oz0 = log(h10/z0)

   gz2oz0 = log(h2/z0)


   ! 2. Calculate the virtual temperature
   ! =====================================

   ! 2.1 Compute Virtual temperature on the lowest half sigma level
   !     ---------------------------------------------------------

   tvs_prime  = ts_prime * (1.0 + 0.608 * qs) + 0.608 * ts * qs_prime
   tvs  = ts * (1.0 + 0.608 * qs)

   ! 2.2 Compute the ground saturated mixing ratio and the ground virtual 
   !     temperature
   !     ----------------------------------------------------------------

   call da_tp_to_qs(tg, psfc, eg, qg)
   call da_tp_to_qs_lin1(tg, psfc, eg, tg_prime, psfc_prime, qg_prime)

   qg_prime = qg_prime * qg

   tvg_prime  = tg_prime * (1.0 + 0.608 * qg) + 0.608 * tg * qg_prime
   tvg  = tg * (1.0 + 0.608 * qg)

   ! 3.  Compute the potential temperature and virtual potential temperature
   ! =======================================================================

   ! 3.1 Potential temperature on the lowest half sigma level
   !     ----------------------------------------------------

   Pi = (100000.0 / ps) ** rcp
   ths_prime  = (ts_prime - ps_prime * rcp * ts/ps) * Pi 
   ths  = ts * Pi

   ! 3.2 Virtual potential temperature on the lowest half sigma level
   !     ------------------------------------------------------------

   thvs_prime  = (tvs_prime - ps_prime * rcp * tvs/ps) * Pi 
   thvs = tvs * Pi

   ! 3.3 Potential temperature at the ground
   !     -----------------------------------

   Pi = (100000.0 / psfc) ** rcp
   thg_prime  = (tg_prime - psfc_prime * rcp * tg/psfc) * Pi 
   thg  = tg * Pi

   ! 3.4 Virtual potential temperature at ground
   !     ---------------------------------------

   thvg_prime  = (tvg_prime - psfc_prime * rcp * tvg/psfc) * Pi
   thvg = tvg * Pi

   ! 4.  BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
   ! =================================================

   ! 4.1 Velocity
   !     --------
   
   ! Wind speed:

   Va2_prime =  2.0*us*us_prime + 2.0*vs*vs_prime
   Va2 =   us*us + vs*vs
    
   ! Convective velocity:

   if (thvg >= thvs) then
      Vc2_prime = 4.0 * (thvg_prime - thvs_prime)
      Vc2 = 4.0 * (thvg - thvs)
   else
      Vc2_prime = 0.0
      Vc2 = 0.0
   end if
 
   ! Calculate Mahrt and Sun low-res correction                    ! Add by Eric Chiang ( July 2010 )
   ! dx is a constant, so vsgd is also a constant, 
   ! the perturnations of vsgd and vsgd2 are zero. (YRG, 09/15/2011)
                                     
   vsgd = 0.32 * (max(dx/5000.-1.,0.))**0.33                       ! Add by Eric Chiang ( July 2010 )
   vsgd2 = vsgd * vsgd                                             ! Add by Eric Chiang ( July 2010 )

   ! V2_prime should be computed before used below. (YRG, 09/15/2011)    
   V2_prime = Va2_prime + Vc2_prime
   V2  = Va2 + Vc2 + vsgd2                                         ! Modified by Eric Chiang ( July 2010 )

   ! 4.2 Bulk richardson number
   !     ----------------------

   Pi = gravity * hs / (ths*V2)
   rib_prime = (thvs_prime - thvg_prime   &
              - (thvs-thvg)/V2  * V2_prime &
              - (thvs-thvg)/ths * ths_prime) * Pi 
   rib = (thvs - thvg) * Pi

   ! 5.  CALCULATE PSI BASED UPON REGIME
   ! =======================================

   iregime = int(regime)

   select case (iregime) 

   ! 5.1 Stable conditions (REGIME 1)
   !     ---------------------------

   case (1);

      psim_prime  = 0.0
      psimz_prime = 0.0
      psim2_prime = 0.0
      psim  = -10.0*gzsoz0
      psimz = -10.0*gz10oz0
      psim2 = -10.0*gz2oz0
      psim  = max(psim,-10.0)
      psimz = max(psimz,-10.0)
      psim2 = max(psim2,-10.0)

      psih_prime  = psim_prime
      psihz_prime = psimz_prime
      psih2_prime = psim2_prime
      psih  = psim
      psihz = psimz
      psih2 = psim2

   ! 5.2 Mechanically driven turbulence (REGIME 2)
   !     ------------------------------------------

   case (2);

      Pi =  - 1.0 / ((1.1 - 5.0*rib)*(1.1 - 5.0*rib))
      psim_prime  = 5.5 * gzsoz0  * rib_prime * Pi 
      psimz_prime = 5.5 * gz10oz0 * rib_prime * Pi
      psim2_prime = 5.5 * gz2oz0  * rib_prime * Pi

      Pi =  (-5.0 * rib) / (1.1 - 5.0*rib)
      psim  = gzsoz0  * Pi
      psimz = gz10oz0 * Pi
      psim2 = gz2oz0  * Pi

      if (psim >= -10.0) then
         psim = psim
         psim_prime = psim_prime
      else
         psim = -10.0
         psim_prime = 0.0
      end if
      if (psimz >= -10.0) then
         psimz = psimz
         psimz_prime = psimz_prime
      else
         psimz = -10.0
         psimz_prime = 0.0
      end if
      if (psim2 >= -10.0) then
         psim2 = psim2
         psim2_prime = psim2_prime
      else
         psim2 = -10.0
         psim2_prime = 0.0
      end if

      psih_prime  = psim_prime
      psihz_prime = psimz_prime
      psih2_prime = psim2_prime
      psih = psim
      psihz = psimz
      psih2 = psim2

   ! 5.3 Unstable Forced convection (REGIME 3)
   !     -------------------------------------

   case (3);

      psim_prime  = 0.0
      psimz_prime = 0.0
      psim2_prime = 0.0

      psim  = 0.0
      psimz = 0.0
      psim2 = 0.0

      psih_prime  = psim_prime
      psihz_prime = psimz_prime
      psih2_prime = psim2_prime

      psih  = psim
      psihz = psimz
      psih2 = psim2


      ! 5.4 Free convection (REGIME 4)
      !     --------------------------

   case (4);

      ! Calculate psi m and pshi h using iteration method

      psim_prime = 0.0
      psih_prime = 0.0
      psim = 0.0
      psih = 0.0
      cc = 2.0 * atan(1.0)

      ! do k = 1 , k_iteration

      ! 5.4.1  Calculate   ust, m/L (mol), h/L (hol)
      !        --------------------------

      ! Friction speed

      ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
      ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust

      ! Heat fux factor

      mol = k_kar * (ths - thg)/(gzsoz0 - psih)
      mol_prime = ((ths_prime - thg_prime) /(ths - thg) + &
                     psih_prime /(gzsoz0 - psih)) * mol

      ! Ratio of PBL height to Monin-Obukhov length

      if (ust .LT. 0.01) then
         hol_prime = rib_prime * gzsoz0
         hol = rib * gzsoz0
      else
         hol = k_kar * gravity * hs * mol / (ths * ust * ust)
         hol_prime = (mol_prime / mol - ths_prime / ths &
                        - 2.0* ust_prime / ust) * hol
      end if

      ! 5.4.2  Calculate n, nz, R, Rz
      !        --------------------------

      if (hol >= 0.0) then
         hol_prime = 0.0
         hol = 0.0
      else
         hol_prime = hol_prime
         hol = hol
      end if
      if (hol >= -10.0) then
         hol_prime = hol_prime
         hol = hol
      else
         hol_prime = 0.0
         hol = -10.0
      end if

      holz_prime = (h10 / hs) * hol_prime
      holz = (h10 / hs) * hol
      if (holz >= 0.0) then
         holz_prime = 0.0
         holz = 0.0
      else
         holz_prime = holz_prime
         holz = holz
      end if
      if (holz >= -10.0) then
         holz_prime = holz_prime
         holz = holz
      else
         holz_prime = 0.0
         holz = -10.0
      end if

      hol2_prime = (h2 / hs) * hol_prime
      hol2 = (h2 / hs) * hol
      if (hol2 >= 0.0) then
         hol2_prime = 0.0
         hol2 = 0.0
      else
         hol2_prime = hol2_prime
         hol2 = hol2
      end if
      if (hol2 >= -10.0) then
         hol2_prime = hol2_prime
         hol2 = hol2
      else
         hol2_prime = 0.0
         hol2 = -10.0
      end if

      ! 5.4.3 Calculate Psim & psih
      !        --------------------------

      ! Using the continuous function:
      xx_prime = -4.0* hol_prime /((1.0 - 16.0 * hol) ** 0.75)
      xx = (1.0 - 16.0 * hol) ** 0.25
      yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
      yy = log((1.0+xx*xx)/2.0)
      psim_prime = 2 * xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime 
      psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
      psih_prime = 2.0 * yy_prime
      psih = 2.0 * yy

      ! Using the continuous function:
      xx_prime = -4.0* holz_prime /((1.0 - 16.0 * holz) ** 0.75)
      xx = (1.0 - 16.0 * holz) ** 0.25
      yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
      yy = log((1.0+xx*xx)/2.0)
      psimz_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
      psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
      psihz_prime = 2.0 * yy_prime
      psihz = 2.0 * yy

      ! Using the continuous function:
      xx_prime = -4.0* hol2_prime /((1.0 - 16.0 * hol2) ** 0.75)
      xx = (1.0 - 16.0 * hol2) ** 0.25
      yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
      yy = log((1.0+xx*xx)/2.0)
      psim2_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
      psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
      psih2_prime = 2.0 * yy_prime
      psih2 = 2.0 * yy

      ! end do 

      ! 5.4.4 Define the limit value for psim & psih
      !        --------------------------

      if (psim <= 0.9*gzsoz0) then
         psim_prime = psim_prime
         psim = psim
      else
         psim = 0.9*gzsoz0
         psim_prime = 0.0
      end if
      if (psimz <= 0.9*gz10oz0) then
         psimz_prime = psimz_prime
         psimz = psimz
      else
         psimz_prime = 0.0
         psimz = 0.9*gz10oz0
      end if
      if (psim2 <= 0.9*gz2oz0) then
         psim2_prime = psim2_prime
         psim2 = psim2
      else
         psim2_prime = 0.0
         psim2 = 0.9*gz2oz0
      end if
      if (psih <= 0.9*gzsoz0) then
         psih_prime = psih_prime
         psih = psih
      else
         psih_prime = 0.0
         psih = 0.9*gzsoz0
      end if
      if (psihz <= 0.9*gz10oz0) then
         psihz_prime = psihz_prime
         psihz = psihz
      else
         psihz_prime = 0.0
         psihz = 0.9*gz10oz0
      end if
      if (psih2 <= 0.9*gz2oz0) then
         psih2_prime = psih2_prime
         psih2 = psih2
      else
         psih2_prime = 0.0
         psih2 = 0.9*gz2oz0
       end if

    case default;
       write(unit=message(1),fmt='(A,I2,A)') &
          "Regime=",iregime," is invalid."
       call da_error(__FILE__,__LINE__,message(1:1))

   end select

   ! 6.  CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
   ! =======================================

   psiw_prime = - psim_prime
   psiw = gzsoz0 - psim
   psiz_prime = - psimz_prime
   psiz = gz10oz0 - psimz
   psit_prime = - psih_prime
   psit = gzsoz0 - psih
   psit2_prime = - psih2_prime
   psit2 = gz2oz0 - psih2

   ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
   ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust

   psiq2_prime = k_kar*hs/(ka*(k_kar*ust*hs/ka + hs / zq0))*ust_prime
   psiq_prime  = psiq2_prime - psih_prime
   psiq2_prime = psiq2_prime - psih2_prime

   psiq  = log(k_kar*ust*hs/ka + hs / zq0) - psih
   psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2

   ! 7.  CALCULATE THE PERTURBATIONS for 10M WinD, 2M TEMPERATURE AND MOISTURE
   ! =======================================

   Pi = psiz / psiw
   u10_prime= (us_prime + us/psiz * psiz_prime - us/psiw * psiw_prime) * Pi
   v10_prime= (vs_prime + vs/psiz * psiz_prime - vs/psiw * psiw_prime) * Pi

   t2_prime = ((1.0-psit2/psit) * thg_prime + (ths_prime + &
                           (ths - thg)/psit2 * psit2_prime - &
                           (ths - thg)/psit  * psit_prime) * psit2/psit &
             + rcp*(thg + (ths - thg)*psit2/psit)/psfc * psfc_prime) &
             * (psfc/100000.0)**rcp

   q2_prime = (1.0-psiq2/psiq) * qg_prime + psiq2/psiq * qs_prime + &
              (qs -qg)*(psiq2/psiq) * (psiq2_prime/psiq2 - psiq_prime/psiq)

   if (trace_use) call da_trace_exit("da_sfc_wtq_lin")

end subroutine da_sfc_wtq_lin