var/da/da_radiance/da_read_kma1dvar.inc

   1 subroutine da_read_kma1dvar (inst,iv, ob, infile)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: read in kma 1dvar innovation output to innovation and obs structure
   5    !
   6    ! METHOD: use F90 sequantial data structure to avoid read file twice
   7    !          so that da_scan_bufrtovs is not necessary any more.
   8    !          1. read radiance data in sequential data structure
   9    !          2. assign sequential data structure to innovation structure
  10    !                and deallocate sequential data structure
  11    !---------------------------------------------------------------------------
  12
  13    implicit none
  14
  15    ! subroutine argument
  16    integer           ,  intent (in)    :: inst
  17    character(20)     ,  intent (in)    :: infile
  18    type (y_type)     ,  intent (inout) :: ob
  19    type (iv_type)    ,  intent (inout) :: iv
  20
  21    ! local variables
  22    integer          :: iost, n, i, j,inunit
  23
  24    ! Declare local variables
  25    logical outside,inside_halo
  26    character(20)  :: instrument
  27
  28    ! pixel information
  29    integer   ::  idate, npass, scanpos
  30    real      ::  rlat, rlon                         !  lat/lon in degrees   for Anfovs
  31    real      ::  satzen      !  scan angles for Anfovs
  32    integer   ::  landsea_mask
  33    real      ::  srf_height
  34    integer   ::  nchanl,knchan                      !  number of channels
  35
  36    real      :: fastem(5)
  37    real , allocatable :: &         !  bright temperatures
  38                          otb(:),oerr(:),inov(:), emiss(:), &
  39                          tb(:), inv(:),err(:)
  40    integer, allocatable :: kchan(:),qc(:)
  41
  42    type (datalink_type), pointer    :: head, p, current
  43
  44    integer                        ::  size, error
  45    type(info_type)                ::  info
  46    type(model_loc_type)           ::  loc
  47
  48    if (trace_use) call da_trace_entry("da_read_kma1dvar")
  49
  50    !**************************************************************************
  51    ! Initialize variables
  52
  53    nchanl = iv%instid(inst)%nchan
  54
  55    allocate ( kchan(nchanl) )
  56    allocate ( otb(nchanl) )
  57    allocate ( oerr(nchanl) )
  58    allocate ( inov(nchanl) )
  59    allocate ( emiss(nchanl) )
  60
  61    allocate ( tb(nchanl) )
  62    allocate ( inv(nchanl) )
  63    allocate ( err(nchanl) )
  64    allocate ( qc(nchanl) )
  65
  66    ! 0.0  Open unit for kma 1dvar file and read file header
  67    !--------------------------------------------------------------
  68
  69    call da_get_unit(inunit)
  70    open(unit=inunit,file=trim(infile),form='formatted', &
  71        iostat = iost, status = 'old')
  72    if (iost /= 0) then
  73       call da_error(__FILE__,__LINE__,&
  74          (/"Cannot open file "//infile/))
  75    end if
  76
  77    read(unit=inunit,fmt='(A,I10,I12)') instrument,idate,npass
  78
  79    ! Loop to read pixel and assign information to a sequential structure
  80    !-------------------------------------------------------------------------
  81
  82    allocate ( head )
  83    nullify  ( head % next )
  84    p => head
  85    size = 0
  86
  87    do j=1,npass
  88       read(inunit,'(2F12.4)')    rlat,rlon
  89       read(inunit,'(I5,20I5)')   knchan,(kchan(i),i=1,knchan)
  90       ! landsea_mask: 0:land ; 1:sea (same as RTTOV)
  91       read(inunit,'(I5,F12.4,I5,F12.4)')  scanpos,satzen,landsea_mask,srf_height
  92       read(inunit,'(20F12.4)') (oerr(i),i=1,knchan)
  93       read(inunit,'(20F12.4)') (emiss(i),i=1,knchan)
  94       read(inunit,'(5F12.4)')  (fastem(i),i=1,5)
  95       read(inunit,'(20F12.4)') (otb(i),i=1,knchan)
  96       read(inunit,'(20F12.4)') (inov(i),i=1,knchan)
  97
  98       ! 1.0 Extract observation location and other required information
  99       !     judge if data is in the domain, read next record if not
 100       !-------------------------------------------------------------------
 101
 102       info%lat  =  rlat
 103       info%lon  =  rlon
 104       call da_llxy (info, loc, outside, inside_halo )
 105       if (outside) cycle
 106
 107       info%elv = srf_height
 108       write(unit=info%date_char, fmt='(i10)') idate
 109
 110       tb(1:nchanl) = missing_r
 111       inv(1:nchanl) = missing_r
 112       err(1:nchanl) = missing_r
 113       qc(1:nchanl) = qc_bad
 114
 115       do i=1,knchan
 116          tb(kchan(i)) = otb(i)
 117          inv(kchan(i)) = inov(i)
 118          err(kchan(i)) = oerr(i)
 119          qc(kchan(i)) = qc_good
 120       end do
 121
 122       !  2.0   assign information to sequential radiance structure
 123       !--------------------------------------------------------------------
 124
 125       allocate (p % tb_inv (1:nchanl))
 126       allocate (p % tb_error (1:nchanl))
 127       allocate (p % tb_qc (1:nchanl))
 128       allocate (p % tb_ob (1:nchanl))
 129       p%info             = info
 130       p%loc              = loc
 131       p%landsea_mask     = landsea_mask
 132       p%scanpos          = scanpos
 133       p%satzen           = satzen
 134       p%satazi           = missing_r
 135       p%solzen           = missing_r
 136       p%tb_inv(1:nchanl)     = inv(1:nchanl)
 137       p%tb_error(1:nchanl)   = err(1:nchanl)
 138       p%tb_qc(1:nchanl)      = qc(1:nchanl)
 139       p%tb_ob(1:nchanl)      = tb(1:nchanl)
 140       p%sensor_index         = inst
 141
 142       size = size + 1
 143       allocate ( p%next, stat=error)   ! add next data
 144       if (error /= 0 ) then
 145          call da_error(__FILE__,__LINE__, &
 146             (/"Cannot allocate radiance structure"/))
 147       end if
 148
 149       p => p%next
 150       nullify (p%next)
 151    end do
 152
 153    iv%total_rad_pixel   = iv%total_rad_pixel + size
 154    iv%total_rad_channel = iv%total_rad_channel + size*nchanl
 155
 156    deallocate (kchan)
 157    deallocate (otb)
 158    deallocate (oerr)
 159    deallocate (inov)
 160    deallocate (emiss)
 161
 162    deallocate (tb)
 163    deallocate (inv)
 164    deallocate (err)
 165    deallocate (qc)
 166
 167    close(inunit)
 168    call da_free_unit(inunit)
 169
 170    !  3.0 allocate innovation and obs radiance structure
 171    !----------------------------------------------------------------
 172    iv%instid(inst)%num_rad = size
 173    ob%instid(inst)%num_rad = size
 174
 175    write(unit=stdout,fmt='(a,i3,2x,a,3x,i10)')  ' allocating space for radiance innov structure', &
 176                              i, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 177
 178    !  4.0 assign sequential structure to innovation structure
 179    !-------------------------------------------------------------
 180    p => head
 181    allocate (iv%instid(inst)%tb_inv(1:nchanl,size))
 182    allocate (iv%instid(inst)%tb_error(1:nchanl,size))
 183    allocate (iv%instid(inst)%tb_qc(1:nchanl,size))
 184    allocate (ob%instid(inst)%tb(1:nchanl,size))
 185    do n = 1, size
 186
 187       iv%instid(i)%info%name(n)  = p%info%name
 188       iv%instid(i)%info%platform(n)  = p%info%platform
 189       iv%instid(i)%info%id(n)    = p%info%id
 190       iv%instid(i)%info%date_char(n) = p%info%date_char
 191       iv%instid(i)%info%levels(n)    = p%info%levels
 192       iv%instid(i)%info%lat(:,n)         = p%info%lat
 193       iv%instid(i)%info%lon(:,n)         = p%info%lon
 194       iv%instid(i)%info%elv(n)   = p%info%elv
 195       iv%instid(i)%info%pstar(n)     = p%info%pstar
 196
 197       iv%instid(inst)%info%i(:,n)    = p%loc%i
 198       iv%instid(inst)%info%j(:,n)    = p%loc%j
 199       iv%instid(inst)%info%dx(:,n)   = p%loc%dx
 200       iv%instid(inst)%info%dy(:,n)   = p%loc%dy
 201       iv%instid(inst)%info%dxm(:,n)  = p%loc%dxm
 202       iv%instid(inst)%info%dym(:,n)  = p%loc%dym
 203
 204       iv%instid(inst)%landsea_mask(n) = p%landsea_mask
 205       iv%instid(inst)%scanpos(n)      = p%scanpos
 206       iv%instid(inst)%satzen(n) = p%satzen
 207       iv%instid(inst)%satazi(n) = p%satazi
 208       iv%instid(inst)%solzen(n) = p%solzen
 209       iv%instid(inst)%tb_inv(1:nchanl,n)   = p%tb_inv(1:nchanl)
 210       iv%instid(inst)%tb_error(1:nchanl,n) = p%tb_error(1:nchanl)
 211       iv%instid(inst)%tb_qc(1:nchanl,n)    = p%tb_qc(1:nchanl)
 212       ob%instid(inst)%tb(1:nchanl,n)       = p%tb_ob(1:nchanl)
 213
 214       ! write(unit=stdout,*) inst, nread, iv%instid(inst)%tb_inv(1:nchanl,n)
 215       current => p
 216       p => p%next
 217
 218       ! free current data
 219       deallocate (current % tb_ob)
 220       deallocate (current % tb_inv)
 221       deallocate (current % tb_error)
 222       deallocate (current % tb_qc)
 223       deallocate (current)
 224    end do
 225
 226    ! check if sequential structure has been freed
 227    !
 228    ! p => head
 229    ! do i = 1, size
 230    !    write (unit=stdout,fmt=*)  i, p%tb(1:nchanl)%inv
 231    !    p => p%next
 232    ! end do
 233
 234    if (trace_use) call da_trace_exit("da_readkma1dvar")
 235
 236 end subroutine da_read_kma1dvar
 237
 238