var/da/da_tools/da_mo_correction.inc

   1 function da_mo_correction (ho, po, to, qo, &
   2                         hm, pm, tm, qm, um ,vm, &
   3                         roughness)   RESULT (correc)
   4
   5    !----------------------------------------------------------------------------
   6    ! Purpose: Compute the correction factor to convert surface wind into 40m
   7    ! wind using similarity laws.
   8    !
   9    ! Reference:
  10    ! ---------
  11    !  A description of the fifth generation Penn State/NCAR Mesoscale Model
  12    !  Grell et al. 1994, page 29-30 and 80-83.
  13    !----------------------------------------------------------------------------
  14
  15    implicit none
  16
  17    real, intent (in)                :: ho, po, to, qo
  18    real, intent (in)                :: hm, pm, tm, qm, um, vm
  19    real, intent (in)                :: roughness
  20
  21    real                             :: correc, winfac
  22
  23    real :: thm , tho, tvm, tvo, thvm, thvo, rcp
  24    real :: za, Vc2, Va2, V2
  25    ! FIX? real :: hdif, rib, rll, hll, zint
  26    real :: hdif, rib, hll, zint
  27
  28    ! Height difference (in m) above wich correction is applied
  29
  30    real, parameter :: hmax = 10.0, hs_max = 40.0
  31
  32    ! Default Roughness length im m (for land, set to 0 if on sea)
  33
  34    real, parameter :: zint0 = 0.2
  35
  36    rcp = gas_constant/cp
  37
  38    ! 0.0  initialize correction factor to 1
  39    ! =====================================
  40
  41    correc = 1.0
  42    winfac = 1.0
  43
  44    ! 1.  height difference and roughness length
  45    ! ==========================================
  46
  47    ! 1.1 Height difference
  48    !     -----------------
  49
  50    hdif = hm - ho
  51
  52    ! 1.2 No correction if height difference is less than hmax.
  53    !     -----------------------------------------------------
  54
  55    if (hdif <= hmax) then
  56       return
  57    end if
  58
  59    ! too many
  60    ! if (trace_use) call da_trace_entry("da_mo_correction")
  61
  62    ! 1.3 Compute the roughness length based upon season and land use
  63    !     -----------------------------------------------------------
  64
  65    zint = roughness
  66
  67    if (zint < 0.0001) zint = 0.0001
  68
  69    ! 2.  potential temperature at model surface and observation location
  70    ! ===================================================================
  71
  72    ! 2.1 potential temperature on model surface
  73    !     --------------------------------------
  74
  75    thm  = tm * (1000.0 / (pm/100.0)) ** rcp
  76
  77    ! 2.2 potential temperature at observation location
  78    !     ---------------------------------------------
  79
  80    tho  = to * (1000.0 / (po/100.0)) ** rcp
  81
  82    ! 3.  virtual temperature at model surface and observation location
  83    ! ===================================================================
  84
  85    ! 3.1 Virtual temperature on model surface
  86    !     -------------------------------------
  87
  88    tvm  = tm * (1.0 + 0.608 * qm)
  89
  90    ! 3.2 Virtual temperature at observation location
  91    !     -------------------------------------------
  92
  93       tvo  = to * (1.0 + 0.608 * qo)
  94
  95    ! 4.  potential virtual temperature at model surface and observation location
  96    ! ===========================================================================
  97
  98    ! 4.1 potential virtual temperature on model surface
  99    !     ----------------------------------------------
 100
 101    thvm = tvm * (1000.0 / (pm/100.0)) ** rcp
 102
 103    ! 4.2 potential virtual temperature at observation location
 104    !     -----------------------------------------------------
 105
 106    thvo = tvo * (1000.0 / (po/100.0)) ** rcp
 107
 108
 109    ! 5.  bulk richardson number and moni-obukov length
 110    ! =================================================
 111
 112    ! 5.1 Pre-calculations
 113    !     ----------------
 114
 115    za  = hm - ho
 116
 117    ! Because the following formula is derived based on
 118    ! the surface layer height is hs_max=40m. Under
 119    ! free convection, we assume that atmospheric state
 120    ! above the surface layer is fully mixed, and the
 121    ! wind at the lowest sigma half level is same as the
 122    ! wind at top of the surface layer.
 123
 124    if (za > hs_max) za = hs_max
 125
 126    Va2 =   um*um + vm*vm
 127    Vc2 = 4.0 * MAX ((tho - thm), 0.0)
 128
 129    V2  = Va2 + Vc2
 130
 131    ! 5.2 Bulk richardson number
 132    !     ----------------------
 133
 134    rib = (gravity * za / thm) * (thvm - thvo) / V2
 135
 136    ! 5.3 Monin-obukov length
 137    !     -------------------
 138
 139       hll = rib * LOG (za/zint)
 140
 141    ! FIX? is this right?
 142    ! rll = za / hll
 143
 144    ! 5.4 Ratio PBL height/Monin Obukov length: za/rll
 145    !     ------------------------------------
 146
 147    hll =  ABS (hll)
 148
 149
 150    ! 6.  CORRECTION FACTOR BASED UPON REGIME
 151    ! =======================================
 152
 153    ! 6.1 Stable conditions (REGIME 1)
 154    !     ---------------------------
 155
 156    ! correc = 1.0      !  rib > 0.2
 157
 158    ! 6.2 Mechanically driven turbulence (REGIME 2)
 159    !     ------------------------------------------
 160
 161    ! correc = 1.0      !  0.0 =< rib <= 0.2
 162
 163    correc = 1.0
 164
 165    if (rib < 0.0) then
 166
 167       ! 6.3 Unstable Forced convection (REGIME 3)
 168       !     -------------------------------------
 169
 170       ! correc = 1.0  !   hll <= 1.5
 171
 172
 173       ! 6.4 Free convection (REGIME 4)
 174       !     --------------------------
 175
 176       if (hll > 1.5) then
 177          if (zint < 0.2) then
 178             correc = 1.000 + 0.320 * zint ** 0.200
 179          else if (zint >= 0.2) then
 180             correc = 1.169 + 0.315 * zint
 181          end if
 182
 183          ! 6.4.1 The factor depended on Za (MM5, param.F)
 184
 185          winfac = 0.5*(log(za/0.05)/log(40.0/0.05) &
 186                        + log(za)/log(40.0))
 187
 188          correc =  correc * winfac
 189       end if
 190    end if
 191
 192    ! if (trace_use) call da_trace_exit("da_mo_correction")
 193
 194 end function da_mo_correction
 195
 196