var/da/da_radiance/da_read_obs_ncgoesimg.inc

   1 subroutine  da_read_obs_ncgoesimg (iv,infile)
   2 !----------------------------------------------------
   3 ! Purpose:  read netcdf format goes imager GVAR data
   4 !           and convert to brightness temperature
   5 !
   6 ! Method: Yang et al., 2017: Impact of assimilating GOES imager clear-sky radiance
   7 !         with a rapid refresh assimilation system for convection-permitting forecast
   8 !         over Mexico. J. Geophys. Res. Atmos., 122, 5472–5490
   9 !----------------------------------------------------
  10    implicit none
  11
  12    character(100),   intent (in) :: infile
  13    type (iv_type),intent (inout) :: iv
  14
  15    real(kind=8)                  ::  obs_time
  16
  17    type(datalink_type), pointer  :: head, p, current, prev
  18    type(info_type)               :: info
  19    type(model_loc_type)          :: loc
  20
  21 ! Number of channels for sensors in NETCDF
  22
  23    integer(i_kind),parameter    :: nchan = 4
  24    integer(i_kind)              :: inst,platform_id,satellite_id,sensor_id
  25    integer(i_kind), allocatable :: ptotal(:)
  26    real(r_kind)                 :: crit
  27    integer(i_kind)              :: ifgat, iout, iobs
  28    logical                      :: outside, outside_all, iuse
  29
  30 ! Variables for netcdf IO
  31    logical                      :: found, head_found, isfile
  32    integer(i_kind)              :: iret,ireadsb,ireadmg,irec,isub,next
  33    character(10)                :: date
  34
  35 ! Work variables for time
  36    integer(i_kind)              :: idate, im, iy, idd, ihh
  37    integer(i_kind)              :: idate5(6)
  38
  39 ! Other work variables
  40
  41    real(r_kind)     :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg
  42    real(r_kind)     :: timedif, pred, crit1, dist1
  43    real(r_kind),allocatable  :: data_all(:)
  44
  45    real,allocatable          :: sat_zen(:,:)
  46    real,allocatable          :: lon(:,:),lat(:,:)
  47
  48    real,allocatable          :: ch_tmp(:,:,:)
  49    integer                   :: bands
  50    integer                   :: crdate,crtime,dd
  51
  52    real,allocatable          :: bt(:,:,:)
  53    integer(i_kind)           :: doy,mlen(12),mday(12),mon,mm
  54
  55    integer(i_kind)  :: iscan
  56    integer(i_kind)  :: num_goesimg_file
  57    integer(i_kind)  :: num_goesimg_local, num_goesimg_global, num_goesimg_used, num_goesimg_thinned
  58    integer(i_kind)  :: num_goesimg_used_tmp
  59    integer(i_kind)  :: i, j, l,nchannel,ichan, iskip, ngoes
  60    integer(i_kind)  :: itx, k, nele, itt, n
  61    integer(i_kind)  :: error_status
  62    character(20)    ::filename
  63    real, allocatable :: in(:), out(:)
  64
  65 ! Set standard parameters
  66    integer          :: start(1:20)       ! Array of dimension starts.
  67
  68 ! use start for ir sensors.
  69    character(len=512)    :: input_file
  70    integer               :: length,ndims     ! Filename length.
  71    integer               :: rcode,io_status,ii,jj           ! NETCDF return code.
  72    integer               :: cdfid,varid,ivtype              ! NETCDF file IDs.
  73    character(len=20)     :: var2d(6)
  74    integer, dimension(3) :: dims ,dimids
  75    integer               :: natts
  76
  77 !------------------convert GVAR to BT-----------------------------
  78    real             :: tff
  79    real,parameter   :: c1 = 1.191066e-05, c2 = 1.438833
  80    real             :: m(4),b(4),va(4),vb(4),aa(4),ab(4), &
  81                        ba(4),bb(4),ra(4),rb(4)
  82    real             :: va_14(4),vb_14(4),aa_14(4),ab_14(4), &
  83                        ba_14(4),bb_14(4),ra_14(4),rb_14(4)
  84    real             :: va_13(4),vb_13(4),aa_13(4),ab_13(4), &
  85                        ba_13(4),bb_13(4),ra_13(4),rb_13(4)
  86    real             :: va_15(4),vb_15(4),aa_15(4),ab_15(4), &
  87                        ba_15(4),bb_15(4),ra_15(4),rb_15(4)
  88    data  m  /  227.3889,   38.8383,     5.2285,     5.5297/
  89    data  b  /   68.2167,   29.1287,    15.6854,    16.5892/
  90
  91 !-------------------goes-14 calibration coefficient----------------------------------
  92    data va_14 / 2577.3518, 1519.3488,  933.98541,  752.88143/
  93    data vb_14 / 2577.3518, 1518.5610,  934.19579,  752.82392/
  94    data aa_14  /-1.5565294,-3.9656363,-0.50944128,-0.16549136/
  95    data ab_14  /-1.5565294,-3.9615475,-0.51352435,-0.16396181/
  96    data ba_14  / 1.0027731, 1.0133248,  1.0029188,  1.0001953/
  97    data bb_14  / 1.0027731, 1.0135737,  1.0027813,  1.0002291/
  98    data ra_14  /-4.0683469e-07,-7.5834376e-06,-3.3213255e-06,9.0998038e-07/
  99    data rb_14  /-4.0683469e-07,-7.9139638e-06,-2.9825801e-06,8.1000947e-07/
 100 !-------------------goes-13 calibration coefficient----------------------------------
 101    data va_13  / 2561.7421, 1522.5182, 937.23449, 751.92589/
 102    data vb_13  / 2561.7421, 1521.6645, 937.27498, 751.92589/
 103    data aa_13  /-1.4755462, -4.1556932, -0.52227011, -0.16048355/
 104    data ab_13  /-1.4755462, -4.1411143, -0.51783545, -0.16048355/
 105    data ba_13  / 1.0025686, 1.0142082, 1.0023802, 1.0006854/
 106    data bb_13  / 1.0028656, 1.0142255, 1.0023789, 1.0006854/
 107    data ra_13  /-5.8203946e-07, -8.0255086e-06, -2.0798856e-06, -3.9399190e-07/
 108    data rb_13  /-5.8203946e-07, -8.0755893e-06, -2.1027609e-06, -3.9399190e-07/
 109 !-------------------goes-15 calibration coefficient----------------------------------
 110    data va_15  / 2562.7905, 1521.1988, 935.89417, 753.72229/
 111    data vb_15  / 2562.7905, 1521.5277, 935.78158, 753.93403/
 112    data aa_15  /-1.5903939, -3.9614910, -0.51798741, -0.16676232/
 113    data ab_15  /-1.5903939, -3.9580619, -0.51371643, -0.17147513/
 114    data ba_15  / 1.0026700, 1.0132094, 1.0025141, 1.0002673/
 115    data bb_15  / 1.0026700, 1.0130975, 1.0025320, 1.0001237/
 116    data ra_15  /-3.1926330e-07, -7.4310172e-06, -2.3893260e-06, 7.3287547e-07/
 117    data rb_15  /-3.1926330e-07, -7.3039584e-06, -2.4517009e-06, 1.1424349e-06/
 118
 119    data  mlen /31,28,31,30,31,30,31,31,30,31,30,31/ ! days of each month
 120
 121    if (trace_use) call da_trace_entry("da_read_obs_ncgoesimg")
 122
 123 !  0.0  Initialize variables
 124 !-----------------------------------
 125    platform_id  = 4   ! GOES series
 126    sensor_id    = 22  ! imager
 127    ifgat        = 1
 128    nchannel     = 4   ! # of channels
 129
 130    if(infile(6:7)=='13') then
 131       satellite_id = 13
 132    elseif(infile(6:7)=='14') then
 133       satellite_id = 14
 134    elseif(infile(6:7)=='15') then
 135       satellite_id = 15
 136    else
 137       write(*,*) 'goes satellite ', satellite_id, ' is not supported'
 138       return
 139    endif
 140
 141 ! determine if sensor triplet is in the sensor list
 142 !-----------------------------------------------------
 143    inst = 0
 144    do ngoes = 1, rtminit_nsensor
 145       if (platform_id  == rtminit_platform(ngoes) &
 146          .and. sensor_id == rtminit_sensor(ngoes) &
 147          .and. satellite_id == rtminit_satid(ngoes)) then
 148          inst = ngoes
 149       else
 150          cycle
 151       end if
 152    end do
 153    if (inst == 0) then
 154       write(unit=message(1),fmt='(A,I2,A)') " goes-",satellite_id,"-imager is not in sensor list"
 155       call da_warning(__FILE__,__LINE__, message(1:1))
 156       return
 157    end if
 158
 159 ! assign corresponding calibration coeffs
 160 !-------------------------------------------
 161    if (satellite_id==13) then  !only use the first detector's parameters
 162       va=va_13
 163       aa=aa_13
 164       ba=ba_13
 165       ra=ra_13
 166    end if
 167    if (satellite_id==14) then
 168       va=va_14
 169       aa=aa_14
 170       ba=ba_14
 171       ra=ra_14
 172    end if
 173    if (satellite_id==15) then
 174       va=va_15
 175       aa=aa_15
 176       ba=ba_15
 177       ra=ra_15
 178    end if
 179
 180    allocate(ptotal(0:num_fgat_time))
 181    ptotal(0:num_fgat_time) = 0
 182    iobs = 0                 ! for thinning, argument is inout
 183
 184 ! 1.0 Read NC files: sepeate file for different channels
 185 !-------------------------------------------------------
 186
 187 ! variables to read from NC file
 188 !--------------------------------
 189       var2d(1)='imageDate'
 190       var2d(2)='imageTime'
 191       var2d(3)='lat'
 192       var2d(4)='lon'
 193       var2d(5)='bands'
 194       var2d(6)='data'
 195
 196 ! determine file existant and open file
 197 !----------------------------------------
 198       write(input_file,fmt='(A,A,I2.2,A)') trim(infile),'-',01,'.nc'
 199       write(*,*) "Reading dimension of ", trim(input_file)
 200       inquire(file=trim(input_file), exist=isfile)
 201       if ( isfile ) then
 202          length = len_trim(input_file)
 203          rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid)
 204       else
 205          write(unit=message(1),fmt='(A,I2,A)') " goes-",satellite_id,"-imager-01.nc not exist"
 206          call da_warning(__FILE__,__LINE__, message(1:1))
 207          return
 208       end if
 209
 210          io_status = nf_inq_varid(cdfid, trim(var2d(6)), varid)
 211          rcode     = nf_inq_var  (cdfid, varid, trim(var2d(6)), ivtype, ndims, dimids, natts)
 212
 213          do i = 1, ndims
 214             rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
 215          end do
 216
 217          rcode = nf_close(cdfid)
 218
 219          write(*,*) "ndims dims for lat lon ",ndims, dims
 220
 221          allocate(  bt(nchan,dims(1), dims(2)))
 222          allocate(       lat(dims(1), dims(2)))
 223          allocate(       lon(dims(1), dims(2)))
 224          allocate(   sat_zen(dims(1), dims(2)))
 225          allocate(  ch_tmp(dims(1), dims(2), dims(3)) )
 226
 227    ncfile_loop:  do ichan=1,nchannel  ! Channel loop
 228 !-------------------------------------------------------
 229
 230       num_goesimg_file    = 0
 231       num_goesimg_local   = 0
 232       num_goesimg_global  = 0
 233       num_goesimg_used    = 0
 234       num_goesimg_thinned = 0
 235
 236 ! determine file existant and open file
 237 !----------------------------------------
 238       write(input_file,fmt='(A,A,I2.2,A)') trim(infile),'-',ichan,'.nc'
 239       write(*,*) "Reading ", trim(input_file)
 240       inquire(file=trim(input_file), exist=isfile)
 241       if ( isfile ) then
 242          length = len_trim(input_file)
 243          rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid)
 244       end if
 245
 246          do n=1,4  ! read the first 4 variables
 247         !----------------------------------------
 248             start     = 1
 249             io_status = nf_inq_varid(cdfid, trim(var2d(n)), varid)
 250             if (io_status /= 0 ) then
 251                write(unit=message(1),fmt='(a,a,a)') "Variable ",trim(var2d(n)), " does not exist"
 252                call da_warning(__FILE__,__LINE__, message(1:1))
 253                cycle
 254             endif
 255             rcode = nf_inq_var( cdfid, varid, trim(var2d(n)), ivtype, ndims, dimids, natts )
 256
 257             do i = 1, ndims
 258                rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
 259             end do
 260
 261             select case(trim(var2d(n)))
 262             case ('imageDate')
 263                rcode = nf_get_vara_int( cdfid, varid, start, dims, crdate)
 264             case ('imageTime')
 265                rcode = nf_get_vara_int( cdfid, varid, start, dims, crtime)
 266             case ('lat')
 267                rcode = nf_get_vara_double( cdfid, varid, start, dims, lat)
 268             case ('lon')
 269                rcode = nf_get_vara_double( cdfid, varid, start, dims, lon)
 270             end select
 271          end do  ! end reading of the first 4 variables
 272         !--------------------------------------------------
 273
 274         ! determine ccyymmddhh
 275         !---------------------------
 276          idate5(1)= int(crdate/1000)
 277          dd = mod(crdate,1000)
 278          if ( mod(idate5(1),4) == 0 ) then
 279             mlen(2) = 29
 280             if ( MOD(idate5(1),100) == 0 ) then
 281                mlen(2) = 28
 282             endif
 283             if ( MOD(idate5(1),400) == 0 ) then
 284                mlen(2) = 29
 285             endif
 286          endif
 287          mm = 0
 288          do mon=1,12
 289             mday(mon) = mm
 290             mm = mm + mlen(mon)
 291          end do
 292
 293          do mon= 2,12
 294             if (dd.gt.mday(mon-1).and.dd.le.mday(mon))then
 295                idate5(2) = mon-1
 296                idate5(3) = dd - mday(mon-1)
 297             end if
 298          end do
 299
 300          idate5(4)=int(crtime/10000)
 301          idate5(5)=int(crtime/100)-int(crtime/10000)*100
 302          idate5(6)=crtime-int(crtime/100)*100
 303
 304 !------------- calculate satelize zenith angle -----------------------
 305          sat_zen=missing_r
 306          do jj=1,dims(2)
 307          do ii=1,dims(1)
 308             call da_get_satzen(lat(ii,jj),lon(ii,jj),satellite_id,sat_zen(ii,jj))
 309             if(sat_zen(ii,jj) > 75.0) sat_zen(ii,jj)=missing_r
 310          end do
 311          end do
 312
 313 !---------------------------------------------------------------------------
 314       do n=5,6  ! read bands/data
 315
 316          io_status = nf_inq_varid(cdfid, trim(var2d(n)), varid)
 317          if (io_status /= 0 ) then
 318             write(unit=message(1),fmt='(a,a,a)') "Variable ",trim(var2d(n)), " does not exist"
 319             call da_warning(__FILE__,__LINE__, message(1:1))
 320             cycle
 321          endif
 322          rcode = nf_inq_var( cdfid, varid, trim(var2d(n)), ivtype, ndims, dimids, natts )
 323          dims = 0
 324
 325          do i = 1, ndims
 326             rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
 327          end do
 328
 329          select case(trim(var2d(n)))
 330          case ('data')
 331             rcode = nf_get_vara_double( cdfid, varid, start, dims, ch_tmp)
 332          case ('bands')
 333             rcode = nf_get_vara_int2( cdfid, varid, start, dims, bands)
 334          end select
 335       end do
 336
 337       if (bands==6) bands=5  ! crtm coeffs index?
 338
 339 !------------------------convert to BT-----------------------
 340
 341       do jj=1,dims(2)
 342       do ii=1,dims(1)
 343          ch_tmp(ii,jj,1) = ch_tmp(ii,jj,1)/32.0  !16 bit to 10 bit
 344          tff  = (ch_tmp(ii,jj,1)-b(bands-1))/m(bands-1)
 345
 346          if ( (tff .gt. 0.0) .and. (lat(ii,jj).lt.100.0) ) then
 347             tff = c2*va(bands-1)/alog(1+c1*va(bands-1)**3/tff)
 348             tff = aa(bands-1)+ba(bands-1)*tff+ra(bands-1)*tff**2
 349             bt(ichan,ii,jj) = tff
 350          else
 351             bt(ichan,ii,jj) = missing_r
 352          end if
 353       end do
 354       end do
 355
 356       rcode = nf_close(cdfid)
 357
 358    end do ncfile_loop ! end channel loop
 359 !---------------------------------
 360
 361    allocate (head)
 362    nullify  (head % next )
 363    p => head
 364
 365    allocate(data_all(nchan))
 366
 367 ! 2.0 loop to read data file
 368 !----------------------------------------------
 369    do jj=1,dims(2)
 370    do ii=1,dims(1)
 371
 372          num_goesimg_file = num_goesimg_file + 1
 373          info%lon  =  lon(ii,jj)  ! longitude
 374          info%lat  =  lat(ii,jj)  ! latitude
 375          call da_llxy (info, loc, outside, outside_all)
 376          if (outside_all) cycle
 377          if (inst == 0) cycle
 378          if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
 379              idate5(2) < 1    .or. idate5(2) >   12 .or. &
 380              idate5(3) < 1    .or. idate5(3) >   31 .or. &
 381              idate5(4) <0     .or. idate5(4) >   24 .or. &
 382              idate5(5) <0     .or. idate5(5) >   60 )then
 383             write(6,*)'READ_goeimg:  ### ERROR IN READING ', 'goeimg', ' BUFR DATA:', &
 384                 ' STRANGE OBS TIME (YMDHM):', idate5(1:5)
 385             exit
 386          end if
 387          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 388          if ( obs_time < time_slots(0) .or.  &
 389             obs_time >= time_slots(num_fgat_time) ) cycle
 390          do ifgat=1,num_fgat_time
 391             if ( obs_time >= time_slots(ifgat-1) .and.  &
 392                   obs_time  < time_slots(ifgat) ) exit
 393          end do
 394          num_goesimg_global = num_goesimg_global + 1
 395          ptotal(ifgat) = ptotal(ifgat) + 1
 396          if (outside)  cycle
 397          num_goesimg_local = num_goesimg_local + 1
 398          write(unit=info%date_char, &
 399             fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 400             idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 401             ':', idate5(5), ':', idate5(6)
 402             info%elv = 0.0  !aquaspot%selv
 403
 404 ! 3.0 Map obs to thinning grid
 405 !-------------------------------------------------------------------
 406          if (thinning) then
 407             dlat_earth = info%lat
 408             dlon_earth = info%lon
 409             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 410             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 411             dlat_earth = dlat_earth*deg2rad
 412             dlon_earth = dlon_earth*deg2rad
 413             timedif = 0.
 414             crit = 1.
 415             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 416             if (.not. iuse) then
 417                num_goesimg_thinned=num_goesimg_thinned+1
 418                cycle
 419             end if
 420          end if
 421
 422          do l=1,nchan
 423             data_all(l)=bt(l,ii,jj)
 424          end do
 425
 426          num_goesimg_used = num_goesimg_used + 1
 427
 428 !  4.0   assign information to sequential radiance structure
 429 !--------------------------------------------------------------------------
 430          allocate ( p % tb_inv (1:nchan) )
 431          p%info                  = info
 432          p%loc                   = loc
 433          p%landsea_mask          = 1  ! ???
 434          p%scanpos               = ii ! ??? "scan" position
 435          p%satzen                = sat_zen(ii,jj)
 436          p%solzen                = 0.0
 437          p%tb_inv(1:nchan)       = data_all(1:nchan)
 438          p%sensor_index          = inst
 439          p%ifgat                 = ifgat
 440          allocate (p%next)   ! add next data
 441          p => p%next
 442          nullify (p%next)
 443
 444    end do
 445    end do
 446
 447    deallocate(data_all) ! Deallocate data arrays
 448    deallocate(bt)
 449    deallocate(sat_zen)
 450    deallocate(lat)
 451    deallocate(lon)
 452    deallocate(ch_tmp)
 453
 454 !------------------------------------------------------
 455
 456    if (thinning .and. num_goesimg_global > 0 ) then
 457 #ifdef DM_PARALLEL
 458
 459       ! Get minimum crit and associated processor index.
 460       j = 0
 461       do ifgat = 1, num_fgat_time
 462             j = j + thinning_grid(inst,ifgat)%itxmax
 463       end do
 464
 465
 466       allocate ( in  (j) )
 467       allocate ( out (j) )
 468       j = 0
 469       do ifgat = 1, num_fgat_time
 470             do i = 1, thinning_grid(inst,ifgat)%itxmax
 471                j = j + 1
 472                in(j) = thinning_grid(inst,ifgat)%score_crit(i)
 473             end do
 474       end do
 475
 476       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 477
 478       call wrf_dm_bcast_real (out, j)
 479
 480       j = 0
 481       do ifgat = 1, num_fgat_time
 482             do i = 1, thinning_grid(inst,ifgat)%itxmax
 483                j = j + 1
 484                if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(inst,ifgat)%ibest_obs(i) = 0
 485             end do
 486       end do
 487       deallocate( in  )
 488       deallocate( out )
 489
 490 #endif
 491       ! Delete the nodes which being thinning out
 492       p => head
 493       prev => head
 494       head_found = .false.
 495       num_goesimg_used_tmp = num_goesimg_used
 496
 497       do j = 1, num_goesimg_used_tmp
 498          n = p%sensor_index
 499          ifgat = p%ifgat
 500          found = .false.
 501
 502          do i = 1, thinning_grid(n,ifgat)%itxmax
 503             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 504                found = .true.
 505                exit
 506             end if
 507          end do
 508
 509          ! free current data
 510          if ( .not. found ) then
 511             current => p
 512             p => p%next
 513             if ( head_found ) then
 514                prev%next => p
 515             else
 516                head => p
 517                prev => p
 518             end if
 519             deallocate ( current % tb_inv )
 520       !      deallocate ( current % cloud_flag )
 521             deallocate ( current )
 522             num_goesimg_thinned = num_goesimg_thinned + 1
 523             num_goesimg_used = num_goesimg_used - 1
 524             continue
 525          end if
 526
 527          if ( found .and. head_found ) then
 528             prev => p
 529             p => p%next
 530             continue
 531          end if
 532          if ( found .and. .not. head_found ) then
 533             head_found = .true.
 534             head => p
 535             prev => p
 536             p => p%next
 537          end if
 538
 539       end do
 540
 541    end if  ! End of thinning
 542 !stop
 543    iv%total_rad_pixel   = iv%total_rad_pixel + num_goesimg_used
 544    iv%total_rad_channel = iv%total_rad_channel + num_goesimg_used*nchan
 545
 546    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_goesimg_used
 547    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_goesimg_global
 548
 549    do i = 1, num_fgat_time
 550       ptotal(i) = ptotal(i) + ptotal(i-1)
 551       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 552    end do
 553    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 554       write(unit=message(1),fmt='(A,I10,A,I10)') &
 555           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 556       call da_warning(__FILE__,__LINE__,message(1:1))
 557    endif
 558
 559    write(unit=stdout,fmt='(a)') 'num_goesimg_file, num_goesimg_global, num_goesimg_local, num_goesimg_used, num_goesimg_thinned'
 560    write(stdout,*) num_goesimg_file, num_goesimg_global, num_goesimg_local, num_goesimg_used, num_goesimg_thinned
 561
 562
 563    !  5.0 allocate innovation radiance structure
 564    !----------------------------------------------------------------
 565
 566
 567    if (num_goesimg_used > 0) then
 568       iv%instid(inst)%num_rad  = num_goesimg_used
 569       iv%instid(inst)%info%nlocal = num_goesimg_used
 570       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 571         'Allocating space for radiance innov structure', &
 572          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 573       call da_allocate_rad_iv (inst, nchan, iv)
 574    end if
 575
 576    !  6.0 assign sequential structure to innovation structure
 577    !-------------------------------------------------------------
 578    p => head
 579    do n = 1, num_goesimg_used
 580       i = p%sensor_index
 581       call da_initialize_rad_iv (i, n, iv, p)
 582       current => p
 583       p => p%next
 584
 585       ! free current data
 586       deallocate ( current % tb_inv )
 587 !      deallocate ( current % cloud_flag )
 588       deallocate ( current )
 589    end do
 590    deallocate ( p )
 591    deallocate (ptotal)
 592
 593    if (trace_use) call da_trace_exit("da_read_obs_ncgoesimg")
 594
 595 end subroutine da_read_obs_ncgoesimg