var/da/da_obs_io/da_read_obs_bufr_satwnd.inc

   1 subroutine da_read_obs_bufr_satwnd(filename, iv, ob)
   2
   3 implicit none
   4
   5 character(len=*), intent(in)    :: filename
   6 type (iv_type),   intent(inout) :: iv
   7 type (y_type),    intent(inout) :: ob
   8
   9 integer, parameter :: i_kind   = selected_int_kind(8)
  10 integer, parameter :: r_double = selected_real_kind(15) ! double precision
  11 integer, parameter :: r_kind   = selected_real_kind(15) ! double precision
  12
  13 integer(i_kind), parameter :: missing_i = -999
  14 real(r_kind),    parameter :: missing_r = -999.0
  15
  16 integer(i_kind), parameter :: nstring   = 50
  17 integer(i_kind), parameter :: nmeta     = 3  ! number of metadata
  18 integer(i_kind), parameter :: nvars     = 2  ! number of obs variables
  19
  20 character(len=nstring), dimension(nmeta) :: name_meta = &
  21    (/                        &
  22       'air_pressure       ', &
  23       'latitude           ', &
  24       'longitude          '  &
  25    /)
  26
  27 character(len=nstring), dimension(nvars) :: name_vars = &
  28    (/                        &
  29       'eastward_wind      ', &
  30       'northward_wind     '  &
  31    /)
  32
  33 integer(i_kind) :: idx_p = 1
  34 integer(i_kind) :: idx_lat = 2
  35 integer(i_kind) :: idx_lon = 3
  36 integer(i_kind) :: idx_u = 1
  37 integer(i_kind) :: idx_v = 2
  38
  39 type datalink_satwnd
  40    character(len=nstring)    :: stid       ! station identifier
  41    character(len=nstring)    :: name       ! instrument name
  42    character(len=nstring)    :: date_char  ! ccyy-mm-dd_hh:mm:ss
  43    integer(i_kind)           :: satid      ! satellite identifier
  44    integer(i_kind)           :: rptype     ! prepbufr report type
  45    integer(i_kind)           :: ifgat      ! index of time slot
  46    integer(i_kind)           :: itg        ! index of thinning grid
  47    integer(i_kind)           :: qm         ! quality marker
  48    real(r_double)            :: tdiff      ! time difference between obs_time and analysis_time
  49    real(r_kind)              :: err        ! ob error
  50    real(r_kind)              :: lat        ! latitude in degree
  51    real(r_kind)              :: lon        ! longitude in degree
  52    real(r_kind)              :: prs        ! pressure
  53    real(r_kind)              :: landsea    ! land sea qualifier 0:land, 1:sea, 2:coast
  54    real(r_kind)              :: satzen     ! satellite zenith angle in degree
  55    real(r_kind)              :: wspd
  56    real(r_kind)              :: wdir
  57    real(r_kind)              :: uwind
  58    real(r_kind)              :: vwind
  59    real(r_kind)              :: cvwd       ! coefficient of variation
  60    real(r_kind)              :: pccf1      ! percent confidence, Quality Index without forecast (qifn)
  61    real(r_kind)              :: pccf2      ! percent confidence, Estimated Error (EE) in m/s converted to percent confidence
  62
  63    type (datalink_satwnd), pointer :: next ! pointer to next data
  64 end type datalink_satwnd
  65
  66 type(datalink_satwnd), pointer :: rhead=>null(), rlink=>null()
  67
  68 real(r_kind), parameter  :: r8bfms = 9.0E08  ! threshold to check for BUFR missing value
  69 real(r_kind), parameter  :: pi = acos(-1.0)
  70
  71 integer(i_kind), parameter :: ntime = 6      ! number of data to read in timestr
  72 integer(i_kind), parameter :: ninfo = 5      ! number of data to read in infostr
  73 integer(i_kind), parameter :: nlalo = 2      ! number of data to read in lalostr
  74 integer(i_kind), parameter :: ndata = 3      ! number of data to read in datastr
  75 integer(i_kind), parameter :: nqc1  = 2      ! number of data to read in qc1str
  76 integer(i_kind), parameter :: nqc2  = 2      ! number of data to read in qc2str
  77
  78 character(len=80) :: timestr, infostr, lalostr, datastr, qc1str, qc2str
  79
  80 real(r_double), dimension(ntime)  :: timedat
  81 real(r_double), dimension(ninfo)  :: infodat
  82 real(r_double), dimension(nlalo)  :: lalodat
  83 real(r_double), dimension(ndata)  :: wdata
  84 real(r_double), dimension(nqc1,2) :: qc1dat
  85 real(r_double), dimension(nqc2,4) :: qc2dat
  86
  87 integer(i_kind), parameter :: nmsgtyp = 10  ! number of message types to process
  88 ! message types that are not processed into prepbufr
  89 character(len=8), dimension(nmsgtyp) :: message_types = &
  90    (/ 'NC005030',  &
  91       'NC005031',  &
  92       'NC005032',  &
  93       'NC005034',  &
  94       'NC005039',  &
  95       'NC005080',  &
  96       'NC005081',  &  ! not tested
  97       'NC005090',  &
  98       'NC005091',  &
  99       'NC005072'   &
 100    /)
 101 character(len=8)  :: subset
 102
 103 integer(i_kind) :: iunit, iost, iret
 104 integer(i_kind) :: k, kx, ilev, itype, ivar
 105 integer(i_kind) :: idate
 106 integer(i_kind) :: num_report_infile
 107 integer(i_kind) :: ireadmg, ireadsb
 108
 109 integer(i_kind) :: iyear, imonth, iday, ihour, imin, isec
 110 real(r_double)  :: obs_time, analysis_time
 111
 112 logical      :: outside, outside_all
 113 logical      :: use_errtable
 114 real(r_kind) :: coef
 115 real(r_kind) :: pob
 116 real(r_kind) :: xlat, xlon, prs
 117 real(r_kind) :: angearth
 118 type(info_type) :: xinfo
 119 type(model_loc_type) :: xloc1
 120
 121 integer,                  allocatable :: itmp(:)
 122 character(len=nstring),   allocatable :: str_stid(:)
 123 character(len=nstring),   allocatable :: str_name(:)
 124 character(len=nstring),   allocatable :: date_char(:)     ! date_time string in wrf format
 125 type(model_loc_type),     allocatable :: xloc(:)
 126 logical,                  allocatable :: inside_domain(:)
 127 logical,                  allocatable :: inside_patch(:)
 128 integer(i_kind),          allocatable :: iwindow(:)
 129 integer(i_kind),          allocatable :: rptype(:)
 130 real(r_kind),    allocatable :: tdiff_avg(:)
 131 real(r_kind),    allocatable :: num_val(:)
 132 real(r_kind),    allocatable :: num_err(:)
 133 real(r_kind),    allocatable :: num_meta(:)
 134 real(r_kind),    allocatable :: sdata_val(:,:)
 135 real(r_kind),    allocatable :: sdata_err(:,:)
 136 real(r_kind),    allocatable :: smeta(:,:)
 137 character(len=14) :: cdate14       ! ccyymmddhhnnss used in da_advance_time
 138 character(len=10) :: adate10       ! ccyymmddhh     used in da_advance_time
 139 character(len=14) :: tdiff_char
 140
 141 integer(i_kind) :: i, ii, kk, n, isize
 142 integer(i_kind) :: num_outside_domain, num_outside_time
 143 integer(i_kind) :: num_qc_fail
 144 integer(i_kind) :: ntotal, nlocal, max_lev, nlev, nlevel
 145 integer(i_kind) :: itotal, ilocal, ifgat
 146 integer(i_kind) :: ntotal_ifgat(0:num_fgat_time)
 147 logical         :: fexist
 148
 149 real(r_kind) :: err_val
 150 real(r_kind) :: uob, vob
 151 real(r_kind) :: u_grid, v_grid
 152
 153 ! for thinning
 154 real(r_kind)    :: dlat_earth, dlon_earth, crit
 155 integer(i_kind) :: nobs, icnt, num_date
 156 integer(i_kind) :: iobs, itt, itx, iout
 157 integer(i_kind) :: num_thinned, num_keep
 158 logical :: iuse
 159 real(r_kind)    :: xval
 160 real(r_kind)    :: c1, c2
 161 ! for superobbing in vertical
 162 real(r_kind)    :: zk
 163 integer(i_kind) :: nplev
 164 real(r_kind), allocatable :: plev(:)
 165
 166 continue ! end of declaration
 167
 168 inquire(file=trim(filename), exist=fexist)
 169 if ( .not. fexist ) then
 170    write(unit=message(1),fmt='(a)') trim(filename)//' does not exist.'
 171    call da_warning(__FILE__,__LINE__,message(1:1))
 172    return
 173 end if
 174
 175 if ( thin_conv_opt(polaramv) == thin_single ) then
 176    thin_conv_opt(polaramv) = thin_multi
 177    write(unit=message(1),fmt='(a,2i5)') 'thin_conv_opt(polaramv)=1 is not implemented, resetting to thin_conv_opt(polaramv)=2'
 178    call da_message(message(1:1))
 179 end if
 180
 181 timestr = 'YEAR MNTH DAYS HOUR MINU SECO'
 182 infostr = 'SAID LSQL SAZA OGCE SWCM'
 183 lalostr = 'CLATH CLONH'
 184 datastr = 'PRLC WDIR WSPD'
 185 qc1str  = 'TCOV CVWD'
 186 qc2str  = 'GNAPS PCCF'
 187
 188 ! oetab(300,33,6) processed in da_setup_obs_structures_bufr.inc
 189 use_errtable = allocated(oetab)
 190
 191 num_report_infile  = 0
 192
 193 allocate (itmp(1)) ! for broadcasting nobs
 194
 195 ! only rootproc reads from file
 196 ! thinning check is also done on rootproc only
 197 if ( rootproc ) then
 198
 199    write(unit=message(1),fmt='(a)') 'Reading from '//trim(filename)
 200    call da_message(message(1:1))
 201
 202    ! open bufr file
 203    call da_get_unit(iunit)
 204    open (unit=iunit, file=trim(filename), &
 205          iostat=iost, form='unformatted', status='old')
 206    if (iost /= 0) then
 207       write(unit=message(1),fmt='(a,i5,a)') &
 208          "Error",iost," opening BUFR obs file "//trim(filename)
 209       call da_warning(__FILE__,__LINE__,message(1:1))
 210       return
 211    end if
 212
 213    call openbf(iunit,'IN',iunit)
 214    call datelen(10)
 215    call readmg(iunit,subset,idate,iret)
 216
 217    if ( iret /= 0 ) then
 218       write(unit=message(1),fmt='(A,I5,A)') &
 219          "Error",iret," reading BUFR obs file "//trim(filename)
 220       call closbf(iunit)
 221       close(iunit)
 222       call da_free_unit(iunit)
 223       call da_warning(__FILE__,__LINE__,message(1:1))
 224       return
 225    end if
 226    rewind(iunit)
 227
 228    write(unit=message(1),fmt='(a,i10)') 'da_read_obs_bufr_satwnd: '//trim(filename)//' file date is: ', idate
 229    call da_message(message(1:1))
 230
 231    read (analysis_date,'(i4,4(1x,i2))') iyear, imonth, iday, ihour, imin
 232    call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time)
 233    adate10 = analysis_date(1:4)//analysis_date(6:7)//analysis_date(9:10)//analysis_date(12:13)
 234
 235    ! initialize counters
 236    num_outside_domain = 0
 237    num_outside_time = 0
 238    num_qc_fail = 0
 239    num_thinned = 0
 240    num_keep = 0
 241
 242    if ( thin_conv_opt(polaramv) > no_thin ) then
 243       call cleangrids_conv(polaramv)
 244       iobs = 0 ! initialize an intent(inout) variable used by map2grids_conv
 245       ! itxmax is calculated in subroutine make3grids called in da_setup_obs_structures_bufr
 246       ! itxmax is horizontal thinning grids when thin_conv_opt/=thin_superob_hv
 247       ! itxmax is horizontal*vertical thinning grids when thin_conv_opt==thin_superob_hv
 248       write(unit=message(1),fmt='(a,i8)') &
 249         'da_read_obs_bufr_satwnd: number of thinning grids = ', thinning_grid_conv(polaramv)%itxmax
 250       call da_message(message(1:1))
 251       if ( thin_conv_opt(polaramv) == thin_superob_hv ) then
 252          !to-do: nplev and plev should be set in make3grids and stored in thinning_grid_conv structure
 253          nplev = int(1200/thin_mesh_vert_conv(polaramv))
 254          allocate (plev(nplev))
 255          do k = 1, nplev
 256             plev(k) = 1200.0 - (k-1)*thin_mesh_vert_conv(polaramv)
 257          end do
 258       end if
 259    end if ! thinning
 260
 261    if ( .not. associated(rhead) ) then
 262       nullify ( rhead )
 263       allocate ( rhead )
 264       nullify ( rhead%next )
 265    end if
 266
 267    if ( .not. associated(rlink) ) then
 268       rlink => rhead
 269    else
 270       allocate ( rlink%next )
 271       rlink => rlink%next
 272       nullify ( rlink%next )
 273    end if
 274
 275    msg_loop: do while (ireadmg(iunit,subset,idate)==0)
 276 !print*,subset
 277       if ( ufo_vars_getindex(message_types, subset) <= 0 ) then
 278          ! skip types other than GOES-16/17, AVHRR (METOP/NOAA), VIIRS (NPP/NOAA-20) AMVs
 279          ! that are included in prepbufr
 280          cycle msg_loop
 281       end if
 282       read_loop_subset: do while (ireadsb(iunit)==0)
 283
 284          num_report_infile = num_report_infile + 1
 285
 286          call fill_datalink(rlink, missing_r, missing_i)
 287
 288          call ufbint(iunit,timedat,ntime,1,iret,timestr)
 289
 290          iyear  = nint(timedat(1))
 291          imonth = nint(timedat(2))
 292          iday   = nint(timedat(3))
 293          ihour  = nint(timedat(4))
 294          imin   = nint(timedat(5))
 295          isec   = min(59, nint(timedat(6))) ! raw BUFR data that has SECO = 60.0 SECOND
 296                                             ! that was probably rounded from 59.x seconds
 297                                             ! reset isec to 59 rather than advancing one minute
 298          write(unit=rlink%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 299             iyear, '-', imonth, '-', iday, '_', ihour, ':', imin, ':', isec
 300          call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 301          if ( obs_time < time_slots(0) .or.  &
 302               obs_time >= time_slots(num_fgat_time) ) then
 303             num_outside_time = num_outside_time + 1
 304             cycle read_loop_subset
 305          end if
 306          ! determine fgat index
 307          do ifgat = 1, num_fgat_time
 308             if ( obs_time >= time_slots(ifgat-1) .and.  &
 309                  obs_time  < time_slots(ifgat) ) then
 310                 exit
 311             end if
 312          end do
 313          rlink%ifgat = ifgat
 314          if ( num_fgat_time > 1 ) then
 315             rlink%tdiff = obs_time + float(isec)/60.0 - &
 316                (time_slots(0) + (ifgat-1)*(time_slots(num_fgat_time)-time_slots(0))/float(num_fgat_time-1))
 317          else
 318             rlink%tdiff = obs_time + float(isec)/60.0 - analysis_time
 319          end if
 320          rlink%tdiff = rlink%tdiff/60.0 ! minute to hour
 321
 322          if ( subset == 'NC005080' .or. &
 323               subset == 'NC005090' ) then
 324             call ufbint(iunit,lalodat,nlalo,1,iret,'CLAT CLON')
 325          else
 326             call ufbint(iunit,lalodat,nlalo,1,iret,lalostr)  ! CLATH CLONH
 327          end if
 328
 329          xlat = lalodat(1)
 330          xlon = lalodat(2)
 331          if ( abs(xlat) > 90.0 .or. abs(xlon) > 360.0 ) cycle read_loop_subset
 332          if ( xlon > 180.0 ) xlon = xlon - 360.0
 333          xinfo%lat = xlat
 334          xinfo%lon = xlon
 335          xinfo%lat = max(xinfo%lat, -89.995)
 336          xinfo%lat = min(xinfo%lat,  89.995)
 337          call da_llxy(xinfo, xloc1, outside, outside_all)
 338          if ( outside_all ) then
 339             num_outside_domain = num_outside_domain + 1
 340             cycle read_loop_subset
 341          end if
 342          rlink % lat = xlat
 343          rlink % lon = xlon
 344
 345          call ufbint(iunit,infodat,ninfo,1,iret,infostr)  ! SAID LSQL SAZA OGCE SWCM
 346          call ufbint(iunit,wdata,ndata,1,iret,datastr)    ! PRLC WDIR WSPD
 347          call ufbrep(iunit,qc1dat,nqc1,2,iret,qc1str)     ! 2 replications TCOV CVWD
 348          call ufbrep(iunit,qc2dat,nqc2,4,iret,qc2str)     ! 4 replications GNAPS PCCF
 349
 350          if ( wdata(1) > r8bfms ) cycle read_loop_subset  ! missing pressure
 351
 352          rlink % satid  = nint(infodat(1))  ! SAID satellite identifier
 353          ! construct rptype and stid from subset and satid
 354          call set_rptype_satwnd(subset, rlink%satid, rlink%rptype, rlink%stid, rlink%name)
 355
 356          if ( wdata(1) < r8bfms ) rlink % prs  = wdata(1)   ! PRLC pressure in Pa
 357          if ( wdata(2) < r8bfms ) rlink % wdir = wdata(2)
 358          if ( wdata(3) < r8bfms ) rlink % wspd = wdata(3)
 359
 360          if ( infodat(4) < r8bfms ) rlink % landsea = infodat(4)  ! LSQL land sea qualifier 0:land, 1:sea, 2:coast
 361          if ( infodat(5) < r8bfms ) rlink % satzen  = infodat(5)  ! SAZA satellite zenith angle (degree)
 362
 363          if ( qc1dat(2,1) < r8bfms ) rlink % cvwd  = qc1dat(2,1)
 364          if ( qc2dat(2,2) < r8bfms ) rlink % pccf1 = qc2dat(2,2)
 365          if ( qc2dat(2,4) < r8bfms ) rlink % pccf2 = qc2dat(2,4)
 366
 367          call filter_obs_satwnd(rlink%qm, rlink%rptype, rlink%prs, &  ! set rlink%qm
 368             rlink%satzen, rlink%landsea, rlink%pccf1, rlink%pccf2, rlink%cvwd, rlink%wspd)
 369
 370          if ( rlink%qm > qmarker_retain ) then
 371             num_qc_fail = num_qc_fail + 1
 372             cycle read_loop_subset
 373          end if
 374
 375          if ( thin_conv_opt(polaramv) > no_thin ) then
 376             !crit = abs(tdiff(ii))
 377             crit = 1.0
 378             dlat_earth = rlink%lat
 379             dlon_earth = rlink%lon
 380             if ( dlon_earth < 0.0 ) dlon_earth = dlon_earth + 360.0
 381             if ( dlon_earth >= 360.0 ) dlon_earth = dlon_earth - 360.0
 382             dlat_earth = dlat_earth * deg2rad
 383             dlon_earth = dlon_earth * deg2rad
 384             if ( thin_conv_opt(polaramv) /= thin_superob_hv ) then
 385                call map2grids_conv(polaramv,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 386             else
 387                zk = get_zk(nplev, plev, rlink%prs*0.01) ! 0.01 for Pa to hPa
 388                call map2grids_conv(polaramv,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse,zk=zk,thin_3d=.true.,nlev=nplev)
 389             end if
 390             if ( thin_conv_opt(polaramv) /= thin_superob .and. &
 391                  thin_conv_opt(polaramv) /= thin_superob_hv ) then
 392                if ( .not. iuse ) then
 393                   num_thinned = num_thinned + 1
 394                   if ( print_detail_obs ) then
 395                      write(unit=message(1),fmt='(a,2(1x,f8.3),1x,f8.3,a)')  &
 396                         rlink%stid, rlink%lat, rlink%lon, rlink%prs*0.01, '  -> thinned'
 397                   end if
 398                   cycle read_loop_subset
 399                end if
 400             end if ! thin_conv_opt /= thin_superob, thin_superob_hv
 401             rlink % itg = itx
 402          end if ! thin_conv_opt(polaramv) > nothin
 403
 404          num_keep = num_keep + 1
 405
 406          if ( rlink%wdir >= 0.0 .and. rlink%wdir <= 360.0 .and. &
 407               rlink%wspd < r8bfms ) then
 408             angearth = rlink%wdir * pi / 180.0
 409             rlink%uwind = -1.0 * rlink%wspd * sin(angearth)
 410             rlink%vwind = -1.0 * rlink%wspd * cos(angearth)
 411          end if
 412
 413          if ( use_errtable ) then
 414             kx = rlink % rptype
 415             pob = rlink % prs * 0.01 ! Pa to hPa
 416             do ilev = 1, 32
 417                if ( pob >= oetab(kx,ilev+1,1) .and. pob <= oetab(kx,ilev,1) ) then
 418                   coef = (pob-oetab(kx,ilev,1))/(oetab(kx,ilev,1)-oetab(kx,ilev+1,1))
 419                   rlink % err = (1.0+coef)*oetab(kx,ilev,4)-coef*oetab(kx,ilev+1,4) !uv
 420                   exit
 421                end if
 422             end do
 423             ! in case errors are not set for some report types, use uv_error_val from namelist
 424             if ( rlink % err > 100.0 ) rlink % err = uv_error_val(polaramv)  ! m/s
 425          else
 426             rlink % err = uv_error_val(polaramv)  ! m/s
 427          end if
 428
 429          allocate ( rlink%next )
 430          rlink => rlink%next
 431          nullify ( rlink%next )
 432
 433       end do read_loop_subset ! ireadsb
 434    end do msg_loop ! ireadmg
 435
 436    if ( thin_conv_opt(polaramv) == thin_superob_hv ) then
 437       deallocate (plev)
 438    end if
 439
 440    call closbf(iunit)
 441    close(iunit)
 442    call da_free_unit(iunit)
 443
 444    write(message(1),'(a,a,a,i10)') 'da_read_obs_bufr_satwnd: num_report_infile ', trim(filename), ' : ', num_report_infile
 445    call da_message(message(1:1))
 446    if ( thin_conv_opt(polaramv) == thin_superob .or. &
 447         thin_conv_opt(polaramv) == thin_superob_hv ) then
 448       write(unit=message(1),fmt='(a, 5i8)') &
 449         'da_read_obs_bufr_satwnd: num_outside_domain, num_outside_time, num_qc_fail, num_superob = ', &
 450          num_outside_domain, num_outside_time, num_qc_fail, iobs
 451       call da_message(message(1:1))
 452    else
 453       write(unit=message(1),fmt='(a, 5i8)') &
 454         'da_read_obs_bufr_satwnd: num_outside_domain, num_outside_time, num_qc_fail, num_thinned = ', &
 455          num_outside_domain, num_outside_time, num_qc_fail, num_thinned
 456       call da_message(message(1:1))
 457    end if
 458    if ( thin_conv_opt(polaramv) == no_thin .or. &
 459         thin_conv_opt(polaramv) == thin_multi ) then
 460       nobs = num_keep
 461    else
 462       nobs = iobs
 463    end if
 464
 465    itmp(1) = nobs
 466 end if ! rootproc
 467
 468 call wrf_dm_bcast_integer(itmp, 1)
 469 nobs = itmp(1)
 470 deallocate (itmp)
 471
 472 !print*, 'nobs = ', nobs
 473 if ( nobs <= 0 ) then
 474    write(unit=message(1),fmt='(a)') ' No valid obs found'
 475    call da_warning(__FILE__,__LINE__,message(1:1))
 476    return
 477 end if
 478
 479 ! superobbing is done on rootproc
 480 if ( rootproc) then
 481
 482    allocate (sdata_val(nobs, nvars))
 483    allocate (sdata_err(nobs, nvars))
 484    allocate (smeta(nobs, nmeta))
 485    allocate (iwindow(nobs))
 486    allocate (date_char(nobs))
 487    allocate (str_stid(nobs))
 488    allocate (str_name(nobs))
 489
 490    if ( thin_conv_opt(polaramv) == no_thin .or. &
 491         thin_conv_opt(polaramv) == thin_multi ) then
 492       allocate (rptype(nobs))
 493       icnt = 0
 494       rlink => rhead
 495       do while ( associated(rlink) )
 496          icnt = icnt + 1
 497          sdata_val(icnt,idx_u) = rlink%uwind
 498          sdata_val(icnt,idx_v) = rlink%vwind
 499          sdata_err(icnt,idx_u) = rlink%err
 500          sdata_err(icnt,idx_v) = rlink%err
 501          smeta(icnt,idx_p) = rlink%prs
 502          smeta(icnt,idx_lat) = rlink%lat
 503          smeta(icnt,idx_lon) = rlink%lon
 504          iwindow(icnt) = rlink%ifgat
 505          rptype(icnt) = rlink%rptype
 506          date_char(icnt) = rlink%date_char
 507          str_stid(icnt) = rlink%stid
 508          str_name(icnt) = rlink%name
 509          rlink => rlink%next
 510       end do
 511    else if ( thin_conv_opt(polaramv) == thin_superob .or. &
 512              thin_conv_opt(polaramv) == thin_superob_hv ) then
 513       sdata_val = missing_r  ! initialize
 514       sdata_err = missing_r  ! initialize
 515       smeta = missing_r  ! initialize
 516       allocate (tdiff_avg(nobs))
 517       tdiff_avg = 0.0
 518       iwindow = -1
 519       iobs = 0 ! re-use this counter
 520       allocate (num_val(nvars))
 521       allocate (num_err(nvars))
 522       allocate (num_meta(nmeta))
 523       ! loop over thinning grids
 524       itx_loop: do i = 1, thinning_grid_conv(polaramv)%itxmax
 525          if ( thinning_grid_conv(polaramv)%icount(i) <= 0 ) cycle itx_loop
 526          ! increment the counter for the thinning grid that has obs
 527          iobs = iobs + 1
 528          ! reset counter for obs within each thinning grid
 529          icnt = 0
 530          num_date = 0
 531          num_val(:) = 0
 532          num_err(:) = 0
 533          num_meta(:) = 0
 534          ! loop over nobs to process the obs within the thinning grid
 535          rlink => rhead
 536          nobs_loop_1: do while ( associated(rlink) )
 537             if ( rlink%itg /= i ) then
 538                rlink => rlink%next
 539                cycle nobs_loop_1
 540             end if
 541             icnt = icnt + 1
 542             if ( icnt == 1 ) then ! first ob in the thinning grid
 543               ifgat = rlink%ifgat
 544               str_stid(iobs) = rlink%stid
 545               str_name(iobs) = rlink%name
 546             end if
 547             ! only average obs in the same time window
 548             if ( rlink%ifgat /= ifgat ) then
 549                rlink => rlink%next
 550                cycle nobs_loop_1
 551             end if
 552             ! average dates
 553             num_date = num_date + 1
 554             c1 = 1.0/num_date
 555             c2 = (num_date-1)*c1
 556             tdiff_avg(iobs) = c2*tdiff_avg(iobs) + rlink%tdiff*c1
 557             do ivar = 1, nvars
 558                if ( trim(name_vars(ivar)) == 'eastward_wind' ) then
 559                   xval = rlink%uwind
 560                else if ( trim(name_vars(ivar)) == 'northward_wind' ) then
 561                   xval = rlink%vwind
 562                end if
 563                if ( abs(xval-missing_r) > 0.1 ) then
 564                   num_val(ivar) = num_val(ivar) + 1
 565                   c1 = 1.0/num_val(ivar)
 566                   c2 = (num_val(ivar)-1)*c1
 567                   sdata_val(iobs, ivar) = c2*sdata_val(iobs, ivar) + xval*c1
 568                end if
 569                xval = rlink%err
 570                if ( abs(xval-missing_r) > 0.1 ) then
 571                   num_err(ivar) = num_err(ivar) + 1
 572                   c1 = 1.0/num_err(ivar)
 573                   c2 = (num_err(ivar)-1)*c1
 574                   sdata_err(iobs, ivar) = c2*sdata_err(iobs, ivar) + xval*c1
 575                end if
 576             end do
 577             do ivar = 1, nmeta
 578                if ( trim(name_meta(ivar)) == 'air_pressure' ) then
 579                   xval = rlink%prs
 580                else if ( trim(name_meta(ivar)) == 'latitude' ) then
 581                   xval = rlink%lat
 582                else if ( trim(name_meta(ivar)) == 'longitude' ) then
 583                   xval = rlink%lon
 584                end if
 585                if ( abs(xval-missing_r) > 0.1 ) then
 586                   num_meta(ivar) = num_meta(ivar) + 1
 587                   c1 = 1.0/num_meta(ivar)
 588                   c2 = (num_meta(ivar)-1)*c1
 589                   smeta(iobs, ivar) = c2*smeta(iobs, ivar) + xval*c1
 590                end if
 591             end do
 592             if ( icnt == thinning_grid_conv(polaramv)%icount(i) ) then
 593                rlink => null()
 594                exit nobs_loop_1
 595             end if
 596             rlink => rlink%next
 597          end do nobs_loop_1
 598          iwindow(iobs) = ifgat
 599          tdiff_char = ''
 600          write(tdiff_char,'(i4,a1)') int(tdiff_avg(iobs)*60.0), 'm'
 601          call da_advance_time(adate10, trim(adjustl(tdiff_char)), cdate14)
 602          date_char(iobs) = cdate14(1:4)//'-'//cdate14(5:6)//'-'//cdate14(7:8)// &
 603                            '_'//cdate14(9:10)//':'//cdate14(11:12)//':'//cdate14(13:14)
 604       end do itx_loop
 605       deallocate (num_meta)
 606       deallocate (num_err)
 607       deallocate (num_val)
 608       deallocate (tdiff_avg)
 609    end if ! thin_conv_opt(polaramv)
 610
 611    ! done with rlink
 612    ! release the linked list
 613    rlink => rhead
 614    do while ( associated(rlink) )
 615       rhead => rlink%next
 616       if ( associated (rlink) ) deallocate (rlink)
 617       rlink => rhead
 618    end do
 619    nullify (rhead)
 620
 621 end if ! rootproc
 622
 623 if ( .not. allocated(sdata_val) ) then
 624    allocate (sdata_val(nobs, nvars))
 625 end if
 626 isize = nobs*nvars
 627 call wrf_dm_bcast_real(sdata_val, isize)
 628 if ( .not. allocated(sdata_err) ) then
 629    allocate (sdata_err(nobs, nvars))
 630 end if
 631 call wrf_dm_bcast_real(sdata_err, isize)
 632 if ( .not. allocated(smeta) ) then
 633    allocate (smeta(nobs, nmeta))
 634 end if
 635 isize = nobs*nmeta
 636 call wrf_dm_bcast_real(smeta, isize)
 637 isize = nobs
 638 if ( .not. allocated(iwindow) ) then
 639    allocate (iwindow(nobs))
 640 end if
 641 call wrf_dm_bcast_integer(iwindow, isize)
 642 if ( thin_conv_opt(polaramv) == no_thin .or. &
 643      thin_conv_opt(polaramv) == thin_multi ) then
 644    if ( .not. allocated(rptype) ) then
 645       allocate (rptype(nobs))
 646    end if
 647    call wrf_dm_bcast_integer(rptype, isize)
 648 end if
 649 if ( .not. allocated(str_stid) ) then
 650    allocate (str_stid(nobs))
 651 end if
 652 if ( .not. allocated(str_name) ) then
 653    allocate (str_name(nobs))
 654 end if
 655 if ( .not. allocated(date_char) ) then
 656    allocate (date_char(nobs))
 657 end if
 658 do i = 1, nobs
 659    call wrf_dm_bcast_string(str_stid(i), nstring)
 660    call wrf_dm_bcast_string(str_name(i), nstring)
 661    call wrf_dm_bcast_string(date_char(i), nstring)
 662 end do
 663
 664 allocate (xloc(nobs))
 665 allocate (inside_domain(nobs))
 666 allocate (inside_patch(nobs))
 667 inside_domain(:) = .false.
 668 inside_patch(:)  = .false.
 669
 670 nlocal = 0
 671 ntotal = 0
 672 ntotal_ifgat(0:num_fgat_time) = 0
 673
 674 check_loop_2: do ii = 1, nobs
 675
 676    ! check loc
 677    if ( abs(smeta(ii,idx_lat)) > 90.0 .or. abs(smeta(ii,idx_lon)) > 360.0 ) then
 678       cycle check_loop_2
 679    end if
 680    if ( smeta(ii,idx_lon) > 180.0 ) smeta(ii,idx_lon) = smeta(ii,idx_lon) - 360.0
 681    xinfo%lat = smeta(ii,idx_lat)
 682    xinfo%lon = smeta(ii,idx_lon)
 683    xinfo%lat = max(xinfo%lat, -89.995)
 684    xinfo%lat = min(xinfo%lat,  89.995)
 685    call da_llxy(xinfo, xloc(ii), outside, outside_all)
 686    if ( outside_all ) then
 687       cycle check_loop_2
 688    end if
 689    inside_domain(ii) = .true.
 690    ifgat = iwindow(ii)
 691    if ( ifgat < 0 ) cycle check_loop_2
 692    ntotal = ntotal + 1
 693    ntotal_ifgat(ifgat) = ntotal_ifgat(ifgat) + 1
 694    if ( outside ) then
 695       cycle check_loop_2
 696    end if
 697    inside_patch(ii) = .true.
 698    nlocal = nlocal + 1
 699
 700 end do check_loop_2
 701
 702 ! transfer data arrays to iv structures
 703
 704 iv%info(polaramv)%ntotal = ntotal
 705 iv%info(polaramv)%nlocal = nlocal
 706 max_lev = 1
 707 iv%info(polaramv)%max_lev = max_lev
 708
 709 !allocate iv and iv%info
 710 if ( nlocal > 0 ) then
 711    allocate (iv%polaramv(1:nlocal))
 712    allocate (iv%info(polaramv)%name(nlocal))
 713    allocate (iv%info(polaramv)%platform(nlocal))
 714    allocate (iv%info(polaramv)%id(nlocal))
 715    allocate (iv%info(polaramv)%date_char(nlocal))
 716    allocate (iv%info(polaramv)%levels(nlocal))
 717    allocate (iv%info(polaramv)%lat(max_lev,nlocal))
 718    allocate (iv%info(polaramv)%lon(max_lev,nlocal))
 719    allocate (iv%info(polaramv)%elv(nlocal))
 720    allocate (iv%info(polaramv)%pstar(nlocal))
 721    allocate (iv%info(polaramv)%slp(nlocal))
 722    allocate (iv%info(polaramv)%pw(nlocal))
 723    allocate (iv%info(polaramv)%x  (kms:kme,nlocal))
 724    allocate (iv%info(polaramv)%y  (kms:kme,nlocal))
 725    allocate (iv%info(polaramv)%i  (kms:kme,nlocal))
 726    allocate (iv%info(polaramv)%j  (kms:kme,nlocal))
 727    allocate (iv%info(polaramv)%dx (kms:kme,nlocal))
 728    allocate (iv%info(polaramv)%dxm(kms:kme,nlocal))
 729    allocate (iv%info(polaramv)%dy (kms:kme,nlocal))
 730    allocate (iv%info(polaramv)%dym(kms:kme,nlocal))
 731    allocate (iv%info(polaramv)%k  (max_lev,nlocal))
 732    allocate (iv%info(polaramv)%dz (max_lev,nlocal))
 733    allocate (iv%info(polaramv)%dzm(max_lev,nlocal))
 734    allocate (iv%info(polaramv)%zk (max_lev,nlocal))
 735    allocate (iv%info(polaramv)%proc_domain(max_lev,nlocal))
 736    allocate (iv%info(polaramv)%thinned(max_lev,nlocal))
 737    allocate (iv%info(polaramv)%obs_global_index(nlocal))
 738    iv%info(polaramv)%elv(:)            = missing_r
 739    iv%info(polaramv)%proc_domain(:,:)  = .false.
 740    iv%info(polaramv)%thinned(:,:)      = .false.
 741    iv%info(polaramv)%zk(:,:)           = missing_r
 742    iv%info(polaramv)%slp(:) = field_type(missing_r, missing_data, xmiss, missing_r, missing_r)
 743    iv%info(polaramv)%pw(:)  = field_type(missing_r, missing_data, xmiss, missing_r, missing_r)
 744
 745    nlev = 1
 746    do ilocal = 1, nlocal
 747       allocate (iv%polaramv(ilocal)%p(1:nlev))
 748       allocate (iv%polaramv(ilocal)%u(1:nlev))
 749       allocate (iv%polaramv(ilocal)%v(1:nlev))
 750       iv%polaramv(ilocal)%p(:) = missing_r
 751       iv%polaramv(ilocal)%u(:) = field_type(missing_r, missing_data, xmiss, missing_r, missing_r)
 752       iv%polaramv(ilocal)%v(:) = field_type(missing_r, missing_data, xmiss, missing_r, missing_r)
 753    end do
 754 end if
 755
 756 ilocal = 0
 757 itotal = 0
 758
 759 do kk = 1, num_fgat_time
 760
 761    iv%info(polaramv)%ptotal(kk)=0
 762
 763    nobs_loop_2: do ii=1, nobs
 764
 765       if ( .not. inside_domain(ii) ) cycle nobs_loop_2
 766       itotal = itotal + 1
 767
 768       if ( .not. inside_patch(ii)  ) cycle nobs_loop_2
 769
 770       if ( iwindow(ii) /= kk ) then
 771          cycle nobs_loop_2
 772       else
 773          ilocal = ilocal+1
 774          iv%info(polaramv)%platform(ilocal) = 'FM-88 SATWND'
 775          if ( thin_conv_opt(polaramv) == thin_superob .or. &
 776               thin_conv_opt(polaramv) == thin_superob_hv ) then
 777             iv%info(polaramv)%name(ilocal) = 'satwnd'
 778             iv%info(polaramv)%id(ilocal)   = 'SUPEROB'
 779          else
 780             write(iv%info(polaramv)%name(ilocal), '(a,a,i3)') &
 781                trim(str_name(ii)),' ',rptype(ii)
 782             iv%info(polaramv)%id(ilocal)        = str_stid(ii)
 783          end if
 784          iv%info(polaramv)%date_char(ilocal) = date_char(ii)
 785          iv%info(polaramv)%levels(ilocal)    = nlev  ! = 1
 786          iv%info(polaramv)%lat(:,ilocal)     = smeta(ii,idx_lat)
 787          iv%info(polaramv)%lon(:,ilocal)     = smeta(ii,idx_lon)
 788
 789          iv%info(polaramv)%x(:,ilocal)       = xloc(ii)%x
 790          iv%info(polaramv)%y(:,ilocal)       = xloc(ii)%y
 791          iv%info(polaramv)%i(:,ilocal)       = xloc(ii)%i
 792          iv%info(polaramv)%j(:,ilocal)       = xloc(ii)%j
 793          iv%info(polaramv)%dx(:,ilocal)      = xloc(ii)%dx
 794          iv%info(polaramv)%dxm(:,ilocal)     = xloc(ii)%dxm
 795          iv%info(polaramv)%dy(:,ilocal)      = xloc(ii)%dy
 796          iv%info(polaramv)%dym(:,ilocal)     = xloc(ii)%dym
 797
 798          iv%info(polaramv)%obs_global_index(ilocal) = itotal
 799
 800          iv%polaramv(ilocal)%p(1) = smeta(ii,idx_p)
 801
 802          uob  = sdata_val(ii,idx_u)
 803          if ( abs(uob - missing_r) > 1.0 ) then
 804             iv%polaramv(ilocal)%u(1)%inv = uob
 805             iv%polaramv(ilocal)%u(1)%qc  = 0
 806             err_val = sdata_err(ii,idx_u)
 807             if ( err_val > 0.0 ) then
 808                iv%polaramv(ilocal)%u(1)%error = err_val
 809             end if
 810          end if
 811          vob  = sdata_val(ii,idx_v)
 812          if ( abs(vob - missing_r) > 1.0 ) then
 813             iv%polaramv(ilocal)%v(1)%inv = vob
 814             iv%polaramv(ilocal)%v(1)%qc  = 0
 815             err_val = sdata_err(ii,idx_v)
 816             if ( err_val > 0.0 ) then
 817                iv%polaramv(ilocal)%v(1)%error = err_val
 818             end if
 819          end if
 820
 821          if ( uv_error_opt(polaramv) == error_opt_nml ) then
 822             iv%polaramv(ilocal)%u(1)%error = uv_error_val(polaramv)
 823             iv%polaramv(ilocal)%v(1)%error = uv_error_val(polaramv)
 824          end if
 825
 826          ! u/v-earth to u/v-grid
 827          if ( abs(uob - missing_r) > 1.0 .and. &
 828               abs(vob - missing_r) > 1.0 ) then
 829             call da_earth_2_model_wind ( uob, vob, u_grid, v_grid, &
 830                                          iv%info(polaramv)%lon(1,ilocal) )
 831             iv%polaramv(ilocal)%u(1)%inv = u_grid
 832             iv%polaramv(ilocal)%v(1)%inv = v_grid
 833          end if
 834
 835       end if ! inside time window
 836    end do nobs_loop_2
 837
 838    ntotal_ifgat(kk) = ntotal_ifgat(kk) + ntotal_ifgat(kk-1)
 839    iv%info(polaramv)%ptotal(kk) = ntotal_ifgat(kk)
 840    iv%info(polaramv)%plocal(kk) = ilocal
 841 end do ! num_fgat_time
 842
 843 deallocate (inside_patch)
 844 deallocate (inside_domain)
 845 deallocate (xloc)
 846 if ( thin_conv_opt(polaramv) == no_thin .or. &
 847      thin_conv_opt(polaramv) == thin_multi ) then
 848    deallocate (rptype)
 849 end if
 850 deallocate (iwindow)
 851 deallocate (smeta)
 852 deallocate (sdata_err)
 853 deallocate (sdata_val)
 854 deallocate (str_stid)
 855 deallocate (str_name)
 856
 857 ! transfer iv to a smaller ob structure
 858
 859 ob % nlocal(polaramv) = iv%info(polaramv)%nlocal
 860 ob % ntotal(polaramv) = iv%info(polaramv)%ntotal
 861 if ( nlocal > 0) then
 862    allocate (ob % polaramv(1:nlocal))
 863    do n = 1, nlocal
 864       nlevel = iv%info(polaramv)%levels(n)
 865       allocate (ob % polaramv(n)%u(1:nlevel))
 866       allocate (ob % polaramv(n)%v(1:nlevel))
 867       do k = 1, nlevel
 868          ob % polaramv(n) % u(k) = iv % polaramv(n) % u(k) % inv
 869          ob % polaramv(n) % v(k) = iv % polaramv(n) % v(k) % inv
 870       end do
 871    end do
 872 end if
 873
 874 contains
 875
 876 !--------------------------------------------------------------
 877
 878 subroutine filter_obs_satwnd(qm, rptype, prs, satzen, landsea, pccf1, pccf2, cvwd, wspd)
 879
 880 ! refer to GSI/read_satwnd.f90
 881
 882 implicit none
 883
 884 integer(i_kind), intent(out) :: qm
 885 integer(i_kind), intent(in)  :: rptype
 886 real(r_kind),    intent(in)  :: prs
 887 real(r_kind),    intent(in)  :: satzen
 888 real(r_kind),    intent(in)  :: landsea
 889 real(r_kind),    intent(in)  :: pccf1
 890 real(r_kind),    intent(in)  :: pccf2
 891 real(r_kind),    intent(in)  :: cvwd
 892 real(r_kind),    intent(in)  :: wspd
 893
 894 real(r_kind) :: experr_norm
 895 logical :: EC_AMV_QC = .true.
 896 integer(i_kind) :: iland = 0
 897
 898 qm = 2  ! neutral or not checked
 899
 900 if ( satzen > 68.0_r_kind ) qm = 15
 901
 902 if ( rptype == 260 ) then
 903    if ( pccf1 < 85.0_r_kind ) qm = 15
 904 end if
 905
 906 if ( rptype == 240 .or. &
 907      rptype == 245 .or. &
 908      rptype == 246 .or. &
 909      rptype == 247 .or. &
 910      rptype == 251 ) then
 911
 912    if ( pccf1 < 80.0_r_kind .or. pccf1 > 100.0_r_kind ) qm = 15 ! reject data with low QI
 913
 914    if ( prs < 12500.0_r_kind ) qm = 15 ! reject data above 125hPa
 915
 916    experr_norm = 10.0_r_kind - 0.1_r_kind * pccf2
 917    if ( wspd > 0.1_r_kind ) then
 918       experr_norm = experr_norm / wspd
 919    else
 920       experr_norm = 100.0_r_kind
 921    end if
 922    if ( experr_norm > 0.9_r_kind ) qm = 15 ! reject data with estimated error/spd>0.9
 923
 924    if ( rptype == 240 .or. &
 925         rptype == 245 .or. &
 926         rptype == 246 .or. &
 927         rptype == 251 ) then
 928       if ( cvwd < 0.04_r_kind ) qm = 15
 929       if ( cvwd > 0.50_r_kind ) qm = 15
 930    end if
 931
 932    if ( EC_AMV_QC ) then
 933       if ( pccf1 < 90_r_kind .or. pccf1 > 100.0_r_kind ) qm = 15 ! stricter QI
 934       if ( prs < 15000.0_r_kind) qm = 15 ! all high level
 935       if ( rptype == 251 .and. prs < 70000.0_r_kind ) qm = 15  ! VIS
 936       if ( rptype == 246 .and. prs > 30000.0_r_kind ) qm = 15  ! WVCA
 937       if ( nint(landsea) == iland .and. prs > 85000.0_r_kind) qm = 15  ! low over land
 938    end if
 939 end if  ! rptype 240, 245, 246, 247, 251
 940
 941 end subroutine filter_obs_satwnd
 942
 943 !--------------------------------------------------------------
 944
 945 subroutine fill_datalink (datalink, rfill, ifill)
 946
 947 implicit none
 948
 949 type (datalink_satwnd), intent(inout) :: datalink
 950 real(r_kind),           intent(in)    :: rfill    ! fill value in real
 951 integer(i_kind),        intent(in)    :: ifill    ! fill value in integer
 952
 953 integer(i_kind) :: i
 954
 955 datalink % lat      = rfill
 956 datalink % lon      = rfill
 957 datalink % satid    = ifill
 958 datalink % itg      = ifill
 959 datalink % ifgat    = ifill
 960 datalink % landsea  = rfill
 961 datalink % satzen   = rfill
 962 datalink % prs      = rfill
 963 datalink % wdir     = rfill
 964 datalink % wspd     = rfill
 965 datalink % cvwd     = rfill
 966 datalink % uwind    = rfill
 967 datalink % vwind    = rfill
 968 datalink % pccf1    = rfill
 969 datalink % pccf2    = rfill
 970 datalink % err      = rfill
 971 datalink % qm       = ifill
 972
 973 end subroutine fill_datalink
 974
 975 !--------------------------------------------------------------
 976
 977 subroutine set_rptype_satwnd(subset, satid, rptype, stid, name)
 978
 979 implicit none
 980 character(len=*), intent(in)  :: subset  ! bufr message subset
 981 integer(i_kind),  intent(in)  :: satid   ! satellite id
 982 integer(i_kind),  intent(out) :: rptype  ! report type
 983 character(len=nstring), intent(out) :: stid    ! station id
 984 character(len=nstring), intent(out) :: name    ! instrument name
 985
 986 character(len=3) :: csatid
 987
 988 rptype = -1
 989 stid   = ''
 990 write(csatid, '(i3.3)') satid
 991
 992 select case ( trim(subset) )
 993 case ( 'NC005030' )  ! GOES IR LW
 994    rptype = 245
 995    stid = 'IR'//csatid
 996    name = 'GOES IR LW'
 997 case ( 'NC005039' )  ! GOES IR SW
 998    rptype = 240
 999    stid = 'IR'//csatid
1000    name = 'GOES IR SW'
1001 case ( 'NC005032' )  ! GOES VIS
1002    rptype = 251
1003    stid = 'VI'//csatid
1004    name = 'GOES VIS'
1005 case ( 'NC005034' )  ! GOES WV cloud top
1006    rptype = 246
1007    stid = 'WV'//csatid
1008    name = 'GOES WV cloud top'
1009 case ( 'NC005031' )  ! GOES WV clear sky/deep layer
1010    rptype = 247
1011    stid = 'WV'//csatid
1012    name = 'GOES WV clear sky/depp layer'
1013 case ( 'NC005080', 'NC005081' )  ! AVHRR (METOP-a/b/c, NOAA-15/18/19)
1014    rptype = 244
1015    stid = 'IR'//csatid
1016    name = 'AVHRR'
1017 case ( 'NC005090', 'NC005091' )  ! VIIRS (NPP, NOAA-20)
1018    rptype = 260
1019    stid = 'IR'//csatid
1020    name = 'VIIRS'
1021 case ( 'NC005072' )  ! LEOGEO (non-specific mixture of geostationary and low earth orbiting satellites)
1022    rptype = 255
1023    stid = 'IR'//csatid  ! satid=854
1024    name = 'LEOGEO'
1025 end select
1026
1027 end subroutine set_rptype_satwnd
1028
1029 !--------------------------------------------------------------
1030
1031 integer function ufo_vars_getindex(vars, varname)
1032 implicit none
1033 character(len=*), intent(in) :: vars(:)
1034 character(len=*), intent(in) :: varname
1035
1036 integer :: ivar
1037
1038 ufo_vars_getindex = -1
1039
1040 do ivar = 1, size(vars)
1041   if (trim(vars(ivar)) == trim(varname)) then
1042     ufo_vars_getindex = ivar
1043     exit
1044   endif
1045 enddo
1046
1047 end function ufo_vars_getindex
1048
1049 !--------------------------------------------------------------
1050
1051 real function get_zk (nlev, levels, val)
1052
1053 implicit none
1054 integer, intent(in)  :: nlev
1055 real,    intent(in)  :: levels(nlev)
1056 real,    intent(in)  :: val
1057 integer              :: k
1058 integer              :: level_index
1059 logical              :: decreasing
1060
1061 if ( levels(nlev) < levels(1) ) then
1062    decreasing = .true.
1063 else
1064    decreasing = .false.
1065 end if
1066
1067 level_index = 0 ! initialize
1068
1069 if ( decreasing ) then
1070    if ( val >= levels(1) ) then
1071       level_index = 1
1072    else if ( val < levels(nlev-1) .and. val >= levels(nlev) ) then
1073       level_index = nlev
1074    else
1075       do k = 2, nlev-1
1076          if ( val >= levels(k) ) then
1077             level_index = k
1078             exit
1079          end if
1080       end do
1081    end if
1082 else
1083    if ( val <= levels(1) ) then
1084       level_index = 1
1085    else if ( val > levels(nlev-1) .and. val <= levels(nlev) ) then
1086       level_index = nlev
1087    else
1088       do k = 2, nlev-1
1089          if ( val <= levels(k) ) then
1090             level_index = k
1091             exit
1092          end if
1093       end do
1094    end if
1095 end if ! decreasing or not
1096
1097 get_zk = real(level_index) + &
1098          (val-levels(level_index))/(levels(level_index+1)-levels(level_index))
1099
1100 end function get_zk
1101
1102 end subroutine da_read_obs_bufr_satwnd