var/da/da_tools/da_obs_sfc_correction.inc

   1 subroutine da_obs_sfc_correction(info, sfc_obs, n, xb)
   2
   3    !--------------------------------------------------------------------
   4    ! Purpose: correct the surface measurements (wind,
   5    ! temperature, and pressure) from the observed height to the WRF
   6    ! model's lowest half-zeta level before going to the minimization.
   7    !
   8    !   Wind       : based on the similarity theory
   9    !   Temperature: Frank Ruggiero's (1996) method
  10    !   Pressure   : Hydrostatic correction
  11    !
  12    ! The order of the vertical index is "kts=1(bottom) and kte(top)".
  13    ! With cv_options=2 and sfc_assi_option=1, this procedure must be
  14    ! gone through, otherwise unrealistic results may be obtained.
  15    !--------------------------------------------------------------------
  16
  17    implicit none
  18
  19    type(infa_type),  intent(in)    :: info
  20    type(synop_type), intent(inout) :: sfc_obs
  21    integer,          intent(in)    :: n
  22    type(xb_type),    intent(in)    :: xb
  23
  24    real    :: roughness, psfc, mslp, dx, dxm, dy, dym, ho, po, to, qo
  25    real    :: hm, pm, tm, qm, um, vm, correc, val
  26
  27    integer :: i, j, k
  28    real    :: t_mdl(kts:kte)
  29    real    :: q_mdl(kts:kte)
  30    real    :: u_mdl(kts:kte)
  31    real    :: v_mdl(kts:kte)
  32    real    :: height(kts:kte)
  33    real    :: pressure(kts:kte)
  34
  35    if (trace_use_dull) call da_trace_entry("da_obs_sfc_correction")
  36
  37    ! 1. Check if it needs to do the surface correction at the first level
  38
  39    ! 1.1 Surface reports located at far below the lowest model level
  40
  41    ! 2. Model profile at OBS site for surface correction
  42
  43    i   = info%i(1,n)
  44    j   = info%j(1,n)
  45    dx  = info%dx(1,n)
  46    dy  = info%dy(1,n)
  47    dxm = info%dxm(1,n)
  48    dym = info%dym(1,n)
  49
  50    ! Model roughness at the obs site
  51
  52    roughness = dym*(dxm*xb%rough(i,j)   + dx*xb%rough(i+1,j)) &
  53       + dy *(dxm*xb%rough(i,j+1) + dx*xb%rough(i+1,j+1))
  54
  55    do k = kts, kte
  56       pressure(k) = dym*(dxm*xb%p(i,j,k) + dx*xb%p(i+1,j,k)) + dy*(dxm*xb%p(i,j+1,k) + dx*xb%p(i+1,j+1,k))
  57       height(k)   = dym*(dxm*xb%h(i,j,k) + dx*xb%h(i+1,j,k)) + dy*(dxm*xb%h(i,j+1,k) + dx*xb%h(i+1,j+1,k))
  58       t_mdl(k)    = dym*(dxm*xb%t(i,j,k) + dx*xb%t(i+1,j,k)) + dy*(dxm*xb%t(i,j+1,k) + dx*xb%t(i+1,j+1,k))
  59       q_mdl(k)    = dym*(dxm*xb%q(i,j,k) + dx*xb%q(i+1,j,k)) + dy*(dxm*xb%q(i,j+1,k) + dx*xb%q(i+1,j+1,k))
  60       u_mdl(k)    = dym*(dxm*xb%u(i,j,k) + dx*xb%u(i+1,j,k)) + dy*(dxm*xb%u(i,j+1,k) + dx*xb%u(i+1,j+1,k))
  61       v_mdl(k)    = dym*(dxm*xb%v(i,j,k) + dx*xb%v(i+1,j,k)) + dy*(dxm*xb%v(i,j+1,k) + dx*xb%v(i+1,j+1,k))
  62    end do
  63
  64    ! 3. OBS data and recover the surface pressure from the
  65    ! mean sea level pressure (mslp)
  66
  67    ho   = sfc_obs % h
  68    po   = sfc_obs % p % inv
  69    to   = sfc_obs % t % inv
  70    qo   = sfc_obs % q % inv
  71
  72    ! 3.1 Compute the surface OBS pressure from mean sea level pressure
  73
  74    if ( psfc_from_slp ) then
  75       mslp = info%slp(n)%inv
  76       if (abs(mslp - missing_r) > 1.0) then
  77          psfc = missing_r
  78          if (abs(ho - missing_r) > 1.0) then
  79             if (abs(to - missing_r) > 1.0) then
  80                call da_sfcprs (kts, kte, pressure, t_mdl, q_mdl, height, psfc, mslp, ho, to)
  81             else
  82                call da_sfcprs (kts, kte, pressure, t_mdl, q_mdl, height, psfc, mslp, ho)
  83             end if
  84          end if
  85          sfc_obs % p % inv = psfc
  86          ! YRG: to allow assmilate the Psfc from mslp:
  87          sfc_obs % p % qc  = 0
  88       end if
  89       po = sfc_obs % p % inv
  90    end if
  91
  92    if (sfc_obs % p % inv < 1.0) then
  93       sfc_obs % p % qc  = missing_data
  94    end if
  95
  96    po = sfc_obs % p % inv
  97
  98    ! 3.2 Check that obs pressure and height are present
  99    !     ----------------------------------------------
 100
 101    if (abs(po - missing_r) < 1.0  .OR. abs(ho - missing_r) < 1.0) then
 102       if (trace_use_dull) call da_trace_exit("da_obs_sfc_correction")
 103       return
 104
 105       ! write(unit=message(1), fmt='(/3(1x,a))') &
 106       !    'MISSinG HEIGHT OR PRESSURE OBSERVATION ID ', &
 107       !    trim (sfc_obs%info % id), trim (sfc_obs%info % name)
 108
 109       ! write(unit=message(2), fmt='(2(A,F12.3,/))') &
 110       !                         ' height   = ',ho,&
 111       !                         ' pressure = ',po
 112       ! call da_error(__FILE__,__LINE__,message(1:2))
 113
 114    end if
 115
 116    ! 4.  Bring surface observation below model levels with good height quality
 117    ! ================================================
 118
 119    if (sfc_obs % h < height(kts)) then
 120
 121       ! 2.3 Make use of local variables for model
 122       !     -------------------------------------
 123
 124       um = u_mdl(kts)
 125       vm = v_mdl(kts)
 126       tm = t_mdl(kts)
 127       pm = pressure (kts)
 128       qm = q_mdl(kts)
 129       hm = height(kts)
 130
 131       ! 3.2 Correction wind based on similarity laws
 132       !     -------------------------------------------
 133
 134       if ((abs(sfc_obs%u%inv - missing_r) > 1.0) .AND. (abs(sfc_obs%v%inv - missing_r) > 1.0)) then
 135          ! 3.2.1 Correction factor
 136
 137          ! temperature and moisture at sigma level:
 138
 139          if (abs(to - missing_r) < 1.0) then
 140             correc = da_mo_correction(ho, po, tm, qo, hm, pm, tm, qm, um ,vm, roughness)
 141          else
 142             correc = da_mo_correction(ho, po, to, qo, hm, pm, tm, qm, um ,vm, roughness)
 143          end if
 144
 145          ! 3.2.2 Wind correction
 146          !       ---------------
 147
 148          !  Correct data and replace any previous wind qcs
 149          !  with surface correction flag
 150
 151          sfc_obs % u % inv = correc * sfc_obs % u % inv
 152          if ( abs(correc-1.0)>0.0 ) sfc_obs % u % qc  = surface_correction
 153
 154          sfc_obs % v % inv = correc * sfc_obs % v % inv
 155          if ( abs(correc-1.0)>0.0 ) sfc_obs % v % qc  = surface_correction
 156       end if
 157
 158       ! 3.4 Correct pressure
 159       !     ----------------
 160
 161       if (sfc_obs % p % qc >= 0) then
 162          !  Correct data
 163          if (abs(to  - missing_r) > 1.0 .and. abs(qo - missing_r) > 1.0) then
 164             call da_intpsfc_prs (val, ho, po, hm, tm, qm, to, qo)
 165          else if (abs(to  - missing_r) > 1.0) then
 166             call da_intpsfc_prs (val, ho, po, hm, tm, qm, to)
 167          else
 168             call da_intpsfc_prs (val, ho, po, hm, tm, qm)
 169          end if
 170
 171          !  Replace any previous pressure qc by the surface correction
 172
 173          sfc_obs % p % inv = val
 174          sfc_obs % p % qc  = surface_correction
 175       end if
 176
 177       ! 3.5 Correct temperature
 178       !     -------------------
 179
 180       if (abs(sfc_obs % t % inv - missing_r) > 1.0) then
 181          !  Correct data
 182          call da_intpsfc_tem (val, ho, po, to, height, pressure, t_mdl, kts, kte)
 183
 184          sfc_obs % t % inv = val
 185
 186          !  Replace any previous temperature qc by the surface correction
 187          sfc_obs % t % qc  = surface_correction
 188       end if
 189
 190       ! 3.6 Assign model lowest level height + 1m to observation
 191       !      -----------------------------------------------------
 192       ! sfc_obs % h = height(kts) + 1.0
 193
 194       ! 3.7 Height QC
 195       !     ---------
 196       ! sfc_obs % height_qc = surface_correction
 197    end if
 198
 199    if (trace_use_dull) call da_trace_exit("da_obs_sfc_correction")
 200
 201 end  subroutine da_obs_sfc_correction
 202
 203