var/da/da_setup_structures/da_setup_obs_structures.inc

   1 subroutine da_setup_obs_structures( grid, ob, iv, j_cost)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Allocate and read in components of observation structure.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (domain),     intent(inout) :: grid   ! Model data
  10    type (y_type),     intent(out)   :: ob     ! Observation structure.
  11    type (iv_type),    intent(out)   :: iv     ! O-B structure.
  12    type (j_type),     intent(inout) :: j_cost ! Cost function.
  13
  14    integer :: i,j,icode
  15    character (len=14)   :: ob_name
  16
  17    ! thinning variables
  18    integer  :: istart,iend,jstart,jend
  19    real     :: rlonlat(4)
  20
  21    if (trace_use) call da_trace_entry("da_setup_obs_structures")
  22
  23 ! Default values for fmt_info, fmt_srfc, and fmt_each(YRG, 02/25/2009):
  24
  25    fmt_info = '(a12,1x,a19,1x,a40,1x,i6,3(f12.3,11x),6x,a40)'
  26    fmt_srfc = '(7(:,f12.3,i4,f7.2))'
  27    fmt_each = '(3(f12.3,i4,f7.2),11x,3(f12.3,i4,f7.2),11x,1(f12.3,i4,f7.2))'
  28
  29    call da_message((/'Set up observations (ob)'/))
  30
  31    ! Adjust obs switches
  32    if (use_pseudo_rad .or. num_pseudo > 0) then
  33       call da_message((/"Single OBS Test:: Turn off all the OBS switches ***"/))
  34       use_synopobs         = .false.
  35       use_shipsobs         = .false.
  36       use_metarobs         = .false.
  37       use_soundobs         = .false.
  38       use_mtgirsobs        = .false.
  39       use_tamdarobs        = .false.
  40       use_bogusobs         = .false.
  41       use_pilotobs         = .false.
  42       use_airepobs         = .false.
  43       use_geoamvobs        = .false.
  44       use_polaramvobs      = .false.
  45       use_buoyobs          = .false.
  46       use_profilerobs      = .false.
  47       use_satemobs         = .false.
  48       use_gpspwobs         = .false.
  49       use_gpsztdobs        = .false.
  50       use_gpsrefobs        = .false.
  51       use_gpsephobs        = .false.
  52       use_ssmiretrievalobs = .false.
  53       use_ssmitbobs        = .false.
  54       use_ssmt1obs         = .false.
  55       use_ssmt2obs         = .false.
  56       use_qscatobs         = .false.
  57       use_hirs2obs         = .false.
  58       use_hirs3obs         = .false.
  59       use_hirs4obs         = .false.
  60       use_mhsobs           = .false.
  61       use_msuobs           = .false.
  62       use_amsuaobs         = .false.
  63       use_amsubobs         = .false.
  64       use_mwtsobs          = .false.
  65       use_mwhsobs          = .false.
  66       use_atmsobs          = .false.
  67       use_airsobs          = .false.
  68       use_eos_amsuaobs     = .false.
  69       use_hsbobs           = .false.
  70       use_obsgts           = .false.
  71       use_rad              = .false.
  72       use_airsretobs       = .false.
  73       use_rainobs          = .false.
  74       use_iasiobs          = .false.
  75       use_radarobs         = .false.
  76       use_radar_rv         = .false.
  77       use_radar_rf         = .false.
  78 #if (WRF_CHEM == 1)
  79       use_chemic_surfobs     = .false.
  80 #endif
  81       ob_format = ob_format_ascii
  82    end if
  83
  84    if (use_synopobs .OR. use_shipsobs .OR. use_metarobs .OR. use_pilotobs .OR. &
  85       use_profilerobs .OR. use_buoyobs .OR. use_soundobs .OR. use_mtgirsobs .OR. use_bogusobs .OR. &
  86       use_satemobs .OR. use_geoamvobs .OR. use_polaramvobs .OR. use_airepobs .OR. use_tamdarobs .OR. &
  87       use_gpspwobs .OR. use_gpsztdobs .OR. use_gpsrefobs .OR. use_ssmiretrievalobs .OR. &
  88       use_ssmitbobs .OR. use_ssmt1obs .OR. use_ssmt2obs .OR. use_qscatobs .OR. &
  89       use_airsretobs .OR. use_lsac .OR. use_gpsephobs) then
  90
  91       use_obsgts = .true.
  92    else
  93       use_obsgts = .false.
  94    end if
  95
  96    if (use_hirs2obs .OR. use_hirs3obs .OR. use_msuobs .OR. use_amsuaobs .OR. &
  97       use_amsubobs .OR. use_airsobs .OR. use_eos_amsuaobs .OR. &
  98       use_hsbobs .OR. use_kma1dvar .OR. use_filtered_rad .OR. &
  99       use_ssmisobs .OR. use_hirs4obs .OR. use_mhsobs .OR. use_pseudo_rad .OR. &
 100       use_mwtsobs  .OR. use_mwhsobs .OR. use_atmsobs .OR. use_simulated_rad .OR. &
 101       use_iasiobs .OR. use_seviriobs .OR. use_amsr2obs .OR. use_goesimgobs .OR. &
 102       use_ahiobs .OR. use_mwhs2obs .OR. use_gmiobs) then
 103       use_rad = .true.
 104    else
 105       use_rad = .false.
 106    end if
 107
 108    ! test_dm_exact can be set to .true. to force DM_PARALLEL runs
 109    ! to produce results that are bitwise-identical regardless of the number of
 110    ! MPI tasks used.  This is useful for validating that multi-processor runs
 111    ! are not a source of bugs.  Runtime will be much longer.  This option is
 112    ! automatically overridden to .false. for serial or 1-MPI-task runs.
 113
 114    if (test_dm_exact) then
 115       if (num_procs == 1) then
 116          test_dm_exact = .false.
 117          write(unit=stdout,fmt='(A)') &
 118             ' test_dm_exact overridden to .false. for serial or 1-MPI-task run'
 119       end if
 120       ! only implmenented for Sound and Synop, so switch other types off
 121       use_shipsobs         = .false.
 122       use_metarobs         = .false.
 123       use_bogusobs         = .false.
 124       use_pilotobs         = .false.
 125       use_airepobs         = .false.
 126       use_geoamvobs        = .false.
 127       use_polaramvobs      = .false.
 128       use_buoyobs          = .false.
 129       use_profilerobs      = .false.
 130       use_satemobs         = .false.
 131       use_gpspwobs         = .false.
 132       use_gpsztdobs        = .false.
 133       use_gpsrefobs        = .false.
 134       use_gpsephobs        = .false.
 135       use_ssmiretrievalobs = .false.
 136       use_ssmitbobs        = .false.
 137       use_ssmt1obs         = .false.
 138       use_ssmt2obs         = .false.
 139       use_qscatobs         = .false.
 140       use_hirs2obs         = .false.
 141       use_hirs3obs         = .false.
 142       use_hirs4obs         = .false.
 143       use_mhsobs           = .false.
 144       use_msuobs           = .false.
 145       use_amsuaobs         = .false.
 146       use_amsubobs         = .false.
 147       use_mwtsobs          = .false.
 148       use_mwhsobs          = .false.
 149       use_atmsobs          = .false.
 150       use_airsobs          = .false.
 151       use_eos_amsuaobs     = .false.
 152       use_hsbobs           = .false.
 153       use_obsgts           = .false.
 154       use_rad              = .false.
 155       use_airsretobs       = .false.
 156       use_rainobs          = .false.
 157       use_iasiobs          = .false.
 158       use_radarobs         = .false.
 159       use_radar_rv         = .false.
 160       use_radar_rf         = .false.
 161 #if (WRF_CHEM == 1)
 162       use_chemic_surfobs     = .false.
 163 #endif
 164    end if
 165
 166    if ( use_pseudo_rad ) then
 167       use_rad = .true.
 168       thinning         = .false.
 169       qc_rad           = .false.
 170       rtminit_nsensor  = 1
 171       rtminit_platform = pseudo_rad_platid
 172       rtminit_satid    = pseudo_rad_satid
 173       rtminit_sensor   = pseudo_rad_senid
 174    end if
 175
 176    obs_use(1:num_ob_indexes) = .false.
 177    if ( use_airepobs ) obs_use(airep) = .true.
 178    if ( use_bogusobs ) obs_use(bogus) = .true.
 179    if ( use_buoyobs ) obs_use(buoy) = .true.
 180    if ( use_geoamvobs ) obs_use(geoamv) = .true.
 181    if ( use_gpsephobs ) obs_use(gpseph) = .true.
 182    if ( use_gpsrefobs ) obs_use(gpsref) = .true.
 183    if ( use_gpsztdobs .or. use_gpspwobs ) obs_use(gpspw) = .true.
 184    if ( use_metarobs ) obs_use(metar) = .true.
 185    if ( use_mtgirsobs ) obs_use(mtgirs) = .true.
 186    if ( use_pilotobs ) obs_use(pilot) = .true.
 187    if ( use_polaramvobs ) obs_use(polaramv) = .true.
 188    if ( use_profilerobs ) obs_use(profiler) = .true.
 189    if ( use_qscatobs ) obs_use(qscat) = .true.
 190    if ( use_radarobs ) obs_use(radar) = .true.
 191    if ( use_rainobs ) obs_use(rain) = .true.
 192    if ( use_satemobs ) obs_use(satem) = .true.
 193    if ( use_shipsobs ) obs_use(ships) = .true.
 194    if ( use_soundobs ) then
 195       obs_use(sound) = .true.
 196       obs_use(sonde_sfc) = .true.
 197    end if
 198    if ( use_ssmiretrievalobs ) obs_use(ssmi_rv) = .true.
 199    if ( use_ssmitbobs ) obs_use(ssmi_tb) = .true.
 200    if ( use_ssmt1obs ) obs_use(ssmt1) = .true.
 201    if ( use_ssmt2obs ) obs_use(ssmt2) = .true.
 202    if ( use_synopobs ) obs_use(synop) = .true.
 203    if ( use_tamdarobs ) then
 204       obs_use(tamdar) = .true.
 205       obs_use(tamdar_sfc) = .true.
 206    end if
 207
 208    if (sfc_assi_options < 1 .OR. sfc_assi_options > 2) then
 209       write(unit=message(1),fmt='(A,I3)') &
 210          'Invalid sfc_assi_option = ', sfc_assi_options
 211       call da_error(__FILE__,__LINE__,message(1:1))
 212    end if
 213
 214    if ( use_gpsrefobs ) then
 215       if ( ob_format == ob_format_bufr ) then
 216           ! when main conv input is from bufr, gpsro has to be from bufr too
 217           ! overwrite the namelist setting
 218           ob_format_gpsro = ob_format_bufr
 219       else
 220          if ( ob_format == ob_format_ascii .and. &
 221               ob_format_gpsro == ob_format_bufr ) then
 222             call da_message((/'ob_format=2 while ob_format_gpsro=1, will be reading GPSRO data from gpsro.bufr'/))
 223          end if
 224       end if
 225    end if
 226
 227    ! convert pseudo_var to be in lowercase
 228    do i = 1, len_trim(pseudo_var)
 229       icode = ichar(pseudo_var(i:i))
 230       if (icode>=65 .and. icode<=90) then
 231          pseudo_var(i:i)=char(icode + 97 - 65)
 232       end if
 233    end do
 234
 235    if ( num_pseudo > 0 .and.                   &
 236         (trim(adjustl(pseudo_var))=='qcw' .or. &
 237          trim(adjustl(pseudo_var))=='qrn' .or. &
 238          trim(adjustl(pseudo_var))=='qci' .or. &
 239          trim(adjustl(pseudo_var))=='qsn' .or. &
 240          trim(adjustl(pseudo_var))=='qgr') ) then
 241       if ( cloud_cv_options < 2 ) then
 242          write(unit=message(1),fmt='(A,A)') &
 243             'cloud_cv_options should be 2 or 3 for pseudo_var = ', trim(pseudo_var)
 244          call da_error(__FILE__,__LINE__,message(1:1))
 245       end if
 246    end if
 247
 248    pseudo_tpw = .false.
 249    pseudo_ztd = .false.
 250    pseudo_ref = .false.
 251    if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='tpw' ) then
 252       pseudo_tpw    = .true.
 253       use_gpspwobs  = .true.
 254    end if
 255    if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='ztd' ) then
 256       pseudo_ztd    = .true.
 257       use_gpsztdobs = .true.
 258    end if
 259    if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='ref' ) then
 260       pseudo_ref    = .true.
 261       use_gpsrefobs = .true.
 262    end if
 263    ! if pseudo ob for conventional u, v, t, p, q
 264    if ( num_pseudo > 0 ) then
 265       pseudo_uvtpq = .true.
 266       if ( pseudo_tpw .or. pseudo_ztd .or. pseudo_ref ) then
 267          pseudo_uvtpq = .false.
 268       end if
 269    end if
 270    if ( num_pseudo > 0 ) then
 271       !to avoid errors when writing out filtered_obs for undefined variables
 272       anal_type_qcobs = .false.
 273    end if
 274
 275    if ( thin_conv .or. thin_conv_ascii ) then
 276      thin_conv_opt(gpseph) = no_thin
 277      thin_conv_opt(radar) = no_thin
 278      thin_conv_opt(radiance) = no_thin
 279      thin_conv_opt(rain) = no_thin
 280      if ( thin_conv .and. ob_format==ob_format_bufr ) then
 281         ! gpsref horizontal thinning is not implemented for bufr input
 282         thin_conv_opt(gpsref) = no_thin
 283      end if
 284      do i = 1, num_ob_indexes
 285         ! resetting thin_conv_opt to no_thin for ob types that are not used
 286         if ( thin_conv_opt(i) > no_thin .and. &
 287              (.not. obs_use(i)) ) then
 288            thin_conv_opt(i) = no_thin
 289         end if
 290      end do
 291      if ( thin_conv .and. ob_format==ob_format_bufr .and. &
 292           all(thin_conv_opt(:) <= no_thin) ) thin_conv = .false.
 293      if ( thin_conv_ascii .and. ob_format==ob_format_ascii .and. &
 294           all(thin_conv_opt(:) <= no_thin) ) thin_conv_ascii = .false.
 295    end if
 296    if ( ((.not. thin_conv) .and. ob_format==ob_format_bufr) .or. &
 297         ((.not. thin_conv_ascii) .and. ob_format==ob_format_ascii) ) then
 298      thin_conv_opt(:) = no_thin
 299    end if
 300    write(message(1),'(a,30i2)') 'thin_conv_opt = ', thin_conv_opt(:)
 301    call da_message(message(1:1))
 302
 303    ! get lat/lon info (rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid)
 304    ! to be used for making thinning grids for conv (ascii/bufr), radiance and rain ob types.
 305    call init_constants_derived
 306    rlat_min =  r999
 307    rlat_max = -r999
 308    rlon_min =  r999
 309    rlon_max = -r999
 310    istart=MINVAL( grid%i_start(1:grid%num_tiles) )
 311    iend  =MAXVAL( grid%i_end  (1:grid%num_tiles) )
 312    jstart=MINVAL( grid%j_start(1:grid%num_tiles) )
 313    jend  =MAXVAL( grid%j_end  (1:grid%num_tiles) )
 314    do i = istart, iend
 315       do j = jstart, jend
 316          rlat_min=min(rlat_min, grid%xb%lat(i,j))
 317          rlat_max=max(rlat_max, grid%xb%lat(i,j))
 318          if( grid%xb%lon(i,j) < 0.0) then
 319            rlon_min=min(rlon_min, (r360+grid%xb%lon(i,j)))
 320            rlon_max=max(rlon_max, (r360+grid%xb%lon(i,j)))
 321          else
 322            rlon_min=min(rlon_min, grid%xb%lon(i,j))
 323            rlon_max=max(rlon_max, grid%xb%lon(i,j))
 324          endif
 325       enddo
 326    enddo
 327 #ifdef DM_PARALLEL
 328    call mpi_reduce(rlat_min, rlonlat(1), 1, true_mpi_real, mpi_min, root, comm, ierr)
 329    call mpi_reduce(rlon_min, rlonlat(2), 1, true_mpi_real, mpi_min, root, comm, ierr)
 330    call mpi_reduce(rlat_max, rlonlat(3), 1, true_mpi_real, mpi_max, root, comm, ierr)
 331    call mpi_reduce(rlon_max, rlonlat(4), 1, true_mpi_real, mpi_max, root, comm, ierr)
 332    CALL mpi_bcast( rlonlat, 4 , true_mpi_real , root , comm, ierr )
 333    rlat_min = rlonlat(1)
 334    rlon_min = rlonlat(2)
 335    rlat_max = rlonlat(3)
 336    rlon_max = rlonlat(4)
 337 #endif
 338    dlat_grid = rlat_max - rlat_min
 339    dlon_grid = rlon_max - rlon_min
 340
 341    !---------------------------------------------------------------------------
 342    ! [1.0] Setup and read in fields from first guess:
 343    !----------------------------------------------------------------------------
 344
 345    iv%missing = missing
 346    ! iv%ptop    = grid%xb%ptop
 347
 348    iv%total_rad_pixel   = 0
 349    iv%total_rad_channel = 0
 350
 351    iv%info(:)%nlocal = 0
 352    iv%info(:)%ntotal = 0
 353    iv%info(:)%thin_nlocal = 0
 354    iv%info(:)%thin_ntotal = 0
 355    do i=1,num_ob_indexes
 356       iv%info(i)%plocal(:) = 0
 357       iv%info(i)%ptotal(:) = 0
 358       iv%info(i)%thin_plocal(:) = 0
 359       iv%info(i)%thin_ptotal(:) = 0
 360    end do
 361    iv%num_inst  = 0
 362
 363    ob%nlocal(:) = 0
 364    ob%ntotal(:) = 0
 365    ob%num_inst  = 0
 366
 367    iv%info(:)%max_lev = 1
 368
 369    ! get FGAT time slots
 370
 371    allocate ( time_slots(0:num_fgat_time) )
 372    call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
 373                           time_window_max, time_slots, ifgat_ana)
 374
 375    if (use_obsgts) then
 376       ! Conventional obs can be in BUFR or ascii format
 377       if (ob_format == ob_format_bufr) then
 378          call da_message((/'Using PREPBUFR format observation input'/))
 379          call da_setup_obs_structures_bufr (grid, ob, iv)
 380       else if (ob_format == ob_format_ascii) then
 381          call da_message((/'Using ASCII format observation input'/))
 382          call da_setup_obs_structures_ascii (ob, iv, grid)
 383       else
 384          write(unit=message(1),fmt='(A,I3)') 'Invalid ob_format = ', ob_format
 385          call da_error(__FILE__,__LINE__,message(1:1))
 386       end if
 387    end if
 388
 389    if (use_radarobs) then
 390       ! Radar obs are read from separate file(s)
 391          call da_message((/'Using ASCII format radar observation input'/))
 392          call da_setup_obs_structures_radar (grid, ob, iv)
 393    end if
 394
 395    if (use_rainobs .and. var4d) then
 396       call da_message((/'Using ASCII format precipitation observation input'/))
 397       call da_setup_obs_structures_rain (grid, ob, iv)
 398    end if
 399
 400 #if (WRF_CHEM == 1)
 401    if ( (use_chemic_surfobs) ) then
 402       call da_message((/'Using ASCII format chem IC observation input'/))
 403       call da_setup_obs_structures_chem_sfc (ob, iv, grid)
 404    end if
 405 #endif
 406
 407    if (use_rad) then
 408 #if defined(CRTM) || defined(RTTOV)
 409       ! Radiance files can only be in BUFR
 410       call da_message((/'Using NCEP BUFR radiance 1b input'/))
 411       if ( use_amsr2obs ) then
 412          call da_message((/'Using AMSR2 radiance input in HDF5 format'/))
 413       end if
 414       if ( use_gmiobs ) then
 415          call da_message((/'Using GMI radiance input in HDF5 format'/))
 416       end if
 417       call da_setup_radiance_structures(grid, ob, iv)
 418 #else
 419       call da_error(__FILE__,__LINE__,(/"options to use radiance data are turned on in wrfvar4 namelist, the code needs to be compiled with CRTM or RTTOV library"/))
 420 #endif
 421    end if
 422
 423    if ( num_pseudo > 0 ) then
 424       call da_setup_pseudo_obs(grid, iv, ob)
 425    end if
 426
 427    ! Summarize observations
 428
 429    write(unit=stdout, fmt='(a)')  'Observation summary'
 430    do i=1,num_fgat_time
 431       write(unit=stdout, fmt='(3x,2(a,i3))') 'ob time', i,' num_ob_indexes=',num_ob_indexes
 432       do j=1,num_ob_indexes
 433          if(j.eq.27) cycle        ! skip tamdar_sfc (HA)
 434          if (use_gpsztdobs .and.obs_names(j) == 'gpspw         ' ) then
 435             ob_name = 'gpsztd        '
 436          else
 437             ob_name = obs_names(j)
 438          endif
 439          if(iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1) > 0) then
 440          if ( thin_conv_ascii .and. ob_format == ob_format_ascii ) then
 441             if ( j == radiance ) cycle !temporary fix to not print incorrect counts for radiance
 442             write(unit=stdout, fmt='(2x,2(a,i10),a,i8,a)') &
 443                ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global, ', &
 444                iv%info(j)%thin_ptotal(i) - iv%info(j)%thin_ptotal(i-1), ' global (post-thinning), ', &
 445                iv%info(j)%thin_plocal(i) - iv%info(j)%thin_plocal(i-1), ' local (post-thinning)'
 446          else
 447 #if (WRF_CHEM == 1)
 448             write(unit=stdout, fmt='(i3,3x,a,i10,a,i8,a,i8,a,i8,a)') &
 449                j,ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', &
 450                iv%info(j)%plocal(i) - iv%info(j)%plocal(i-1), ' local', &
 451                iv%info(chemic_surf)%ntotal, ' chem total', iv%info(chemic_surf)%nlocal, ' chem local'
 452 #else
 453             write(unit=stdout, fmt='(6x,a,i10,a,i8,a)') &
 454                ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', &
 455                iv%info(j)%plocal(i) - iv%info(j)%plocal(i-1), ' local'
 456 #endif
 457          end if ! thin_conv_ascii
 458          end if
 459       end do
 460    end do
 461    write(unit=stdout, fmt='(a)') ' '
 462
 463    ! Get horizontal interpolation weights.
 464
 465    call da_setup_obs_interp_wts (iv)
 466
 467    if (trace_use) call da_trace_exit("da_setup_obs_structures")
 468
 469 end subroutine da_setup_obs_structures
 470
 471