var/da/da_setup_structures/da_setup_pseudo_obs.inc

   1 subroutine da_setup_pseudo_obs(grid, iv, ob)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose: Sets up pseudo ob part of observation structure.
   5    !-------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(domain),     intent(in)    :: grid ! Model data
  10    type(iv_type),    intent(inout) :: iv   ! Obs and header structure.
  11    type(y_type),     intent(inout) :: ob   ! (Smaller) observation structure.
  12
  13    integer :: i, j
  14    logical :: outside
  15    real    :: dx, dy, dxm, dym
  16
  17    if ( num_pseudo <= 0 ) return
  18
  19    if (trace_use_dull) call da_trace_entry("da_setup_pseudo_obs")
  20
  21    write(unit=message(1),fmt='(a,a)')'==> pseudo OBS test:'
  22    write(unit=message(2),fmt='(a,a10)')    'pseudo_var = ', trim(pseudo_var)
  23    write(unit=message(3),fmt='(a,f10.5)') 'pseudo_val = ', pseudo_val
  24    write(unit=message(4),fmt='(a,f10.5)') 'pseudo_err = ', pseudo_err
  25    write(unit=message(5),fmt='(a,f10.5)') 'pseudo_x   = ', pseudo_x
  26    write(unit=message(6),fmt='(a,f10.5)') 'pseudo_y   = ', pseudo_y
  27    write(unit=message(7),fmt='(a,f10.5)') 'pseudo_z   = ', pseudo_z
  28    call da_message(message(1:7))
  29
  30    if ( pseudo_tpw .or. pseudo_ztd .or. pseudo_ref ) then
  31       ob%nlocal(pseudo)               = 0
  32       iv%info(pseudo)%nlocal          = 0
  33       iv%info(pseudo)%ntotal          = 0
  34       iv%info(pseudo)%plocal(:)       = 0
  35       iv%info(pseudo)%ptotal(:)       = 0
  36    end if
  37
  38    if ( var4d ) then
  39       ! make sure pseudo_time is within valid range
  40       pseudo_time = min(num_fgat_time,max(1,pseudo_time))
  41       iv%time = pseudo_time
  42    else
  43       iv%time = 1
  44    end if
  45
  46    ! Find out if pseudo ob is local
  47    outside = .false.
  48    i = int(pseudo_x)
  49    j = int(pseudo_y)
  50    if (fg_format == fg_format_kma_global) then
  51       if ((j < jts-1) .or. (j > jte)) outside = .true.
  52    else
  53       if ((i < ids)   .or. (i >= ide) .or. (j < jds)   .or. (j >= jde)) outside = .true.
  54       if ((i < its-1) .or. (i >  ite) .or. (j < jts-1) .or. (j >  jte)) outside = .true.
  55    end if
  56
  57    ! setup gpsref structure
  58    if ( pseudo_ref ) then
  59       iv%info(gpsref)%ptotal(0:iv%time-1)           = 0
  60       iv%info(gpsref)%ptotal(iv%time:max_fgat_time) = num_pseudo
  61       iv%info(gpsref)%ntotal    = num_pseudo
  62       if (outside) then
  63          ob%nlocal(gpsref)               = 0
  64          iv%info(gpsref)%nlocal          = 0
  65          iv%info(gpsref)%plocal(:)       = 0
  66       else
  67          ob%nlocal(gpsref)         = num_pseudo
  68          iv%info(gpsref)%nlocal    = num_pseudo
  69          iv%info(gpsref)%plocal(0:iv%time-1)           = 0
  70          iv%info(gpsref)%plocal(iv%time:max_fgat_time) = num_pseudo
  71          iv%info(gpsref)%max_lev   = 1
  72          allocate(iv%gpsref(1:iv%info(gpsref)%nlocal))
  73          call da_allocate_obs_info(iv, gpsref)
  74          iv%info(gpsref)%proc_domain(:,:) = .true.
  75          iv%info(gpsref)%obs_global_index(iv%info(gpsref)%nlocal) = 1
  76          allocate(iv%gpsref(num_pseudo)%ref(1:1))
  77          allocate(iv%gpsref(num_pseudo)%  h(1:1))
  78          allocate(iv%gpsref(num_pseudo)%  p(1:1))
  79          allocate(iv%gpsref(num_pseudo)%  t(1:1))
  80          allocate(iv%gpsref(num_pseudo)%  q(1:1))
  81          allocate(ob%gpsref(1:num_pseudo))
  82          allocate(ob%gpsref(num_pseudo)%ref(1:num_pseudo))
  83          iv%info(gpsref)%levels(1) = 1
  84          iv%info(gpsref)%x(:,1)   = pseudo_x
  85          iv%info(gpsref)%y(:,1)   = pseudo_y
  86          iv%info(gpsref)%zk(:,1)  = pseudo_z
  87          iv%info(gpsref)%i(:,1)   = int(pseudo_x)
  88          iv%info(gpsref)%j(:,1)   = int(pseudo_y)
  89          iv%info(gpsref)%k(:,:)   = int(pseudo_z)
  90          iv % gpsref(1) %  h(1)   = pseudo_z
  91          iv%info(gpsref)%dx(:,1)  = pseudo_x-real(iv%info(gpsref)%i(1,1))
  92          iv%info(gpsref)%dy(:,1)  = pseudo_y-real(iv%info(gpsref)%j(1,1))
  93          iv%info(gpsref)%dxm(:,1) = 1.0-iv%info(gpsref)%dx(1,1)
  94          iv%info(gpsref)%dym(:,1) = 1.0-iv%info(gpsref)%dy(1,1)
  95          i = iv%info(gpsref)%i(1,1)
  96          j = iv%info(gpsref)%j(1,1)
  97          iv%info(gpsref)%lat(:,1) = grid%xb%lat(i,j)
  98          iv%info(gpsref)%lon(:,1) = grid%xb%lon(i,j)
  99          iv % gpsref(1) %ref(1) % inv   = pseudo_val
 100          iv % gpsref(1) %ref(1) % error = pseudo_err
 101          iv % gpsref(1) %ref(1) % qc    = 0
 102       end if
 103    end if
 104
 105    ! setup gpspw structure
 106    if ( pseudo_tpw .or. pseudo_ztd ) then
 107       iv%info(gpspw)%ptotal(0:iv%time-1)           = 0
 108       iv%info(gpspw)%ptotal(iv%time:max_fgat_time) = num_pseudo
 109       iv%info(gpspw)%ntotal    = num_pseudo
 110       if (outside) then
 111          ob%nlocal(gpspw)               = 0
 112          iv%info(gpspw)%nlocal          = 0
 113          iv%info(gpspw)%plocal(:)       = 0
 114       else
 115          ob%nlocal(gpspw)         = num_pseudo
 116          iv%info(gpspw)%nlocal    = num_pseudo
 117          iv%info(gpspw)%plocal(0:iv%time-1)           = 0
 118          iv%info(gpspw)%plocal(iv%time:max_fgat_time) = num_pseudo
 119          allocate(iv%gpspw(1:iv%info(gpspw)%nlocal))
 120          call da_allocate_obs_info(iv, gpspw)
 121          iv%info(gpspw)%proc_domain(:,:) = .true.
 122          iv%info(gpspw)%obs_global_index(iv%info(gpspw)%nlocal) = 1
 123          allocate (ob % gpspw (1:num_pseudo))
 124          ob % gpspw(1) % tpw   = 0.0
 125          iv%info(gpspw)%x(:,1)   = pseudo_x
 126          iv%info(gpspw)%y(:,1)   = pseudo_y
 127          iv%info(gpspw)%zk(:,1)  = pseudo_z
 128          iv%info(gpspw)%i(:,1)   = int(pseudo_x)
 129          iv%info(gpspw)%j(:,1)   = int(pseudo_y)
 130          iv%info(gpspw)%k(:,1)   = int(pseudo_z)
 131          iv%info(gpspw)%dx(:,1)  = pseudo_x-real(iv%info(gpspw)%i(1,1))
 132          iv%info(gpspw)%dy(:,1)  = pseudo_y-real(iv%info(gpspw)%j(1,1))
 133          iv%info(gpspw)%dxm(:,1) = 1.0 - iv%info(gpspw)%dx(1,1)
 134          iv%info(gpspw)%dym(:,1) = 1.0 - iv%info(gpspw)%dy(1,1)
 135          iv % gpspw(1) % tpw % inv   = pseudo_val
 136          iv % gpspw(1) % tpw % error = pseudo_err
 137          iv % gpspw(1) % tpw % qc    = 0
 138          i = iv%info(gpspw)%i(1,1)
 139          j = iv%info(gpspw)%j(1,1)
 140          iv%info(gpspw)%lat(:,1) = grid%xb%lat(i,j)
 141          iv%info(gpspw)%lon(:,1) = grid%xb%lon(i,j)
 142          if ( pseudo_elv > -999.0 ) then
 143             ! pseudo_elv is set in the namelist
 144             iv%info(gpspw)%elv      = pseudo_elv
 145          else
 146             ! assign model terrain to ob elv
 147             dx  = iv%info(gpspw)%dx(1,1)
 148             dy  = iv%info(gpspw)%dy(1,1)
 149             dxm = iv%info(gpspw)%dxm(1,1)
 150             dym = iv%info(gpspw)%dym(1,1)
 151             iv%info(gpspw)%elv(1) =  &
 152                dym*(dxm*grid%xb%terr(i,j)   + dx*grid%xb%terr(i+1,j)) + &
 153                dy *(dxm*grid%xb%terr(i,j+1) + dx*grid%xb%terr(i+1,j+1))
 154          end if
 155       end if
 156    end if
 157
 158    if ( pseudo_uvtpq ) then
 159       iv%info(pseudo)%ntotal          = num_pseudo
 160       iv%info(pseudo)%ptotal(0:iv%time-1)           = 0
 161       iv%info(pseudo)%ptotal(iv%time:max_fgat_time) = num_pseudo
 162       if (outside) then
 163          ob%nlocal(pseudo)               = 0
 164          iv%info(pseudo)%nlocal          = 0
 165          iv%info(pseudo)%plocal(:)       = 0
 166       else
 167          ob%nlocal(pseudo)               = num_pseudo
 168          iv%info(pseudo)%nlocal          = num_pseudo
 169          iv%info(pseudo)%plocal(0:iv%time-1)           = 0
 170          iv%info(pseudo)%plocal(iv%time:max_fgat_time) = num_pseudo
 171          iv%info(pseudo)%max_lev         = 1
 172          allocate (iv%pseudo(1:iv%info(pseudo)%nlocal))
 173          call da_allocate_obs_info(iv, pseudo)
 174          iv%info(pseudo)%proc_domain(:,:) = .true.
 175          allocate (ob % pseudo (1:num_pseudo))
 176          iv%pseudo(:) % u % inv = missing_r
 177          iv%pseudo(:) % v % inv = missing_r
 178          iv%pseudo(:) % t % inv = missing_r
 179          iv%pseudo(:) % p % inv = missing_r
 180          iv%pseudo(:) % q % inv = missing_r
 181
 182          iv%pseudo(:) % u % error = missing_r
 183          iv%pseudo(:) % v % error = missing_r
 184          iv%pseudo(:) % t % error = missing_r
 185          iv%pseudo(:) % p % error = missing_r
 186          iv%pseudo(:) % q % error = missing_r
 187
 188          iv%pseudo(:) % u % qc  = missing_data
 189          iv%pseudo(:) % v % qc  = missing_data
 190          iv%pseudo(:) % t % qc  = missing_data
 191          iv%pseudo(:) % p % qc  = missing_data
 192          iv%pseudo(:) % q % qc  = missing_data
 193
 194          ob%pseudo(:) % u = missing_r
 195          ob%pseudo(:) % v = missing_r
 196          ob%pseudo(:) % t = missing_r
 197          ob%pseudo(:) % p = missing_r
 198          ob%pseudo(:) % q = missing_r
 199
 200          iv%info(pseudo)%x(:,:)  = pseudo_x
 201          iv%info(pseudo)%y(:,:)  = pseudo_y
 202          iv%info(pseudo)%zk(:,:) = pseudo_z
 203
 204          iv%info(pseudo)%i(:,:) = int(pseudo_x)
 205          iv%info(pseudo)%j(:,:) = int(pseudo_y)
 206          iv%info(pseudo)%k(:,:) = int(pseudo_z)
 207
 208          iv%info(pseudo)%dx(:,:) = pseudo_x-real(iv%info(pseudo)%i(:,:))
 209          iv%info(pseudo)%dy(:,:) = pseudo_y-real(iv%info(pseudo)%j(:,:))
 210          iv%info(pseudo)%dxm(:,:)=1.0-iv%info(pseudo)%dx(:,:)
 211          iv%info(pseudo)%dym(:,:)=1.0-iv%info(pseudo)%dy(:,:)
 212          iv%info(pseudo)%levels(:) = 1
 213
 214          if ( trim(adjustl(pseudo_var)) == 'u' ) then
 215             iv%pseudo(:) % u % inv = pseudo_val
 216             iv%pseudo(:) % u % error = pseudo_err
 217             iv%pseudo(:) % u % qc = 0
 218          else if ( trim(adjustl(pseudo_var)) == 'v' ) then
 219             iv%pseudo(:) % v % inv = pseudo_val
 220             iv%pseudo(:) % v % error = pseudo_err
 221             iv%pseudo(:) % v % qc = 0
 222          else if ( trim(adjustl(pseudo_var)) == 't' ) then
 223             iv%pseudo(:) % t % inv = pseudo_val
 224             iv%pseudo(:) % t % error = pseudo_err
 225             iv%pseudo(:) % t % qc = 0
 226          else if ( trim(adjustl(pseudo_var)) == 'p' ) then
 227             iv%pseudo(:) % p % inv = pseudo_val
 228             iv%pseudo(:) % p % error = pseudo_err
 229             iv%pseudo(:) % p % qc = 0
 230          !else if ( trim(adjustl(pseudo_var)) == 'q' ) then
 231          else if ( pseudo_var(1:1) == 'q' ) then
 232             iv%pseudo(:) % q % inv = pseudo_val
 233             iv%pseudo(:) % q % error = pseudo_err
 234             iv%pseudo(:) % q % qc = 0
 235          end if
 236       end if
 237    end if
 238
 239    call da_allocate_y (iv, ob)
 240
 241    if (trace_use_dull) call da_trace_exit("da_setup_pseudo_obs")
 242
 243 end subroutine da_setup_pseudo_obs
 244
 245