var/da/da_radiance/da_setup_radiance_structures.inc

   1 subroutine da_setup_radiance_structures( grid, ob, iv )
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Define, allocate and read of tovs raidance observation structure.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (domain) , intent(inout)   :: grid       ! model data
  10    type ( y_type), intent(inout)   :: ob         ! Observation structure.
  11    type (iv_type), intent(inout)   :: iv         ! O-B structure.
  12
  13    character(len=200)          :: filename
  14    integer                     :: i, j, n, ios, ifgat
  15    logical                     :: lprinttovs
  16
  17    ! crtm_cloud
  18    integer  :: n1,n2,k,its,ite,jts,jte,kts,kte,inst
  19    integer  :: data_format, iret ! AHI radiance related
  20    integer  :: ahi_info_unit
  21
  22    if (trace_use) call da_trace_entry("da_setup_radiance_structures")
  23
  24    !-------------------------------------------------------------------
  25    ! [1.0] Initialize RTTOV coefs and innovations vector for radiance
  26    !-------------------------------------------------------------------
  27
  28     call da_radiance_init(iv, ob)
  29
  30     do n = 1, rtminit_nsensor
  31        iv%instid(n)%rad_monitoring = rad_monitoring(n)
  32     enddo
  33
  34    !-------------------------------
  35    ! 1.1 Make thinning grids
  36    !------------------------------
  37
  38    if (thinning) then
  39       allocate(thinning_grid(iv%num_inst,num_fgat_time))
  40       do ifgat=1,num_fgat_time
  41          do n=1,iv%num_inst
  42             call makegrids (n,thinning_mesh(n),ifgat)
  43          end do
  44       end do
  45    end if
  46
  47    !-------------------------------------------------------------------
  48    ! [2.0] Read NCEP bufr tovs data in radiance innovations vector
  49    !-------------------------------------------------------------------
  50
  51    if ( (.not. use_filtered_rad) .and. (.not. use_pseudo_rad) .and. (.not. use_simulated_rad) ) then
  52
  53       if (use_hirs2obs) then
  54          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from hirs2.bufr'
  55          filename = 'hirs2 '
  56          call da_read_obs_bufrtovs ('hirs2', iv, filename)
  57       end if
  58
  59       if (use_msuobs) then
  60          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from msu.bufr'
  61          filename = 'msu'
  62          call da_read_obs_bufrtovs ('msu  ', iv, filename)
  63       end if
  64
  65       if (use_hirs3obs) then
  66          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from hirs3.bufr'
  67          filename = 'hirs3'
  68          call da_read_obs_bufrtovs('hirs3', iv, filename)
  69       end if
  70
  71       if (use_amsuaobs) then
  72          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from amsua.bufr'
  73          filename = 'amsua'
  74          call da_read_obs_bufrtovs ('amsua', iv, filename)
  75       end if
  76
  77       if (use_amsubobs) then
  78          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from amsub.bufr'
  79          filename = 'amsub'
  80          call da_read_obs_bufrtovs ('amsub', iv, filename)
  81       end if
  82
  83       if (use_hirs4obs) then
  84          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from hirs4.bufr'
  85          filename = 'hirs4'
  86          call da_read_obs_bufrtovs('hirs4', iv, filename)
  87       end if
  88
  89       if (use_mhsobs) then
  90          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from mhs.bufr'
  91          filename = 'mhs'
  92          call da_read_obs_bufrtovs('mhs  ', iv, filename)
  93       end if
  94
  95       if (use_mwtsobs) then
  96          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from mwtsa.dat and mwtsb.dat'
  97          filename = 'mwtsa'
  98          call da_read_obs_fy3('mwts ', iv, filename)
  99          filename = 'mwtsb'
 100          call da_read_obs_fy3('mwts ', iv, filename)
 101       end if
 102
 103       if (use_mwhsobs) then
 104          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from mwhsa.dat and mwhsb.dat'
 105          filename = 'mwhsa'
 106          call da_read_obs_fy3('mwhs ', iv, filename)
 107          filename = 'mwhsb'
 108          call da_read_obs_fy3('mwhs ', iv, filename)
 109       end if
 110
 111       if (use_mwhs2obs) then
 112          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from mwhs2.hdf'
 113          filename = 'mwhs2'
 114          call da_read_obs_hdf5mwhs2(iv, filename)
 115       end if
 116
 117       if (use_atmsobs) then
 118          write(unit=stdout,fmt='(a)') 'Reading radiance 1b data from atms.bufr'
 119          filename = 'atms'
 120          call da_read_obs_bufratms('atms ', iv, filename)
 121       end if
 122
 123       if (use_airsobs) then
 124          write(unit=stdout,fmt='(a)') 'Reading airs 1b data from airs.bufr'
 125          filename = 'airs'
 126          call da_read_obs_bufrairs ('airs     ',iv, filename)
 127       end if
 128
 129       if (use_eos_amsuaobs) then
 130          write(unit=stdout,fmt='(a)') 'Reading eos_amsua 1b data from airs.bufr'
 131          filename = 'airs'
 132          call da_read_obs_bufrairs ('eos_amsua',iv, filename)
 133       end if
 134
 135       if (use_hsbobs) then
 136          write(unit=stdout,fmt='(a)') 'Reading hsb 1b data from airs.bufr'
 137          filename = 'airs'
 138          call da_read_obs_bufrairs ('hsb      ',iv, filename)
 139       end if
 140
 141       if (use_ssmisobs) then
 142          write(unit=stdout,fmt='(a)') 'Reading ssmis data from ssmis.bufr'
 143          filename = 'ssmis'
 144          call da_read_obs_bufrssmis ('ssmis    ',iv, filename)
 145       end if
 146       if (use_iasiobs) then
 147          write(unit=stdout,fmt='(a)') 'Reading iasi data from iasi.bufr'
 148          filename = 'iasi'
 149          call da_read_obs_bufriasi ('iasi     ',iv, filename)
 150       end if
 151       if (use_seviriobs) then
 152          write(unit=stdout,fmt='(a)') 'Reading seviri data from seviri.bufr'
 153          filename = 'seviri'
 154          call da_read_obs_bufrseviri ('seviri   ',iv, filename)
 155       end if
 156       if (use_amsr2obs) then
 157 #if defined(HDF5)
 158          write(unit=stdout,fmt='(a)') 'Reading AMSR2 data in HDF5 format'
 159          call da_read_obs_hdf5amsr2 (iv, 'L1SGRTBR', 'L2SGCLWLD')
 160 #else
 161          message(1)='To read AMSR2 data, WRFDA must be compiled with HDF5'
 162          call da_error(__FILE__,__LINE__,message(1:1))
 163 #endif
 164       end if
 165       if (use_ahiobs) then
 166          call da_get_unit(ahi_info_unit)
 167          open(unit=ahi_info_unit,file='ahi_info',status='old',iostat=iret)
 168          if(iret /= 0)then
 169              call da_error(__FILE__,__LINE__,(/"Read ahi_info error: no such file"/))
 170          end if
 171
 172         !read ahi information
 173          read(ahi_info_unit,*)
 174          read(ahi_info_unit,*) data_format
 175          close(ahi_info_unit)
 176          call da_free_unit(ahi_info_unit)
 177
 178          if (data_format==1) then
 179 #if defined(HDF5)
 180             write(unit=stdout,fmt='(a)') 'Reading AHI data from cma hdf5'
 181             call da_read_obs_hdf5ahi (iv, 'L1AHITBR', 'L2AHICLP')
 182 #else
 183             call da_error(__FILE__,__LINE__,(/"To read AHI data, WRFDA must be compiled with HDF5"/))
 184 #endif
 185          end if
 186          if (data_format==2) then
 187             write(unit=stdout,fmt='(a)') 'Reading AHI data from geocat NETCDF4'
 188             call da_read_obs_netcdf4ahi_geocat (iv, 'L1AHITBR', 'L2AHICLP')
 189          end if
 190          if (data_format==3) then
 191             write(unit=stdout,fmt='(a)') 'Reading AHI data from JAXA NETCDF4'
 192 !            call da_read_obs_netcdf4ahi_jaxa (iv, 'L1AHITBR', 'L2AHICLP')
 193                         call da_read_obs_netcdf4ahi_jaxa (iv, 'L1AHITBR', 'L1AHITBRP')
 194          end if
 195          !if (data_format==4) then
 196             !filename = 'ahi'
 197             !call da_read_obs_bufrahi ('ahi   ',iv, filename)
 198          !end if
 199       end if
 200       if (use_goesimgobs) then
 201          write(unit=stdout,fmt='(a)') 'Reading netcdf goes imager GVAR data'
 202          !if (imager_format==1)then
 203          !   write(unit=stdout,fmt='(a)') 'Reading goesimg data from HDF4 files'
 204          !   filename = 'HDF4IMAGER-13'
 205          !   call da_read_obs_hdf4imager(iv, filename)
 206          !   filename = 'HDF4IMAGER-15'
 207          !   call da_read_obs_hdf4imager(iv, filename)
 208          !end if
 209          !if (imager_format==2)then
 210          !   write(unit=stdout,fmt='(a)') 'Reading goesimg data from NETCDF files'
 211             filename = 'goes-13-imager'
 212             call da_read_obs_ncgoesimg(iv, filename)
 213             filename = 'goes-14-imager'
 214             call da_read_obs_ncgoesimg(iv, filename)
 215             filename = 'goes-15-imager'
 216             call da_read_obs_ncgoesimg(iv, filename)
 217          !end if
 218          !write(unit=stdout,fmt='(a)') 'Finish reading goesimg data'
 219       end if
 220       if (use_gmiobs) then
 221 #if defined(HDF5)
 222          write(unit=stdout,fmt='(a)') 'Reading GMI data in HDF5 format'
 223          call da_read_obs_hdf5gmi (iv, '1B.GPM.GMI','2A.GPM.GMI')
 224 #else
 225          message(1)='To read GMI data, WRFDA must be compiled with HDF5'
 226          call da_error(__FILE__,__LINE__,message(1:1))
 227 #endif
 228       end if
 229    end if
 230
 231    if ( use_filtered_rad ) then
 232       write(unit=stdout,fmt='(a)') 'Reading filtered radiance'
 233       call da_read_filtered_rad (iv)
 234    end if
 235
 236    if ( use_simulated_rad ) then
 237       write(unit=stdout,fmt='(a)') 'Reading simulated radiance'
 238       call da_read_simulated_rad (iv)
 239    end if
 240
 241    if ( use_pseudo_rad ) then
 242       write(unit=stdout,fmt='(a)') 'Reading pseudo radiance from namelist'
 243       call da_read_pseudo_rad (iv)
 244    end if
 245
 246    if (use_kma1dvar) then
 247       do i=1,rtminit_nsensor
 248          filename = ' '
 249          filename='kma1dvar-'//trim(iv%instid(i)%rttovid_string)//'.inv'
 250          write(unit=stdout,fmt='(a,a)')  ' Reading KMA 1dvar innovation from  ', filename
 251          call da_read_kma1dvar (i,iv, ob, filename)
 252       end do
 253    end if
 254
 255    if (thinning) then
 256       do ifgat=1,num_fgat_time
 257            do n=1,iv%num_inst
 258               call destroygrids (n,ifgat)
 259            end do
 260       end do
 261       deallocate(thinning_grid)
 262    end if
 263
 264    ! sorting obs into FGAT time bins
 265
 266    call da_sort_rad(iv)
 267
 268    !-----------------------------------------------------------------------------
 269    ! [3.0] create (smaller) ob structure:
 270    !-----------------------------------------------------------------------------
 271
 272    if (.not. use_kma1dvar) then
 273       do i = 1,  ob % num_inst
 274          ob % instid(i) % num_rad = iv % instid(i) % num_rad
 275          if (ob % instid(i) % num_rad < 1) cycle
 276          allocate (ob % instid(i) % tb(ob % instid(i) % nchan,ob % instid(i)%num_rad))
 277          ob % instid(i) % tb(:,:) = iv % instid(i) % tb_inv(:,:)
 278       end do
 279    end if
 280
 281 ! Calculate DT for Cloudy Radiance DA
 282
 283    if (use_rad .and. crtm_cloud .and. .not. DT_cloud_model) then
 284       its = grid%xp % its
 285       ite = grid%xp % ite
 286       jts = grid%xp % jts
 287       jte = grid%xp % jte
 288       kts = grid%xp % kts
 289       kte = grid%xp % kte
 290
 291       grid%xb%delt(its:ite,jts:jte,kts:kte) = 0.0
 292
 293       do inst= 1, iv % num_inst
 294          do n=1,iv%instid(inst)%num_rad
 295              i = int(iv%instid(inst)%info%i(1,n))
 296              j = int(iv%instid(inst)%info%j(1,n))
 297              grid%xb%delt(i:i+1, j:j+1, kts:kte) = 1800.0
 298          end do
 299       end do
 300    endif
 301
 302    if (trace_use) call da_trace_exit("da_setup_radiance_structures")
 303
 304 end subroutine da_setup_radiance_structures
 305