var/da/da_obs_io/da_read_lsac_util.inc

   1 subroutine da_read_lsac_wrfinput(iv, onlyscan)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: extract pseduo-observations from wrfinput and use them as bogus
   5    ! observations to constrain the large-scale pattern of WRFDA analysis
   6    !
   7    !    Author: Xiaowen Tang,  NJU,  Date: 5/2/2016
   8    !    Update 7/15/2017: fix the problem of staggering variables
   9    !---------------------------------------------------------------------
  10
  11 implicit none
  12
  13 type (iv_type), intent(inout)       ::  iv
  14 logical       , intent(in)          ::  onlyscan
  15
  16 type (multi_level_type)             ::  platform
  17 logical                             ::  outside, outside_all
  18
  19 integer                             ::  i, j, k, ki, ndims, nrecs, nlocal, iunit, nlevels, ilevel
  20 integer                             ::  u_qc, v_qc, t_qc, q_qc
  21 character(len=512)                  ::  lsac_wrfinput
  22 integer, dimension(4)               ::  dims_u, dims_v, dims_t, dims_p, dims_pb, dims_q
  23 integer, dimension(4)               ::  dims_lat, dims_lon, dims_alt, dims_phb, dims_ph
  24 real, allocatable, dimension(:,:,:) ::  u_lsac, v_lsac, w_lsac, theta_lsac, p_lsac, q_lsac, tk_lsac
  25 real, allocatable, dimension(:,:,:) ::  pb_lsac, ph_lsac, phb_lsac, height_lsac, press
  26 real, allocatable, dimension(:,:)   ::  lat_lsac, lon_lsac
  27 logical                             ::  debug
  28 logical, external                   ::  wrf_dm_on_monitor
  29 logical                             ::  isfile
  30
  31
  32 if (trace_use) call da_trace_entry("da_read_lsac_util")
  33
  34 lsac_wrfinput  = 'fg_l'
  35 debug=.false.
  36
  37 inquire(file=trim(lsac_wrfinput), exist=isfile)
  38 if ( .not. isfile ) then
  39    write(unit=message(1),fmt='(a,a,a)') 'File ',trim(lsac_wrfinput),' for LSAC is missing.'
  40    call da_error(__FILE__,__LINE__,message(1:1))
  41 endif
  42
  43 if (onlyscan) then
  44    if (rootproc) then
  45       call da_get_dims_cdf( lsac_wrfinput, 'XLAT',   dims_lat,  ndims, debug)
  46       call da_get_dims_cdf( lsac_wrfinput, 'XLONG',  dims_lon,  ndims, debug)
  47       call da_get_dims_cdf( lsac_wrfinput, 'T',      dims_t,    ndims, debug)
  48    endif
  49    call wrf_dm_bcast_integer(dims_lon, 4)
  50    call wrf_dm_bcast_integer(dims_lat, 4)
  51    call wrf_dm_bcast_integer(dims_t,   4)
  52 #ifdef DEBUG_LSAC
  53    print *, '****SCAN LSAC lat_dims: ', myproc, dims_lat(1), dims_lat(2), dims_lat(3), dims_lat(4)
  54    print *, '****SCAN LSAC lon_dims: ', myproc, dims_lon(1), dims_lon(2), dims_lon(3), dims_lon(4)
  55    print *, '****SCAN LSAC t_dims: ', myproc, dims_t(1), dims_t(2), dims_t(3), dims_t(4)
  56 #endif
  57    allocate(lat_lsac(dims_lat(1), dims_lat(2)))
  58    allocate(lon_lsac(dims_lon(1), dims_lon(2)))
  59
  60    nlevels = (dims_t(3)-lsac_nv_start+1)/lsac_nv_step
  61
  62    !---------------------------------------------------------
  63    ! Reading data from WRF Input file
  64    !---------------------------------------------------------
  65    if (rootproc) then
  66       call da_get_var_2d_real_cdf( lsac_wrfinput, 'XLAT',   lat_lsac, dims_lat(1), dims_lat(2), 1, debug)
  67       call da_get_var_2d_real_cdf( lsac_wrfinput, 'XLONG',  lon_lsac, dims_lon(1), dims_lon(2), 1, debug)
  68    endif
  69    call wrf_dm_bcast_real(lat_lsac, dims_lat(1)*dims_lat(2))
  70    call wrf_dm_bcast_real(lon_lsac, dims_lon(1)*dims_lon(2))
  71
  72    !Assigning max_lev and counts in the iv type in onlyscan mode
  73    nlocal=0
  74    nrecs =0
  75    do i=1, dims_lon(1), lsac_nh_step
  76       do j=1, dims_lat(2), lsac_nh_step
  77          platform%info%lat      = lat_lsac(i,j)
  78          platform%info%lon      = lon_lsac(i,j)
  79          platform%info%elv      = 0.0
  80          platform%info%name     = 'LSAC'
  81          platform%info%platform = 'FM-???  LSAC'
  82          platform%info%id       = '?????'
  83          platform%info%date_char= '????-??-??_??:??:??'
  84          platform%info%pstar    = 0.D0
  85          platform%info%levels   = nlevels
  86          if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
  87          if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
  88             platform%info%lon = 0.0
  89          endif
  90
  91          call da_llxy (platform%info, platform%loc, outside, outside_all)
  92          if (.not.outside) then
  93             nlocal = nlocal+1
  94          endif
  95          if (.not.outside_all) then
  96             nrecs  = nrecs+1
  97          endif
  98       enddo
  99    enddo
 100 #ifdef DEBUG_LSAC
 101    print *, '******SCAN LSAC: ', myproc, nlevels, nlocal, nrecs
 102 #endif
 103    iv%info(bogus)%max_lev = nlevels
 104    iv%info(bogus)%nlocal  = nlocal
 105    iv%info(bogus)%ntotal  = nrecs
 106    deallocate(lat_lsac)
 107    deallocate(lon_lsac)
 108 else
 109    !---------------------------------------------------------
 110    ! Getting information from NETCDF files (WRF Input file)
 111    !---------------------------------------------------------
 112    if (rootproc) then
 113       call da_get_dims_cdf( lsac_wrfinput, 'XLAT',   dims_lat,  ndims, debug)
 114       call da_get_dims_cdf( lsac_wrfinput, 'XLONG',  dims_lon,  ndims, debug)
 115       call da_get_dims_cdf( lsac_wrfinput, 'U',      dims_u,    ndims, debug)
 116       call da_get_dims_cdf( lsac_wrfinput, 'V',      dims_v,    ndims, debug)
 117       call da_get_dims_cdf( lsac_wrfinput, 'T',      dims_t,    ndims, debug)
 118       call da_get_dims_cdf( lsac_wrfinput, 'PB',     dims_pb,   ndims, debug)
 119       call da_get_dims_cdf( lsac_wrfinput, 'P',      dims_p,    ndims, debug)
 120       call da_get_dims_cdf( lsac_wrfinput, 'QVAPOR', dims_q,    ndims, debug)
 121       call da_get_dims_cdf( lsac_wrfinput, 'PHB',    dims_phb,  ndims, debug)
 122       call da_get_dims_cdf( lsac_wrfinput, 'PH',     dims_ph,   ndims, debug)
 123    endif
 124
 125 #ifdef DEBUG_LSAC
 126    print *, 'lat_dims: ', myproc, dims_lat(1), dims_lat(2), dims_lat(3), dims_lat(4)
 127    print *, 'lon_dims: ', myproc, dims_lon(1), dims_lon(2), dims_lon(3), dims_lon(4)
 128    print *, 'u_dims: ', myproc, dims_u(1), dims_u(2), dims_u(3), dims_u(4)
 129    print *, 'v_dims: ', myproc, dims_v(1), dims_v(2), dims_v(3), dims_v(4)
 130    print *, 't_dims: ', myproc, dims_t(1), dims_t(2), dims_t(3), dims_t(4)
 131    print *, 'p_dims: ', myproc, dims_p(1), dims_p(2), dims_p(3), dims_p(4)
 132    print *, 'pb_dims: ', myproc, dims_pb(1), dims_pb(2), dims_pb(3), dims_pb(4)
 133    print *, 'q_dims: ', myproc, dims_q(1), dims_q(2), dims_q(3), dims_q(4)
 134    print *, 'ph_dims: ', myproc, dims_ph(1), dims_ph(2), dims_ph(3), dims_ph(4)
 135    print *, 'phb_dims: ', myproc, dims_phb(1), dims_phb(2), dims_phb(3), dims_phb(4)
 136 #endif
 137
 138    call wrf_dm_bcast_integer(dims_u,   4)
 139    call wrf_dm_bcast_integer(dims_v,   4)
 140    call wrf_dm_bcast_integer(dims_t,   4)
 141    call wrf_dm_bcast_integer(dims_p,   4)
 142    call wrf_dm_bcast_integer(dims_pb,  4)
 143    call wrf_dm_bcast_integer(dims_q,   4)
 144    call wrf_dm_bcast_integer(dims_lat, 4)
 145    call wrf_dm_bcast_integer(dims_lon, 4)
 146    call wrf_dm_bcast_integer(dims_phb, 4)
 147    call wrf_dm_bcast_integer(dims_ph,  4)
 148
 149    nrecs   = (1 + (dims_lat(1) - 1)/lsac_nh_step) * (1 + (dims_lat(2) - 1) / lsac_nh_step)
 150    nlevels = (dims_t(3) - lsac_nv_start + 1) / lsac_nv_step
 151
 152    !---------------------------------------------------------
 153    ! Allocating memory
 154    !---------------------------------------------------------
 155    allocate(lat_lsac(dims_lat(1),   dims_lat(2)))
 156    allocate(lon_lsac(dims_lon(1),   dims_lon(2)))
 157    allocate(u_lsac(dims_u(1),       dims_u(2),   dims_u(3)   ))
 158    allocate(v_lsac(dims_v(1),       dims_v(2),   dims_v(3)   ))
 159    allocate(theta_lsac(dims_t(1),   dims_t(2),   dims_t(3)   ))
 160    allocate(p_lsac(dims_p(1),       dims_p(2),   dims_p(3)   ))
 161    allocate(pb_lsac(dims_pb(1),     dims_pb(2),  dims_pb(3)  ))
 162    allocate(q_lsac(dims_q(1),       dims_q(2),   dims_q(3)   ))
 163    allocate(phb_lsac(dims_phb(1),   dims_phb(2), dims_phb(3) ))
 164    allocate(ph_lsac(dims_ph(1),     dims_ph(2),  dims_ph(3)  ))
 165    allocate(height_lsac(dims_ph(1), dims_ph(2),  dims_ph(3)  ))
 166    allocate(press(dims_p(1),        dims_p(2),   dims_p(3)   ))
 167    allocate(tk_lsac(dims_t(1),      dims_t(2),   dims_t(3)   ))
 168
 169    !---------------------------------------------------------
 170    ! Reading data from WRF Input file
 171    !---------------------------------------------------------
 172    if (rootproc) then
 173       call da_get_var_2d_real_cdf( lsac_wrfinput, 'XLAT',   lat_lsac,  dims_lat(1), dims_lat(2), 1, debug)
 174       call da_get_var_2d_real_cdf( lsac_wrfinput, 'XLONG',  lon_lsac,  dims_lon(1), dims_lon(2), 1, debug)
 175       call da_get_var_3d_real_cdf( lsac_wrfinput, 'U',      u_lsac,    dims_u(1),   dims_u(2),   dims_u(3),   1, debug)
 176       call da_get_var_3d_real_cdf( lsac_wrfinput, 'V',      v_lsac,    dims_v(1),   dims_v(2),   dims_v(3),   1, debug)
 177       call da_get_var_3d_real_cdf( lsac_wrfinput, 'T',      theta_lsac,dims_t(1),   dims_t(2),   dims_t(3),   1, debug)
 178       call da_get_var_3d_real_cdf( lsac_wrfinput, 'P',      p_lsac,    dims_p(1),   dims_p(2),   dims_p(3),   1, debug)
 179       call da_get_var_3d_real_cdf( lsac_wrfinput, 'PB',     pb_lsac,   dims_p(1),   dims_p(2),   dims_p(3),   1, debug)
 180       call da_get_var_3d_real_cdf( lsac_wrfinput, 'QVAPOR', q_lsac,    dims_q(1),   dims_q(2),   dims_q(3),   1, debug)
 181       call da_get_var_3d_real_cdf( lsac_wrfinput, 'PHB',    phb_lsac,  dims_phb(1), dims_phb(2), dims_phb(3), 1, debug)
 182       call da_get_var_3d_real_cdf( lsac_wrfinput, 'PH',     ph_lsac,   dims_ph(1),  dims_ph(2),  dims_ph(3),  1, debug)
 183    endif
 184    call wrf_dm_bcast_real(lat_lsac,   PRODUCT(dims_lat(1:2)))
 185    call wrf_dm_bcast_real(lon_lsac,   PRODUCT(dims_lon(1:2)))
 186    call wrf_dm_bcast_real(u_lsac,     PRODUCT(dims_u(1:3))  )
 187    call wrf_dm_bcast_real(v_lsac,     PRODUCT(dims_v(1:3))  )
 188    call wrf_dm_bcast_real(theta_lsac, PRODUCT(dims_t(1:3))  )
 189    call wrf_dm_bcast_real(p_lsac,     PRODUCT(dims_p(1:3))  )
 190    call wrf_dm_bcast_real(pb_lsac,    PRODUCT(dims_p(1:3))  )
 191    call wrf_dm_bcast_real(q_lsac,     PRODUCT(dims_q(1:3))  )
 192    call wrf_dm_bcast_real(phb_lsac,   PRODUCT(dims_phb(1:3)))
 193    call wrf_dm_bcast_real(ph_lsac,    PRODUCT(dims_ph(1:3)) )
 194
 195    !---------------------------------------------------------
 196    !Calculating the height
 197    !---------------------------------------------------------
 198    height_lsac = (phb_lsac + ph_lsac) / gravity
 199    press = (p_lsac + pb_lsac) * 0.01
 200
 201    !Temperature from potential temperature
 202    tk_lsac = (t0 + theta_lsac) * ((press / 1000.0) ** (gas_constant / cp))
 203
 204    if (lsac_print_details .and. rootproc) then
 205       call da_get_unit(iunit)
 206       open(iunit, file='lsac_details')
 207    endif
 208
 209    ! Assigning errors, heights and inv in the iv type
 210    nlocal=0
 211    do i=1, dims_lon(1), lsac_nh_step
 212       do j=1, dims_lat(2), lsac_nh_step
 213          ilevel = 0
 214          do k=lsac_nv_start, dims_t(3), lsac_nv_step
 215             ilevel = ilevel+1
 216             if (lsac_use_u) then
 217                u_qc   = 0
 218             else
 219                u_qc   = missing_data
 220             endif
 221             if (lsac_use_v) then
 222                v_qc   = 0
 223             else
 224                v_qc   = missing_data
 225             endif
 226             if (lsac_use_t) then
 227                t_qc   = 0
 228             else
 229                t_qc   = missing_data
 230             endif
 231             if (lsac_use_q) then
 232                q_qc   = 0
 233             else
 234                q_qc   = missing_data
 235             endif
 236
 237             platform%each(ilevel)%height= (height_lsac(i,j,k)+height_lsac(i,j,k+1))/2.
 238
 239             platform%each(ilevel)%u%inv= (u_lsac(i,j,k)+u_lsac(i+1,j,k))/2.
 240             platform%each(ilevel)%u%error=lsac_u_error
 241             platform%each(ilevel)%u%qc=u_qc
 242
 243             platform%each(ilevel)%v%inv= (v_lsac(i,j,k)+v_lsac(i,j+1,k))/2.
 244             platform%each(ilevel)%v%error=lsac_v_error
 245             platform%each(ilevel)%v%qc=v_qc
 246
 247             platform%each(ilevel)%t%inv=tk_lsac(i,j,k)
 248             platform%each(ilevel)%t%error=lsac_t_error
 249             platform%each(ilevel)%t%qc=t_qc
 250
 251             platform%each(ilevel)%q%inv=q_lsac(i,j,k)
 252             platform%each(ilevel)%q%error=lsac_q_error
 253             platform%each(ilevel)%q%qc=q_qc
 254
 255             if(lsac_print_details .and. rootproc) then
 256                write(iunit,'(3I5,3f10.3,x,4(f10.3,x,f10.3,x,i4))') i, j, k, &
 257                  (height_lsac(i,j,k) + height_lsac(i,j,k+1))/2., lat_lsac(i,j), lon_lsac(i,j), &
 258                  (u_lsac(i,j,k)+u_lsac(i+1,j,k))/2.,   lsac_u_error     ,   u_qc, &
 259                  (v_lsac(i,j,k)+v_lsac(i,j+1,k))/2.,   lsac_v_error     ,   v_qc, &
 260                  tk_lsac(i,j,k)                    ,   lsac_t_error     ,   t_qc, &
 261                  q_lsac(i,j,k)*1000.               ,   lsac_q_error*1000.,  q_qc
 262             endif
 263          enddo
 264          platform%info%lat      = lat_lsac(i,j)
 265          platform%info%lon      = lon_lsac(i,j)
 266          platform%info%elv      = 0.
 267          platform%info%name     = 'LSAC'
 268          platform%info%platform = 'FM-???  LSAC'
 269          platform%info%id       = '?????'
 270          platform%info%date_char= '????-??-??_??:??:??'
 271          platform%info%pstar    = 0.D0
 272          platform%info%levels   = nlevels
 273          if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 274          if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 275             platform%info%lon = 0.0
 276          endif
 277
 278          call da_llxy (platform%info, platform%loc, outside, outside_all)
 279          if (outside) then
 280             cycle
 281          endif
 282          nlocal = nlocal+1
 283
 284          if (nlevels > 0) then
 285             allocate (iv%bogus(nlocal)%h (1:nlevels))
 286             allocate (iv%bogus(nlocal)%p (1:nlevels))
 287             allocate (iv%bogus(nlocal)%u (1:nlevels))
 288             allocate (iv%bogus(nlocal)%v (1:nlevels))
 289             allocate (iv%bogus(nlocal)%t (1:nlevels))
 290             allocate (iv%bogus(nlocal)%q (1:nlevels))
 291             do ki = 1, nlevels
 292                iv%bogus(nlocal)%h(ki) = platform%each(ki)%height
 293                iv%bogus(nlocal)%p(ki) = missing_r
 294                iv%bogus(nlocal)%u(ki) = platform%each(ki)%u
 295                iv%bogus(nlocal)%v(ki) = platform%each(ki)%v
 296                iv%bogus(nlocal)%t(ki) = platform%each(ki)%t
 297                iv%bogus(nlocal)%q(ki) = platform%each(ki)%q
 298                iv%bogus(nlocal)%slp%inv   = 0.
 299                iv%bogus(nlocal)%slp%qc    = missing_data
 300                iv%bogus(nlocal)%slp%error = missing_r
 301             end do
 302          else
 303             nlevels = 1
 304             allocate (iv%bogus(nlocal)%h (1:nlevels))
 305             allocate (iv%bogus(nlocal)%p (1:nlevels))
 306             allocate (iv%bogus(nlocal)%u (1:nlevels))
 307             allocate (iv%bogus(nlocal)%v (1:nlevels))
 308             allocate (iv%bogus(nlocal)%t (1:nlevels))
 309             allocate (iv%bogus(nlocal)%q (1:nlevels))
 310             iv%bogus(nlocal)%h       = missing_r
 311             iv%bogus(nlocal)%p       = missing_r
 312             iv%bogus(nlocal)%u%inv   = missing_r
 313             iv%bogus(nlocal)%u%qc    = missing
 314             iv%bogus(nlocal)%u%error = abs(missing_r)
 315             iv%bogus(nlocal)%v       = iv%bogus(nlocal)%u
 316             iv%bogus(nlocal)%t       = iv%bogus(nlocal)%u
 317             iv%bogus(nlocal)%q       = iv%bogus(nlocal)%u
 318          end if
 319
 320          iv%info(bogus)%name(nlocal)          = platform%info%name
 321          iv%info(bogus)%platform(nlocal)      = platform%info%platform
 322          iv%info(bogus)%id(nlocal)            = platform%info%id
 323          iv%info(bogus)%date_char(nlocal)     = platform%info%date_char
 324          iv%info(bogus)%levels(nlocal)        = platform%info%levels
 325          iv%info(bogus)%lat(:,nlocal)         = platform%info%lat
 326          iv%info(bogus)%lon(:,nlocal)         = platform%info%lon
 327          iv%info(bogus)%elv(nlocal)           = platform%info%elv
 328          iv%info(bogus)%pstar(nlocal)         = platform%info%pstar
 329          iv%info(bogus)%max_lev               = platform%info%levels
 330
 331          iv%info(bogus)%slp(nlocal)           = platform%loc%slp
 332          iv%info(bogus)%pw(nlocal)            = platform%loc%pw
 333          iv%info(bogus)%x(:,nlocal)           = platform%loc%x
 334          iv%info(bogus)%y(:,nlocal)           = platform%loc%y
 335          iv%info(bogus)%i(:,nlocal)           = platform%loc%i
 336          iv%info(bogus)%j(:,nlocal)           = platform%loc%j
 337          iv%info(bogus)%dx(:,nlocal)          = platform%loc%dx
 338          iv%info(bogus)%dxm(:,nlocal)         = platform%loc%dxm
 339          iv%info(bogus)%dy(:,nlocal)          = platform%loc%dy
 340          iv%info(bogus)%dym(:,nlocal)         = platform%loc%dym
 341          iv%info(bogus)%proc_domain(:,nlocal) = platform%loc%proc_domain
 342          ! iv%info(bogus)%proc_domain(:,nlocal) = .true.
 343          ! iv%info(bogus)%proc_domain  = .true.
 344          ! iv%info(bogus)%proc_domain(1,1) = .true.
 345
 346          iv%info(bogus)%obs_global_index(nlocal) = nlocal
 347       enddo
 348    enddo
 349    deallocate(u_lsac)
 350    deallocate(v_lsac)
 351    deallocate(theta_lsac)
 352    deallocate(tk_lsac)
 353    deallocate(p_lsac)
 354    deallocate(pb_lsac)
 355    deallocate(q_lsac)
 356    deallocate(lat_lsac)
 357    deallocate(lon_lsac)
 358    deallocate(phb_lsac)
 359    deallocate(ph_lsac)
 360    deallocate(height_lsac)
 361    deallocate(press)
 362
 363    if (lsac_print_details .and. rootproc) then
 364       close(iunit)
 365       call da_free_unit(iunit)
 366    endif
 367
 368 endif !onlyscan
 369
 370 if (trace_use) call da_trace_exit("da_read_lsac_util")
 371
 372 end subroutine da_read_lsac_wrfinput