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