var/da/da_gpseph/da_gpseph_create_ob.inc

   1 SUBROUTINE da_gpseph_create_ob ( nlevel, h, ref, lat, lon, azim, pseudo_ob, lowest_level )
   2
   3 !---------------------------------------------------------------------------------------
   4 !  Purpose: calculate pseudo obs (refractivity) at the mean height of each model level.
   5 !           (Chen et al., 2009, doi:10.3319/TAO.2007.11.29.01(F3C))
   6 !           It is called by da_read_obs_bufrgpsro_eph (when ob_format_gpsro=1)
   7 !           or da_read_obs_ascii (when ob_format_gpsro=2).
   8 !---------------------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    type (ob_in_mean_h), intent(inout) :: pseudo_ob
  13    integer,             intent(out)   :: lowest_level
  14    integer,             intent(in)    :: nlevel
  15    real,                intent(in)    :: h(nlevel)
  16    real,                intent(in)    :: ref(nlevel)
  17    real,                intent(in)    :: lat(nlevel)
  18    real,                intent(in)    :: lon(nlevel)
  19    real,                intent(in)    :: azim(nlevel)
  20
  21    integer :: i, j, k
  22    integer :: nlevel_valid, icount
  23    real    :: nstart, nstop
  24
  25    real, dimension(kds:kde) :: mean_h
  26
  27    real, allocatable :: log_ref(:)   !nlevel_valid
  28    real, allocatable :: h_km(:)      !nlevel_valid
  29    real, allocatable :: lat_used(:)  !nlevel_valid
  30    real, allocatable :: lon_used(:)  !nlevel_valid
  31    real, allocatable :: azim_used(:) !nlevel_valid
  32
  33    real, dimension(interpolate_level) :: tmp_z   ! obs. levels are interpolated
  34                                                  ! into levels of interpolate_level
  35    real, dimension(interpolate_level) :: tmp_ref
  36    real, dimension(interpolate_level) :: tmp_lat
  37    real, dimension(interpolate_level) :: tmp_lon
  38    real, dimension(interpolate_level) :: tmp_azim
  39
  40    lowest_level = -1
  41
  42    pseudo_ob%ref(:) = 0.0
  43    pseudo_ob%lat(:) = 0.0
  44    pseudo_ob%lon(:) = 0.0
  45    pseudo_ob%azim(:)= 0.0
  46
  47    ! exclude the levels with invalid ref and h data
  48    ! first find out the nlevel_valid for allocating the working arrays
  49    icount = 0
  50    do k = 1, nlevel
  51       if ( (ref(k) > 0.0) .and. (h(k) > 0.0) .and. &
  52            (abs(lat(k)) <= 90.0)  .and.            &
  53            (abs(lon(k)) <= 180.0) .and.            &
  54            (abs(azim(k)) <= 180.0) ) then
  55          icount = icount + 1
  56       end if
  57    end do
  58    nlevel_valid = icount
  59    if ( nlevel_valid < 1 ) return
  60    allocate (log_ref  (nlevel_valid))
  61    allocate (h_km     (nlevel_valid))
  62    allocate (lat_used (nlevel_valid))
  63    allocate (lon_used (nlevel_valid))
  64    allocate (azim_used(nlevel_valid))
  65
  66    ! copy the valid levels to the working arrays
  67    icount = 0
  68    do k = 1, nlevel
  69       if ( (ref(k) > 0.0) .and. (h(k) > 0.0) .and. &
  70            (abs(lat(k)) <= 90.0)  .and.            &
  71            (abs(lon(k)) <= 180.0) .and.            &
  72            (abs(azim(k)) <= 180.0) ) then
  73          icount = icount + 1
  74          log_ref(icount)   = log(ref(k))
  75          h_km(icount)      = h(k)*0.001  !m to km, observation height
  76          lat_used(icount)  = lat(k)
  77          lon_used(icount)  = lon(k)
  78          azim_used(icount) = azim(k)
  79       end if
  80    end do
  81
  82    mean_h(:) = global_h_mean(:) !km, model mean height
  83
  84    do k=kts,kte
  85       if ( mean_h(k) > h_km(1) ) exit
  86    end do
  87    lowest_level = k
  88
  89    ! interpolate_level =2000, tmp_z(k) increase with constant interval(0.01 km) in the vertical direction
  90    !hcl-todo make 0.01 and 20km_top (used to derive interpolate_level) namelist options
  91    do k=1,interpolate_level
  92       tmp_z(k) = h_km(1)+0.01*(k-1)
  93    end do
  94
  95    ! interpolate variables and smooth them on the standard height grid
  96    CALL lintp( nlevel_valid, h_km, log_ref,          &
  97                interpolate_level, tmp_z, tmp_ref )
  98    CALL lintp( nlevel_valid, h_km, lat_used,         &
  99                interpolate_level, tmp_z, tmp_lat )
 100    CALL lintp( nlevel_valid, h_km, lon_used,         &
 101                interpolate_level, tmp_z, tmp_lon )
 102    CALL lintp( nlevel_valid, h_km, azim_used,        &
 103                interpolate_level, tmp_z, tmp_azim )
 104
 105    do i = lowest_level+1, kte-1
 106
 107       do j = 1, interpolate_level
 108          if (tmp_z(j) < (mean_h(i)+mean_h(i-1))/2.0) nstart = j
 109          if (tmp_z(j) < (mean_h(i)+mean_h(i+1))/2.0) nstop  = j
 110       end do
 111
 112       pseudo_ob%ref(i) = 0.0
 113       pseudo_ob%lat(i) = 0.0
 114       pseudo_ob%lon(i) = 0.0
 115       pseudo_ob%azim(i)= 0.0
 116
 117       do j = nstart, nstop
 118          pseudo_ob%ref(i)  = pseudo_ob%ref(i)  + tmp_ref(j)
 119          pseudo_ob%lat(i)  = pseudo_ob%lat(i)  + tmp_lat(j)
 120          pseudo_ob%lon(i)  = pseudo_ob%lon(i)  + tmp_lon(j)
 121          pseudo_ob%azim(i) = pseudo_ob%azim(i) + tmp_azim(j)
 122       end do
 123
 124       pseudo_ob%ref(i)  = pseudo_ob%ref(i) / (nstop-nstart+1)
 125       pseudo_ob%lat(i)  = pseudo_ob%lat(i) / (nstop-nstart+1)
 126       pseudo_ob%lon(i)  = pseudo_ob%lon(i) / (nstop-nstart+1)
 127       pseudo_ob%azim(i) = pseudo_ob%azim(i)/ (nstop-nstart+1)
 128
 129    end do
 130
 131    do k = lowest_level+1, kte-1
 132       pseudo_ob%ref(k) = exp(pseudo_ob%ref(k))
 133       if ( pseudo_ob%lon(k) < 0.0 ) then
 134          pseudo_ob%lon(k) = pseudo_ob%lon(k) + 360.0
 135       end if
 136    end do
 137
 138    deallocate (log_ref  )
 139    deallocate (h_km     )
 140    deallocate (lat_used )
 141    deallocate (lon_used )
 142    deallocate (azim_used)
 143
 144 END SUBROUTINE da_gpseph_create_ob
 145
 146 ! ................................................................
 147
 148 SUBROUTINE lintp ( n1, x1, y1, n2, x2, y2 )
 149
 150 ! **********************************************************
 151 !  This subroutine provides linear interpolation of function
 152 !  y(x) from one grid x1(i), onto another grid x2(i)
 153 !
 154 !  y1(i), i=1,n1 are input grid function values
 155 !  y2(i), i=1,n2 are output grid function values
 156 ! **********************************************************
 157 !  both x1(i) and x2(i) should increase with 'i'
 158 ! **********************************************************
 159       implicit none
 160       integer :: i, j, k, l, m, n
 161       integer, intent(in)::n1
 162       integer, intent(in)::n2
 163       real,dimension(n1),intent(in)::x1
 164       real,dimension(n2),intent(in)::x2
 165       real,dimension(n1),intent(in)::y1
 166       real,dimension(n2),intent(out)::y2
 167
 168       do i=1,n2
 169       if (x2(i).le.x1(1)) then
 170       y2(i)=y1(1)+(x2(i)-x1(1))*(y1(2)-y1(1))/(x1(2)-x1(1))
 171       else if (x2(i).ge.x1(n1)) then
 172       y2(i)=y1(n1)+(x2(i)-x1(n1))*(y1(n1)-y1(n1-1))/(x1(n1)-x1(n1-1))
 173       else
 174
 175       do j=1,n1
 176       k=j
 177       if ((x1(j).le.x2(i)).and.(x1(j+1).ge.x2(i))) goto 1
 178       end do
 179
 180 1     y2(i)=y1(k)+(x2(i)-x1(k))*(y1(k+1)-y1(k))/(x1(k+1)-x1(k))
 181
 182       end if
 183       end do
 184
 185 END SUBROUTINE lintp
 186