var/da/da_physics/da_sfc_wtq.inc

   1 subroutine da_sfc_wtq (psfc, tg, ps, ts, qs, us, vs, &
   2    hs, roughness, xland, dx, u10, v10, t2, q2, regime, &
   3    has_lsm, regime_wrf, qsfc_wrf, znt_wrf, ust_wrf, mol_wrf, hfx, qfx, pblh)
   4
   5    !---------------------------------------------------------------------------
   6    ! Purpose: Calculate the  10m wind, 2m temperature and moisture based on the
   7    ! similarity theory/
   8    !
   9    !  The unit for pressure   : psfc, ps          is Pa.
  10    !  The unit for temperature: tg, ts, t2        is K.
  11    !  The unit for moisture   : qs, q2            is kg/kg.
  12    !  The unit for wind       : us, vs, u10, v10  is m/s.
  13    !  The unit for height     : hs, roughness     is m.
  14    !  xland and regime are dimensionless.
  15    !
  16    ! History: Nov 2010 - improve calculation consistency with WRF model (Eric Chiang)
  17    !          Jul 2015 - further improvement on consistency
  18    !
  19    ! Reference:
  20    ! ---------
  21    !
  22    !  input Variables:
  23    !
  24    !   psfc, tg               : surface pressure and ground temperature
  25    !   ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
  26    !   dx  (m)                : horizontal resolution
  27    !
  28    !
  29    !  Constants:
  30    !
  31    !   hs                     : height at the lowest half sigma level
  32    !   roughness              : roughness
  33    !   xland                  : land-water-mask (=2 water, =1 land)
  34    !
  35    !  output Variables:
  36    !
  37    !   regime                 : PBL regime
  38    !   u10, v10               : 10-m high observed wind components
  39    !   t2 , q2                : 2-m high observed temperature and mixing ratio
  40    !
  41    !---------------------------------------------------------------------------
  42    !
  43    !                      psim  : mechanical psi at lowlest sigma level
  44    !                      psim2 : mechanical psi at 2m
  45    !                      psimz : mechanical psi at 10m
  46    !
  47    !---------------------------------------------------------------------------
  48
  49    implicit none
  50
  51    real, intent (in)  :: ps , ts , qs , us, vs
  52    real, intent (in)  :: psfc, tg
  53    real, intent (in)  :: hs, roughness, xland
  54    real, intent (out) :: regime
  55    real, intent (out) :: u10, v10, t2, q2
  56    logical, intent(in), optional :: has_lsm
  57    real,    intent(in), optional :: regime_wrf, qsfc_wrf, znt_wrf, ust_wrf, mol_wrf
  58    real,    intent(in), optional :: hfx, qfx, pblh
  59
  60    logical :: use_table = .true.
  61    logical :: use_ust_wrf = .false.
  62    logical :: vconv_wrf
  63    integer :: nn, nz, n2
  64    real    :: rr, rz, r2
  65    real    :: cqs2, chs2, rho, rhox, fluxc, visc, restar, z0t, z0q
  66
  67    ! h10 is the height of 10m where the wind observed
  68    ! h2  is the height of 2m where the temperature and
  69    !        moisture observed.
  70
  71    real, parameter :: h10 = 10.0, h2 = 2.0
  72
  73    ! Default roughness over the land
  74
  75    real, parameter :: zint0 = 0.01
  76
  77    ! Von Karman constant
  78
  79    real, parameter :: k_kar = 0.4
  80
  81    ! Working variables
  82
  83    real :: Vc2, Va2, V2, vc, wspd
  84    real :: rib, rcp, xx, yy, cc
  85    real :: psiw, psiz, mol, ust, hol, holz, hol2
  86    real :: psim, psimz, psim2, psih, psihz, psih2
  87    real :: psit, psit2, psiq, psiq2
  88    real :: gzsoz0, gz10oz0, gz2oz0
  89    real :: eg, qg, tvg, tvs, tvs2
  90    real :: ths, thg, thvs, thvg, thvs2, vsgd, vsgd2, dx
  91    real :: zq0, z0
  92
  93    real, parameter :: ka = 2.4E-5
  94
  95    if (trace_use_dull) call da_trace_entry("da_sfc_wtq")
  96
  97    rcp = gas_constant/cp
  98
  99    ! 1 Compute the roughness length based upon season and land use
 100
 101    ! 1.1 Define the roughness length
 102
 103    z0 = roughness
 104
 105    if (z0 < 0.0001) z0 = 0.0001
 106
 107    if ( present(znt_wrf) ) then
 108       if ( znt_wrf > 0.0 ) then
 109          z0 = znt_wrf
 110       end if
 111    end if
 112
 113    ! 1.2 Define the rouhgness length for moisture
 114
 115    if (xland .ge. 1.5) then
 116       zq0 = z0
 117    else
 118       zq0 =  zint0
 119    end if
 120
 121    ! 1.3 Define the some constant variable for psi
 122
 123    gzsoz0 = log(hs/z0)
 124
 125    gz10oz0 = log(h10/z0)
 126
 127    gz2oz0 = log(h2/z0)
 128
 129
 130    ! 2. Calculate the virtual temperature
 131
 132    ! 2.1 Compute Virtual temperature on the lowest half sigma level
 133
 134    tvs  = ts * (1.0 + 0.608 * qs)
 135
 136    ! 2.2 Convert ground virtual temperature assuming it's saturated
 137
 138    call da_tp_to_qs(tg, psfc, eg, qg) !output qg is specific humidity
 139    qg = qg*(1.0-qg) !hcl convert to mixing ratio
 140    if ( present(qsfc_wrf) ) then
 141       if ( qsfc_wrf > 0.0 ) then
 142          qg = qsfc_wrf
 143       end if
 144    endif
 145
 146    tvg  = tg * (1.0 + 0.608 * qg)
 147
 148    ! 3.  Compute the potential temperature
 149
 150    ! 3.1 Potential temperature on the lowest half sigma level
 151
 152    ths  = ts * (1000.0 / (ps/100.0)) ** rcp
 153
 154    ! 3.2 Potential temperature at the ground
 155
 156    thg  = tg * (1000.0 / (psfc/100.0)) ** rcp
 157
 158    ! 4. Virtual potential temperature
 159
 160    ! 4.1 Virtual potential temperature on the lowest half sigma level
 161
 162    thvs = tvs * (1000.0 / (ps/100.0)) ** rcp
 163
 164    ! 4.2 Virtual potential temperature at ground
 165
 166    thvg = tvg * (1000.0 / (psfc/100.0)) ** rcp
 167
 168
 169    ! 5.  BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
 170
 171    ! 5.1 Velocity
 172
 173    !     Wind speed:
 174
 175    Va2 =   us*us + vs*vs
 176    !
 177    !     Convective velocity:
 178
 179    vconv_wrf = .false.
 180    if ( present(hfx) .and. present(qfx) .and. present(pblh) ) then
 181       ! calculating vconv over land following wrf method
 182       if ( pblh > 0.0 ) then
 183          vconv_wrf = .true.
 184       end if
 185    end if
 186
 187    if (thvg >= thvs) then
 188       ! prior to V3.7, Vc2 = 4.0 * (thvg - thvs)
 189       Vc2 = thvg - thvs
 190    else
 191       Vc2 = 0.0
 192    end if
 193    if ( xland < 1.5 ) then !land
 194       if ( vconv_wrf ) then
 195          ! following the calculation as in module_sf_sfclay.F
 196          rhox = psfc/(gas_constant*tvg)
 197          fluxc = max(hfx/rhox/cp+0.608*tvg*qfx/rhox, 0.0)
 198          vc = (gravity/tg*pblh*fluxc)**0.33
 199          vc2 = vc*vc
 200       end if
 201    end if
 202
 203    ! Calculate Mahrt and Sun low-res correction         ! Add by Eric Chiang ( July 2010 )
 204    vsgd = 0.32 * (max(dx/5000.-1.,0.))**0.33            ! Add by Eric Chiang ( July 2010 )
 205    vsgd2 = vsgd * vsgd                                  ! Add by Eric Chiang ( July 2010 )
 206
 207    V2  = Va2 + Vc2 + vsgd2                              ! Add by Eric Chiang ( July 2010 )
 208    wspd = sqrt(v2)
 209    wspd = max(wspd,0.1)
 210    v2 = wspd*wspd
 211
 212    ! 5.2 Bulk richardson number
 213
 214    rib = (gravity * hs / ths) * (thvs - thvg) / V2
 215
 216    ! if previously unstable, do not let into regime 1 and 2
 217    if ( present(mol_wrf) ) then
 218       if ( mol_wrf < 0.0 ) rib = min(rib, 0.0)
 219    end if
 220
 221    !  Calculate   ust, m/L (mol), h/L (hol)
 222
 223    psim = 0.0
 224    psih = 0.0
 225
 226    ! Friction speed
 227
 228    if ( present(ust_wrf) ) then
 229       if ( ust_wrf > 0.0 ) then
 230          use_ust_wrf = .true.
 231          ust = ust_wrf
 232       end if
 233    end if
 234    if ( .not. use_ust_wrf ) then
 235       !ust = 0.0001  !init value as in phys/module_physics_init.F
 236       ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
 237    end if
 238
 239    ! Heat flux factor
 240
 241    if ( present(mol_wrf) ) then
 242       mol = mol_wrf
 243    else
 244       mol = k_kar * (ths - thg)/(gzsoz0 - psih)
 245       !mol = 0.0
 246    end if
 247
 248    ! set regimes based on rib
 249    if (rib .GE. 0.2) then
 250       ! Stable conditions (REGIME 1)
 251       regime = 1.1
 252    else if ((rib .LT. 0.2) .AND. (rib .GT. 0.0)) then
 253       ! Mechanically driven turbulence (REGIME 2)
 254       regime = 2.1
 255    else if (rib .EQ. 0.0) then
 256       ! Unstable Forced convection (REGIME 3)
 257       regime = 3.1
 258    else
 259       ! Free convection (REGIME 4)
 260       regime = 4.1
 261    end if
 262
 263    if ( present(regime_wrf) ) then
 264       if ( regime_wrf > 0.0 ) then
 265          regime = regime_wrf
 266       end if
 267    end if
 268
 269    ! 6.  CALCULATE PSI BASED UPON REGIME
 270
 271    !if (rib .GE. 0.2) then
 272    if ( nint(regime) == 1 ) then
 273       ! 6.1 Stable conditions (REGIME 1)
 274       !     ---------------------------
 275       regime = 1.1
 276       psim = -10.0*gzsoz0
 277       psim = max(psim,-10.0)
 278       psimz = h10/hs * psim
 279       psimz = max(psimz,-10.0)
 280       psim2 = h2 /hs * psim
 281       psim2 = max(psim2,-10.0)
 282       psih = psim
 283       psihz = psimz
 284       psih2 = psim2
 285
 286    !else if ((rib .LT. 0.2) .AND. (rib .GT. 0.0)) then
 287    else if ( nint(regime) == 2 ) then
 288
 289       ! 6.2 Mechanically driven turbulence (REGIME 2)
 290
 291       regime = 2.1
 292       psim = (-5.0 * rib) * gzsoz0 / (1.1 - 5.0*rib)
 293       psim = max(psim,-10.0)
 294       psimz = h10/hs * psim
 295       psimz = max(psimz,-10.0)
 296       psim2 = h2 /hs * psim
 297       psim2 = max(psim2,-10.0)
 298       psih = psim
 299       psihz = psimz
 300       psih2 = psim2
 301
 302    !else if (rib .EQ. 0.0) then
 303    else if ( nint(regime) == 3 ) then
 304       ! 6.3 Unstable Forced convection (REGIME 3)
 305
 306       regime = 3.1
 307       psim = 0.0
 308       psimz = 0.0
 309       psim2 = 0.0
 310       psih = psim
 311       psihz = psimz
 312       psih2 = psim2
 313
 314    else
 315       ! 6.4 Free convection (REGIME 4)
 316       regime = 4.1
 317
 318       cc = 2.0 * atan(1.0)
 319
 320       ! Ratio of PBL height to Monin-Obukhov length
 321
 322       if (ust .LT. 0.01) then
 323          hol = rib * gzsoz0
 324       else
 325          hol = k_kar * gravity * hs * mol / (ths * ust * ust)
 326       end if
 327
 328       ! 6.4.2  Calculate n, nz, R, Rz
 329
 330       holz = (h10 / hs) * hol
 331       hol2 = (h2 / hs) * hol
 332
 333       hol = min(hol,0.0)
 334       hol = max(hol,-9.9999)
 335
 336       holz = min(holz,0.0)
 337       holz = max(holz,-9.9999)
 338
 339       hol2 = min(hol2,0.0)
 340       hol2 = max(hol2,-9.9999)
 341
 342       ! 6.4.3 Calculate Psim & psih
 343
 344       if ( use_table ) then
 345          ! Using the look-up table:
 346          nn = int(-100.0 * hol)
 347          rr = (-100.0 * hol) - nn
 348          psim = psimtb(nn) + rr * (psimtb(nn+1) - psimtb(nn))
 349          psih = psihtb(nn) + rr * (psihtb(nn+1) - psihtb(nn))
 350       else
 351          ! Using the continuous function:
 352          xx = (1.0 - 16.0 * hol) ** 0.25
 353          yy = log((1.0+xx*xx)/2.0)
 354          psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
 355          psih = 2.0 * yy
 356       end if
 357
 358       if ( use_table ) then
 359          ! Using the look-up table:
 360          nz = int(-100.0 * holz)
 361          rz = (-100.0 * holz) - nz
 362          psimz = psimtb(nz) + rz * (psimtb(nz+1) - psimtb(nz))
 363          psihz = psihtb(nz) + rz * (psihtb(nz+1) - psihtb(nz))
 364       else
 365          ! Using the continuous function:
 366          xx = (1.0 - 16.0 * holz) ** 0.25
 367          yy = log((1.0+xx*xx)/2.0)
 368          psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
 369          psihz = 2.0 * yy
 370       end if
 371
 372       if ( use_table ) then
 373          ! Using the look-up table:
 374          n2 = int(-100.0 * hol2)
 375          r2 = (-100.0 * hol2) - n2
 376          psim2 = psimtb(n2) + r2 * (psimtb(n2+1) - psimtb(n2))
 377          psih2 = psihtb(n2) + r2 * (psihtb(n2+1) - psihtb(n2))
 378       else
 379          ! Using the continuous function:
 380          xx = (1.0 - 16.0 * hol2) ** 0.25
 381          yy = log((1.0+xx*xx)/2.0)
 382          psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
 383          psih2 = 2.0 * yy
 384       end if
 385
 386       ! 6.4.4 Define the limit value for psim & psih
 387
 388       psim = min(psim,0.9*gzsoz0)
 389       psimz = min(psimz,0.9*gz10oz0)
 390       psim2 = min(psim2,0.9*gz2oz0)
 391       psih = min(psih,0.9*gzsoz0)
 392       psihz = min(psihz,0.9*gz10oz0)
 393       psih2 = min(psih2,0.9*gz2oz0)
 394    end if  ! Regime
 395
 396    ! 7.  Calculate psi for wind, temperature and moisture
 397
 398    psiw = gzsoz0 - psim
 399    psiz = gz10oz0 - psimz
 400    psit = max(gzsoz0-psih, 2.0)
 401    psit2 = gz2oz0 - psih2
 402
 403    if ( .not. use_ust_wrf ) then
 404       ! re-calculate ust since psim is now available
 405       ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
 406    end if
 407
 408    psiq  = log(k_kar*ust*hs/ka + hs / zq0) - psih
 409    psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
 410
 411    !V3.7, as in module_sf_sfclay.F
 412    if ( xland >= 1.5 ) then !water
 413       visc = (1.32+0.009*(ts-273.15))*1.e-5
 414       restar = ust*z0/visc
 415       z0t = (5.5e-5)*(restar**(-0.60))
 416       z0t = min(z0t,1.0e-4)
 417       z0t = max(z0t,2.0e-9)
 418       z0q = z0t
 419       psiq  = max(log((hs+z0q)/z0q)-psih,  2.)
 420       psit  = max(log((hs+z0t)/z0t)-psih,  2.)
 421       psiq2 = max(log((2.+z0q)/z0q)-psih2, 2.)
 422       psit2 = max(log((2.+z0t)/z0t)-psih2, 2.)
 423    end if
 424
 425    ! 8.  Calculate 10m wind, 2m temperature and moisture
 426
 427    u10 = us * psiz / psiw
 428    v10 = vs * psiz / psiw
 429    t2 = (thg + (ths - thg)*psit2/psit)*((psfc/100.0)/1000.0)**rcp
 430    q2 = qg + (qs - qg)*psiq2/psiq
 431
 432    if ( present(has_lsm) ) then
 433       if ( has_lsm ) then
 434          !cqs2: 2m surface exchange coefficient for moisture
 435          !chs2: 2m surface exchange coefficient for heat
 436          cqs2 = ust*k_kar/psiq2
 437          if (xland .ge. 1.5) then
 438             !water
 439             chs2 = ust*k_kar/psit2
 440          else
 441             !land
 442             chs2 = cqs2 !as in subroutine lsm in phys/module_sf_noahdrv.F
 443          end if
 444
 445          !re-calculate T2/Q2 as in module_sf_sfcdiags.F
 446          rho  = psfc/(gas_constant*tg)
 447          if ( cqs2 < 1.e-5 ) then
 448             q2 = qg
 449          else
 450             if ( present(qfx) ) then
 451                q2 = qg - qfx/(rho*cqs2)
 452             end if
 453          end if
 454          if ( chs2 < 1.e-5 ) then
 455             t2 = tg
 456          else
 457             if ( present(hfx) ) then
 458                t2 = tg - hfx/(rho*cp*chs2)
 459             end if
 460          end if
 461       end if
 462    end if
 463
 464    if (trace_use_dull) call da_trace_exit("da_sfc_wtq")
 465
 466 end subroutine da_sfc_wtq
 467
 468