var/da/da_tools/da_sfc_hori_interp_weights.inc

   1 subroutine da_sfc_hori_interp_weights(n, info, sfc_ob, xb, ob_land_type)
   2
   3    implicit none
   4
   5    integer        ,  intent(in)    :: n
   6    type(infa_type),  intent(inout) :: info
   7    type(synop_type), intent(inout) :: sfc_ob
   8    type(xb_type),    intent(in)    :: xb
   9    integer        ,  intent(in)    :: ob_land_type !1:land, 0:water
  10
  11    integer :: i, j, iloc, ii, jj
  12    integer :: iloc_min_hdiff
  13    integer :: landmask(4)
  14    real    :: dx, dxm, dy, dym
  15    real    :: hh(4), w(4), hdiff(4)
  16
  17 ! 1: (i,   j)
  18 ! 2: (i+1, j)
  19 ! 3: (i,   j+1)
  20 ! 4: (i+1, j+1)
  21 !------------------------------------------------------
  22 !  (i,j+1)
  23 !    -----------------(i+1,j+1)
  24 !    |   w2          |
  25 !    |dym        w1  |
  26 !    |      obs      |
  27 !    |------- *      |
  28 !    |dy  w4  |  w3  |
  29 !    |   dx   | dxm  |
  30 !    |----------------
  31 !   (i,j)            (i+1,j)
  32 !
  33 !--------------------------------------------------------
  34
  35    ! initialize
  36    ii = nint(info%x(1,n))
  37    jj = nint(info%y(1,n))
  38
  39    i   = info%i(1,n)
  40    j   = info%j(1,n)
  41    dx  = info%dx(1,n)
  42    dy  = info%dy(1,n)
  43    dxm = info%dxm(1,n)
  44    dym = info%dym(1,n)
  45
  46    w(1) = dym*dxm   ! weight for point (i,  j)
  47    w(2) = dym*dx    ! weight for point (i+1,j)
  48    w(3) = dy *dxm   ! weight for point (i,  j+1)
  49    w(4) = dy *dx    ! weight for point (i+1,j+1)
  50
  51    ! lowest half model level height
  52    hh(1) = xb%h(i,   j,   kts)
  53    hh(2) = xb%h(i+1, j,   kts)
  54    hh(3) = xb%h(i,   j+1, kts)
  55    hh(4) = xb%h(i+1, j+1, kts)
  56
  57    ! landmask: 1-land, 0-water
  58    landmask(1) = int(xb%landmask(i,   j))
  59    landmask(2) = int(xb%landmask(i+1, j))
  60    landmask(3) = int(xb%landmask(i,   j+1))
  61    landmask(4) = int(xb%landmask(i+1, j+1))
  62
  63    do iloc = 1, 4
  64       hdiff(iloc) = abs(hh(iloc) - sfc_ob%h)
  65       if ( landmask(iloc) /= ob_land_type ) then
  66          hdiff(iloc) = 8888.0 !set a large value to not use it
  67       end if
  68    end do
  69    if ( sum(hdiff)/4.0 == 8888.0 ) then
  70       ! do not use the ob, because the model land type is not
  71       ! consistent with the ob land type
  72       !print*, 'Rejecting ', trim(info%id(n)), ' due to land type mismatch'
  73       sfc_ob%p%qc = missing_data
  74       sfc_ob%t%qc = missing_data
  75       sfc_ob%q%qc = missing_data
  76       sfc_ob%u%qc = missing_data
  77       sfc_ob%v%qc = missing_data
  78       return
  79    else
  80       iloc_min_hdiff = 1 !initialize
  81       do iloc = 2, 4
  82          if ( hdiff(iloc) > hdiff(iloc_min_hdiff) ) then
  83             cycle
  84          else
  85             if ( hdiff(iloc) < hdiff(iloc_min_hdiff) ) then
  86                iloc_min_hdiff = iloc
  87             else
  88                if ( w(4-iloc+1) < w(4-iloc_min_hdiff+1) ) then
  89                   iloc_min_hdiff = iloc
  90                end if
  91             end if
  92          end if
  93       end do
  94       select case ( iloc_min_hdiff )
  95          case ( 1 )
  96             ii = i
  97             jj = j
  98          case ( 2 )
  99             ii = i+1
 100             jj = j
 101          case ( 3 )
 102             ii = i
 103             jj = j+1
 104          case ( 4 )
 105             ii = i+1
 106             jj = j+1
 107       end select
 108    end if
 109
 110    if ( ii == i .and. jj == j ) then
 111       dx  = 0.0
 112       dxm = 1.0
 113       dy  = 0.0
 114       dym = 1.0
 115    else if ( ii == i+1 .and. jj == j ) then
 116       dx  = 1.0
 117       dxm = 0.0
 118       dy  = 0.0
 119       dym = 1.0
 120    else if ( ii == i   .and. jj == j+1 ) then
 121       dx  = 0.0
 122       dxm = 1.0
 123       dy  = 1.0
 124       dym = 0.0
 125    else if ( ii == i+1 .and. jj == j+1 ) then
 126       dx  = 1.0
 127       dxm = 0.0
 128       dy  = 1.0
 129       dym = 0.0
 130    end if
 131
 132    !print*, 'Resetting i, j for : ', trim(info%id(n)), ii, jj
 133    ! re-assign the weights, lowest level only
 134    info%dx(1,n)  = dx
 135    info%dy(1,n)  = dy
 136    info%dxm(1,n) = dxm
 137    info%dym(1,n) = dym
 138
 139 !         call da_detsurtyp ( grid%xb%snow, grid%xb%xice, grid%xb%landmask,  &
 140 !            grid%xb%ivgtyp, grid%xb%isltyp, &
 141 !            ims, ime, jms, jme, &
 142 !            i, j, dx, dy, dxm, dym, &
 143 !            model_isflg,model_ivgtyp, model_isltyp, &
 144 !            Surface(1)%Water_Coverage, Surface(1)%Ice_Coverage, &
 145 !            Surface(1)%Land_Coverage, Surface(1)%Snow_Coverage )
 146
 147 end subroutine da_sfc_hori_interp_weights