var/da/da_radiance/da_radiance_init.inc

   1 subroutine da_radiance_init(iv,ob)
   2 !------------------------------------------------------------------------------
   3 !  PURPOSE: subroutine to initialize radiances.
   4 !
   5 !  METHOD:
   6 !  1.0 Set up from namelist parameter
   7 !  2.0 Set up some common variables for innovation/observation
   8 !  3.0 Initialize RTTOV / CRTM
   9 !  4.0 Set up bias correction
  10 !  5.0 Read error factor file
  11 !  6.0 Get FGAT time slots
  12 !
  13 !  HISTORY: 10/24/2007 Created from da_crtm_init            Tom Auligne
  14 !  HISTORY: 12/15/2008 getting FGAT time slots is moved to
  15 !                      da_setup_obs_structures.inc.         Hui-Chuan Lin
  16 !------------------------------------------------------------------------------
  17
  18  implicit none
  19
  20  type (iv_type), intent (inout) :: iv
  21  type (y_type) , intent (inout) :: ob
  22
  23 !
  24 !  local arguments
  25 !-------------------
  26  integer   :: n, j, ichan, iret
  27  integer :: nsensor, unit_factor_rad
  28  integer     :: error
  29  integer, allocatable   ::  nscan(:), nchanl(:)
  30
  31 ! local variables
  32 !----------------
  33  integer             :: idum, wmo_sensor_id, sensor_type, iost
  34  integer             :: iunit
  35  character(len=filename_len)  :: filename
  36  character(len=20)   :: cdum
  37  character(len=12)   :: cdum12
  38  real                :: error_cld
  39
  40 ! local variables for tuning error factor
  41 !----------------------------------------
  42  character(len=20)   ::  rttovid_string
  43  integer             ::  num_tot
  44  real                ::  joa, jo, trace, factor
  45
  46   call da_trace_entry("da_radiance_init")
  47
  48 !--------------------------------------------------------------
  49 !  1.0 setup from namelist parameter
  50 !--------------------------------------------------------------
  51   nsensor = rtminit_nsensor
  52   allocate (nscan(nsensor))
  53   allocate (nchanl(nsensor))
  54
  55 !----------------------------------------------------------------
  56 !  2.0 set up some common variables for innovation/observation structure
  57 !----------------------------------------------------------------
  58   iv % num_inst = nsensor
  59   ob % num_inst = nsensor
  60
  61   allocate (iv%instid(1:nsensor))
  62   allocate (ob%instid(1:nsensor))
  63   allocate (satinfo(1:nsensor))
  64
  65   iv%instid(1:nsensor)%num_rad = 0
  66   ob%instid(1:nsensor)%num_rad = 0
  67
  68   loop_sensor: do n = 1, nsensor
  69
  70    iv%instid(n)%platform_id  = rtminit_platform(n)
  71    iv%instid(n)%satellite_id = rtminit_satid(n)
  72    iv%instid(n)%sensor_id    = rtminit_sensor(n)
  73    if ( rtminit_satid(n) < 10 ) then
  74       write(iv%instid(n)%rttovid_string, '(a,i1,a)')  &
  75              trim( rttov_platform_name(rtminit_platform(n)) )//'-',  &
  76              rtminit_satid(n),     &
  77              '-'//trim( rttov_inst_name(rtminit_sensor(n)) )
  78       write(iv%instid(n)%rttovid_string_coef, '(a,i1,a)')  &
  79              trim( rttov_platform_name(rtminit_platform(n)) )//'_',  &
  80              rtminit_satid(n),     &
  81              '_'//trim( rttov_inst_name(rtminit_sensor(n)) )
  82    else
  83       write(iv%instid(n)%rttovid_string, '(a,i2.2,a)')  &
  84              trim( rttov_platform_name(rtminit_platform(n)) )//'-',  &
  85              rtminit_satid(n),     &
  86              '-'//trim( rttov_inst_name(rtminit_sensor(n)) )
  87       write(iv%instid(n)%rttovid_string_coef, '(a,i2.2,a)')  &
  88              trim( rttov_platform_name(rtminit_platform(n)) )//'_',  &
  89              rtminit_satid(n),     &
  90              '_'//trim( rttov_inst_name(rtminit_sensor(n)) )
  91    end if
  92
  93    if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'msu' ) then
  94       nchanl(n)  = 4
  95       nscan(n)   = 11
  96    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'hirs' ) then
  97       nchanl(n)  = 19
  98       nscan(n)   = 56
  99    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'amsua' ) then
 100       nchanl(n)  = 15
 101       nscan(n)   = 30
 102    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'amsub' ) then
 103       nchanl(n)  = 5
 104       nscan(n)   = 90
 105    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'airs' ) then
 106       nchanl(n)  = 281
 107       nscan(n)   = 90
 108    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'hsb' ) then
 109       nchanl(n)  = 4
 110       nscan(n)   = 90
 111    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'mhs' ) then
 112       nchanl(n)  = 5
 113       nscan(n)   = 90
 114    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'ssmis' ) then
 115       nchanl(n)  = 24
 116       nscan(n)   = 60
 117    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'mwts' ) then
 118       nchanl(n)  = 4
 119       nscan(n)   = 15
 120    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'mwhs' ) then
 121       nchanl(n)  = 5
 122       nscan(n)   = 98
 123    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'mwhs2' ) then
 124       nchanl(n)  = 15
 125       nscan(n)   = 98
 126    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'atms' ) then
 127       nchanl(n)  = 22
 128       nscan(n)   = 96
 129    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'iasi' ) then
 130      nchanl(n)  = 616
 131       nscan(n)   = 60
 132    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'seviri' ) then
 133      nchanl(n)  =  8
 134      nscan(n)   = 90
 135    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'amsr2' ) then
 136       nchanl(n)  = 14
 137       nscan(n)   = 486
 138    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'ahi' ) then
 139      !open the ahi info file
 140       open(unit=1990,file='ahi_info',status='old',iostat=iret)
 141       if(iret /= 0)then
 142          call da_error(__FILE__,__LINE__,(/"Read ahi_info error: no such file"/))
 143       end if
 144      !read ahi information
 145       read(1990,*)
 146       read(1990,*)
 147       read(1990,*)
 148       read(1990,*) nscan(n)
 149       close(1990)
 150       write(*,*) nscan(n)
 151       nchanl(n)  = 10
 152    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'imgr' ) then
 153       nchanl(n)  = 4
 154       nscan(n)   = 60
 155    else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'gmi' ) then
 156       nchanl(n)  = 13
 157       nscan(n)   = 221
 158    else
 159       write(unit=message(1),fmt='(A)') "Unrecognized instrument: "
 160       write(unit=message(2),fmt='(A,I4)') "rtminit_platform = ",rtminit_platform(n)
 161       write(unit=message(3),fmt='(A,I4)') "rtminit_satid    = ",rtminit_satid(n)
 162       write(unit=message(4),fmt='(A,I4)') "rtminit_sensor   = ",rtminit_sensor(n)
 163       write(unit=message(5),fmt='(A)') "Check your namelist settings"
 164       call da_error(__FILE__,__LINE__,message(1:5))
 165    end if
 166
 167    iv%instid(n)%nchan  = nchanl(n)
 168    ob%instid(n)%nchan  = nchanl(n)
 169
 170    allocate ( iv%instid(n)%ichan(1:nchanl(n)), stat = error )
 171    if( error /= 0 ) then
 172       call da_error(__FILE__,__LINE__, &
 173          (/"Memory allocation error to iv%instid(n)%ichan"/))
 174    end if
 175
 176    allocate ( ob%instid(n)%ichan(1:nchanl(n)), stat = error )
 177    if( error /= 0 ) then
 178       call da_error(__FILE__,__LINE__, &
 179          (/"Memory allocation error to ob%instid(n)%ichan"/))
 180    end if
 181
 182    call da_get_unit(iunit)
 183    filename='radiance_info/'//trim(adjustl(iv%instid(n)%rttovid_string))//'.info'
 184    open(unit=iunit,file=filename, form='formatted',iostat = iost, status='old')
 185
 186    if (iost /= 0) then
 187       message(1)="Cannot open radiance info file "//adjustl(filename)
 188       call da_error(__FILE__,__LINE__,message(1:1))
 189    end if
 190
 191    allocate ( satinfo(n) % ichan(nchanl(n)) )
 192    allocate ( satinfo(n) % iuse (nchanl(n)) )
 193    allocate ( satinfo(n) % error(nchanl(n)) )
 194    allocate ( satinfo(n) % error_cld(nchanl(n)) )
 195    allocate ( satinfo(n) % polar(nchanl(n)) )
 196
 197    satinfo(n) % error_cld(:) = 500.0 !initialize
 198
 199    ! Allocate additional fields for AHI
 200    if ( index(iv%instid(n)%rttovid_string, 'ahi') > 0 ) then
 201       allocate ( satinfo(n) % BTLim(nchanl(n)) )
 202       allocate ( satinfo(n) % ca1(nchanl(n)) )
 203       allocate ( satinfo(n) % ca2(nchanl(n)) )
 204       allocate ( satinfo(n) % clearSkyBias(nchanl(n)) )
 205    endif
 206
 207    read(iunit,*)
 208    do j = 1, nchanl(n)
 209       read(iunit,'(1x,5i5,2e18.10,a20)') &
 210                      wmo_sensor_id, &
 211                satinfo(n)%ichan(j), &
 212                        sensor_type, &
 213                satinfo(n)%iuse(j) , &
 214                               idum, &
 215                satinfo(n)%error(j), &
 216                satinfo(n)%polar(j), &
 217                cdum
 218       !in the current radiance info files, the last column
 219       !can be either sensor_id_string or blank
 220       if ( len_trim(cdum) > 0 .and. index(cdum,'-') == 0 ) then
 221          ! read the line again to get error_cld when it is available
 222          backspace(iunit)
 223          read(iunit,'(1x,5i5,2e18.10,f10.5)')  &
 224                      wmo_sensor_id, &
 225                satinfo(n)%ichan(j), &
 226                        sensor_type, &
 227                satinfo(n)%iuse(j) , &
 228                               idum, &
 229                satinfo(n)%error(j), &
 230                satinfo(n)%polar(j), &
 231                error_cld
 232          if ( error_cld > 0.0 ) then
 233             satinfo(n)%error_cld(j) = error_cld
 234          end if
 235       end if
 236
 237       ! If AHI, read some extra things
 238       ! Unfortunately, we need to read everything again...
 239       if ( index(iv%instid(n)%rttovid_string, 'ahi') > 0 ) then
 240          backspace(iunit)
 241          read(iunit,'(1x,5i5,2e18.10,a12,f8.2,2f6.2,f9.3)') &
 242                      wmo_sensor_id, &
 243                satinfo(n)%ichan(j), &
 244                        sensor_type, &
 245                satinfo(n)%iuse(j) , &
 246                               idum, &
 247                satinfo(n)%error(j), &
 248                satinfo(n)%polar(j), &
 249                cdum12,              &
 250                satinfo(n)%BTLim(j), &
 251                satinfo(n)%ca1(j)  , &
 252                satinfo(n)%ca2(j)  , &
 253                satinfo(n)%clearSkyBias(j)
 254          if ( j == 1 ) then
 255             write(*,*)'Reading extra data for AHI'
 256             write(*,*)'Channel       BTLim      ca1     ca2     clearSkyBias'
 257          endif
 258          write(*,fmt='(i7,6x,4f9.3)') satinfo(n)%ichan(j), satinfo(n)%BTLim(j), satinfo(n)%ca1(j), satinfo(n)%ca2(j), satinfo(n)%clearSkyBias(j)
 259       endif
 260
 261      iv%instid(n)%ichan(j) = satinfo(n)%ichan(j)
 262      ob%instid(n)%ichan(j) = satinfo(n)%ichan(j)
 263    end do
 264    call da_free_unit(iunit)
 265
 266    if ( use_blacklist_rad ) then
 267       call da_blacklist_rad(trim(rttov_platform_name(rtminit_platform(n))), &
 268                             rtminit_satid(n), &
 269                             trim(rttov_inst_name(rtminit_sensor(n))), &
 270                             nchanl(n), &
 271                             satinfo(n)%iuse )
 272    end if
 273
 274   end do loop_sensor
 275
 276 !---------------------------------------------------------------------
 277 ! 3.0 Interface to the initialization subroutine of RTTOV and CRTM
 278 !---------------------------------------------------------------------
 279
 280     if (rtm_option == rtm_option_rttov) then
 281 #ifdef RTTOV
 282        call da_rttov_init(iv,ob,nsensor,nchanl)
 283 #else
 284        call da_error(__FILE__,__LINE__, &
 285           (/"Must compile with $RTTOV option for radiances"/))
 286 #endif
 287     end if
 288
 289     if (rtm_option == rtm_option_crtm) then
 290 #ifdef CRTM
 291        call da_crtm_init(iv,ob,nsensor)
 292 #else
 293        call da_error(__FILE__,__LINE__, &
 294           (/"Must compile with $CRTM option for radiances"/))
 295 #endif
 296     end if
 297
 298 !-------------------------------------------------------
 299 !  4.0 read bias correction coefs files
 300 !-------------------------------------------------------
 301
 302  loop_sensor2: do n = 1, nsensor
 303
 304    allocate ( satinfo(n) % scanbias  (nchanl(n),nscan(n)) )
 305    allocate ( satinfo(n) % scanbias_b(nchanl(n),nscan(n),18) )
 306    allocate ( satinfo(n) % bcoef     (nchanl(n),4) )
 307    allocate ( satinfo(n) % bcoef0    (nchanl(n)) )
 308    allocate ( satinfo(n) % error_std (nchanl(n)) )
 309
 310    satinfo(n) % error_std(:) = 500.0
 311    satinfo(n) % scanbias(:,:) = 0.0
 312    satinfo(n) % scanbias_b(:,:,:) = 0.0
 313    satinfo(n) % bcoef(:,:) = 0.0
 314    satinfo(n) % bcoef0(:) = 0.0
 315
 316   if (read_biascoef) then
 317   !  new bias coefs files
 318   !  use o-b standard deviation statistics from Harris and Kelly method as obs errors
 319   !----------------------------------
 320
 321     if ( index(iv%instid(n)%rttovid_string,'eos')  > 0 ) cycle   ! not implemented
 322     if ( index(iv%instid(n)%rttovid_string,'hirs') > 0 ) cycle   ! not implemented
 323
 324     call da_read_biascoef(iv%instid(n)%rttovid_string, &
 325                       nchanl(n),nscan(n),18,4,global, &
 326                       satinfo(n)%scanbias, &
 327                       satinfo(n)%scanbias_b, &
 328                       satinfo(n)%bcoef, &
 329                       satinfo(n)%bcoef0, &
 330                       satinfo(n)%error_std)
 331   else
 332     ! use values specified in radiance_info files as obs errors
 333     satinfo(n)%error_std = satinfo(n)%error
 334   end if
 335  end do loop_sensor2
 336
 337 !-------------------------------------------------------
 338 !  5.0 read error factor file
 339 !-------------------------------------------------------
 340  if (use_error_factor_rad) then
 341
 342     do n=1, rtminit_nsensor
 343        allocate ( satinfo(n)%error_factor(1:nchanl(n)) )
 344        satinfo(n)%error_factor(:) = 1.0
 345     end do
 346
 347     call da_get_unit(unit_factor_rad)
 348     open(unit_factor_rad, file='radiance_error.factor', &
 349          form='formatted',iostat = iost, status='old')
 350
 351     if (iost /= 0) then
 352        call da_error(__FILE__,__LINE__, &
 353          (/"Cannot open radiance error factor file: radiance_error.factor"/))
 354     end if
 355
 356     read(unit_factor_rad, *)
 357     do
 358       read(unit_factor_rad,fmt='(a15,i8,i8,3f15.5,f8.3)',iostat=iost)   &
 359           rttovid_string, ichan, num_tot, joa,jo,trace,factor
 360       if ( iost == 0 ) then
 361         do n=1, rtminit_nsensor
 362           if ( index(rttovid_string,trim(iv%instid(n)%rttovid_string))>0 ) then
 363              satinfo(n)%error_factor(ichan) = factor
 364              write(6,'(a,i5,a,f10.3)') trim(rttovid_string)//' Channel ', ichan, '  Error Factor = ', factor
 365              exit
 366           end if
 367         end do
 368       else
 369          exit
 370       end if
 371     end do
 372     close(unit_factor_rad)
 373     call da_free_unit(unit_factor_rad)
 374
 375  end if
 376
 377   deallocate(nscan)
 378   deallocate(nchanl)
 379
 380   call da_trace_exit("da_radiance_init")
 381 end subroutine da_radiance_init