var/da/da_radiance/da_get_innov_vector_rttov.inc

   1 subroutine da_get_innov_vector_rttov (it, grid, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculate innovation vector for radiance data.
   5    !
   6    !  METHOD:  d = y - H(x)
   7    !       1. interpolate grid%xb to obs location
   8    !       2. call forward RTM to get simulated bright temperature
   9    !       3. obs BT - simulated BT
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    integer,           intent(in)    :: it       ! External iteration.
  15    type (domain),     intent(in)    :: grid     ! first guess state.
  16    type (y_type),     intent(inout) :: ob       ! Observation structure.
  17    type (iv_type),    intent(inout) :: iv       ! O-B structure.
  18
  19 #if defined(RTTOV)
  20
  21    integer  :: n        ! Loop counter.
  22    integer  :: i, j, k  ! Index dimension.
  23    integer  :: nlevels  ! Number of obs levels.
  24    integer  :: nchanprof, errorstatus
  25
  26    character(len=256) :: atlas_path
  27    real*8                           :: seap, icep, lndp, snop
  28    real, allocatable :: v_p(:,:)
  29
  30    integer :: inst, nchan
  31    real,    allocatable   :: pres(:)
  32
  33    integer :: n1,n2,n1n2
  34
  35 ! FIX?
  36 real,allocatable :: temp(:), temp2(:), temp3(:,:)
  37
  38    type(rttov_emissivity), allocatable :: emissivity(:)
  39
  40    type(con_vars_type),  allocatable :: con_vars(:)
  41    type(aux_vars_type),  allocatable :: aux_vars(:)
  42
  43    type(rttov_chanprof), allocatable :: chanprof(:)
  44    type(rttov_profile), allocatable   :: profiles(:)
  45
  46    ! variables for computing clwp
  47    real, allocatable :: dpf(:,:), clw(:,:), pf(:,:)
  48
  49    real, allocatable :: em_mspps(:)    ! emissivity caluclated using MSPPS algorithm
  50    real              :: ts_mspps       ! surface temperature calcualted using MSPPS algorithm
  51
  52    !For Zhuge and Zou cloud detection
  53    real, allocatable :: geoht_full(:,:,:)
  54    real              :: geoht_pixel(kts:min(kte,kme-1))
  55    real              :: tt_pixel(kts:min(kte,kme-1))
  56    real              :: pp_pixel(kts:min(kte,kme-1))
  57
  58    if (trace_use) call da_trace_entry("da_get_innov_vector_rttov")
  59
  60    !------------------------------------------------------
  61    ! [1.0] calculate the background bright temperature
  62    !-------------------------------------------------------
  63
  64    if ( use_clddet_zz ) then
  65       allocate ( geoht_full(ims:ime,jms:jme,kms:kme-1) )
  66       do k = kms, kme-1
  67          do j = jms, jme
  68             do i = ims, ime
  69                geoht_full(i,j,k) = 0.5 * ( grid%ph_2(i,j,k)   + grid%phb(i,j,k) + &
  70                                            grid%ph_2(i,j,k+1) + grid%phb(i,j,k+1) ) / gravity
  71             end do
  72          end do
  73       end do
  74    end if
  75
  76    do inst = 1, iv%num_inst                 ! loop for sensor
  77       if ( iv%instid(inst)%num_rad < 1 ) cycle
  78       nlevels  = iv%instid(inst)%nlevels
  79       nchan    = iv%instid(inst)%nchan
  80
  81       if (iv%instid(inst)%info%n2 < iv%instid(inst)%info%n1) cycle
  82       n1 = iv%instid(inst)%info%n1
  83       n2 = iv%instid(inst)%info%n2
  84
  85       n1n2=n2-n1+1
  86
  87       allocate (pres(1:nlevels))
  88       allocate (con_vars(n1:n2))
  89       allocate (aux_vars(n1:n2))
  90
  91       pres(1:nlevels) = coefs(inst) % coef % ref_prfl_p(1:nlevels)
  92
  93       allocate(v_p(kms:kme,n1:n2))
  94       v_p(:,:)=0.0
  95       allocate(clw(kms:kme,n1:n2))
  96       allocate(dpf(kms:kme,n1:n2))
  97       allocate(pf(kms:kme+1,n1:n2))
  98
  99       ! horizontal interpolate grid%xb pressure to ob position for every grid%xb levels
 100       do n=n1,n2
 101          do k=kts,kte
 102             ! convert to mb
 103             v_p(k,n) = 0.01*(iv%instid(inst)%info%dym(k,n)*( &
 104                iv%instid(inst)%info%dxm(k,n)*grid%xb%p(iv%instid(inst)%info%i(k,n),  iv%instid(inst)%info%j(k,n),k) + &
 105                iv%instid(inst)%info%dx(k,n) *grid%xb%p(iv%instid(inst)%info%i(k,n)+1,iv%instid(inst)%info%j(k,n),k)) + &
 106                iv%instid(inst)%info%dy(k,n) *( &
 107                iv%instid(inst)%info%dxm(k,n)*grid%xb%p(iv%instid(inst)%info%i(k,n),  iv%instid(inst)%info%j(k,n)+1,k) + &
 108                iv%instid(inst)%info%dx(k,n) *grid%xb%p(iv%instid(inst)%info%i(k,n)+1,iv%instid(inst)%info%j(k,n)+1,k)))
 109          end do
 110       end do
 111
 112       call da_to_zk_new(pres, v_p(:,n1:n2), v_interp_p, n1n2,nlevels,iv%instid(inst)%info%zk(:,n1:n2))
 113
 114       call da_convert_zk (iv%instid(inst)%info)
 115
 116       ! [1.2] Interpolate horizontally to ob:
 117       call da_interp_lin_3d (grid%xb%t, iv%instid(inst)%info, iv%instid(inst)%t (:,n1:n2))
 118       call da_interp_lin_3d (grid%xb%q, iv%instid(inst)%info, iv%instid(inst)%mr(:,n1:n2))
 119
 120       do n= n1,n2
 121          do k=1, nlevels
 122             if (iv%instid(inst)%info%zk(k,n) <= 0.0) then
 123                iv%instid(inst)%t(k,n) = coefs(inst) % coef % ref_prfl_t(k,gas_id_watervapour)  ! outside model level
 124                iv%instid(inst)%mr(k,n) = coefs(inst) % coef % ref_prfl_mr(k,gas_id_watervapour)
 125             else
 126                iv%instid(inst)%mr(k,n) = iv%instid(inst)%mr(k,n) * q2ppmv
 127             end if
 128             if (pres(k) < 100.0) iv%instid(inst)%mr(k,n) = coefs(inst) % coef % ref_prfl_mr(k,gas_id_watervapour)
 129          end do
 130
 131          ! determine surface type of obs location
 132          !-----------------------------------------
 133          call da_detsurtyp( grid%xb%snow, grid%xb%xice, grid%xb%landmask,  &
 134             grid%xb%ivgtyp, grid%xb%isltyp, &
 135             ims, ime, jms, jme, &
 136             iv%instid(inst)%info%i(1,n), iv%instid(inst)%info%j(1,n), &
 137             iv%instid(inst)%info%dx(1,n), iv%instid(inst)%info%dy(1,n), &
 138             iv%instid(inst)%info%dxm(1,n), iv%instid(inst)%info%dym(1,n), &
 139             iv%instid(inst)%isflg(n),iv%instid(inst)%vegtyp(n), iv%instid(inst)%soiltyp(n), &
 140             seap, icep, lndp, snop )
 141
 142          iv%instid(inst)%snow_frac(n) = snop   ! snow coverage fraction 0-1
 143
 144          if ( iv%instid(inst)%isflg(n) == 0 .or. iv%instid(inst)%isflg(n) == 4 ) then       ! sea
 145             iv%instid(inst)%surftype(n) = 1
 146          else if ( iv%instid(inst)%isflg(n) == 1 .or. iv%instid(inst)%isflg(n) == 5 ) then  ! sea-ice with snow
 147             iv%instid(inst)%surftype(n) = 2
 148          else
 149             iv%instid(inst)%surftype(n) = 0
 150          end if
 151
 152          if ( use_clddet_zz ) then
 153             ! Find tropopause temperature for Zhuge and Zou Cloud Detection
 154             do k = kts, min(kte,kme-1)
 155                call da_interp_2d_partial ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
 156                call da_interp_2d_partial ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
 157                call da_interp_2d_partial ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
 158
 159 !               call da_interp_lin_2d ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
 160 !               call da_interp_lin_2d ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
 161 !               call da_interp_lin_2d ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
 162             end do
 163             call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
 164          end if
 165       end do
 166
 167       call da_interp_lin_2d (grid%xb % u10,  iv%instid(inst)%info, 1, iv%instid(inst)%u10(n1:n2))
 168       call da_interp_lin_2d (grid%xb % v10,  iv%instid(inst)%info, 1, iv%instid(inst)%v10(n1:n2))
 169       call da_interp_lin_2d (grid%xb % t2,   iv%instid(inst)%info, 1, iv%instid(inst)%t2m(n1:n2))
 170       call da_interp_lin_2d (grid%xb % q2,   iv%instid(inst)%info, 1, iv%instid(inst)%q2m(n1:n2))
 171       call da_interp_lin_2d (grid%xb % psfc, iv%instid(inst)%info, 1, iv%instid(inst)%ps (n1:n2))
 172       call da_interp_lin_2d (grid%xb % tsk,  iv%instid(inst)%info, 1, iv%instid(inst)%ts (n1:n2))
 173       call da_interp_lin_2d (grid%xb % terr, iv%instid(inst)%info, 1, iv%instid(inst)%elevation(n1:n2))
 174
 175       if ( use_mspps_ts(inst) ) then
 176          ! only for AMSU-A over land
 177          if ( trim(rttov_inst_name(rtminit_sensor(inst))) == 'amsua' ) then
 178             do n = n1, n2
 179                if ( iv%instid(inst)%surftype(n) == 0 ) then
 180                   call da_mspps_ts(ob%instid(inst)%tb(1:nchan,n), nchan,  &
 181                                    iv%instid(inst)%satzen(n), ts_mspps)
 182                   ! ts_mspps is initilaized as negative values in the
 183                   ! da_mspps_ts subroutine.  Apply only valid values here.
 184                   if ( ts_mspps > 0.0 ) then
 185                      iv%instid(inst)%ts(n) = ts_mspps
 186                   end if
 187                end if
 188             end do
 189          end if
 190       end if
 191
 192       ! variables for calculation of cloud affected radiance
 193       !-------------------------------------------------------
 194       do k=kts,kte
 195          call da_interp_lin_2d (grid%xb%t  (:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%tm(k,:))
 196          call da_interp_lin_2d (grid%xb%q  (:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qm(k,:))
 197          call da_interp_lin_2d (grid%xb%qrn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qrn(k,:))
 198          call da_interp_lin_2d (grid%xb%qcw(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qcw(k,:))
 199        ! call da_interp_lin_2d (grid%xb%qci(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qci(k,:))
 200        ! call da_interp_lin_2d (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qsn(k,:))
 201        ! call da_interp_lin_2d (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qgr(k,:))
 202       end do
 203
 204       iv%instid(inst)%pm(:,n1:n2)  = v_p(:,n1:n2)
 205
 206       iv%instid(inst)%ps(n1:n2) = 0.01 * iv%instid(inst)%ps(n1:n2)       ! hPa
 207       iv%instid(inst)%mr2m(n1:n2) = iv%instid(inst)%q2m(n1:n2) * q2ppmv  ! ppmv
 208
 209       ! ADD for computing cloud liquid water path (mm) from guess
 210
 211       pf(kts,n1:n2) = 100.0*iv%instid(inst)%ps(n1:n2)
 212       do k = 2,kte
 213          pf(k,n1:n2) = 50.0*(v_p(k-1,n1:n2)+v_p(k,n1:n2))
 214       end do
 215       pf(kte+1,n1:n2)= 50.0*v_p(kte,n1:n2)
 216
 217       iv%instid(inst)%clwp(n1:n2) = 0.0
 218       do k = kts,kte
 219          dpf(k,n1:n2) = pf(k,n1:n2) - pf(k+1,n1:n2)
 220          clw(k,n1:n2) = iv%instid(inst)%qcw(k,n1:n2)*dpf(k,n1:n2)/gravity
 221          where (v_p(k,n1:n2)<100.0) clw (k,n1:n2) = 0.0
 222          iv%instid(inst)%clwp(n1:n2)  =iv%instid(inst)%clwp(n1:n2) + clw(k,n1:n2)
 223       end do
 224
 225       ! surface emissivity
 226       !-------------------------------------------
 227
 228       nchanprof = nchan*n1n2
 229       allocate(emissivity(nchanprof))
 230       emissivity(:)%emis_in = 0.0
 231
 232       if ( rttov_emis_atlas_ir > 0 .or. rttov_emis_atlas_mw > 0 ) then
 233          ! set up emissivity atlas
 234          atlas_path = 'emis_data/'
 235          write(unit=message(1),fmt='(A,A)')  &
 236             'Setting up emissivity atlas for instrument ', trim(iv%instid(inst)%rttovid_string)
 237          call da_message(message(1:1))
 238          call rttov_setup_emis_atlas(                      &
 239             errorstatus,                                   & ! out
 240             opts(inst),                                    & ! in
 241             grid%start_month,                              & ! in
 242             atlas_type(inst),                              & ! in
 243             atlas,                                         & ! inout
 244             atlas_id(inst),                                & ! in, optional
 245             path = trim(atlas_path),                       & ! in, optional
 246             coefs = coefs(inst))                             ! in
 247
 248          if ( errorstatus /= errorstatus_success ) then
 249             call da_error(__FILE__,__LINE__, &
 250                (/"failure in setting up emissivity atlas"/))
 251          end if
 252
 253          ! Generate the chanprof array
 254          allocate(chanprof(nchanprof))
 255          do n = n1, n2
 256             chanprof((n-n1)*nchan+1:(n-n1+1)*nchan)%prof = n-n1+1
 257             chanprof((n-n1)*nchan+1:(n-n1+1)*nchan)%chan = (/ (j, j=1,nchan) /)
 258          end do
 259
 260          allocate(profiles(n2-n1+1))
 261          do n = n1, n2
 262             ! latitude, longitude, surftype are used for retreiving emis from atlas
 263             ! zenangle is used by MW emmisivity atlas
 264             ! snow_frac is used only by IR emmisivity atlas
 265             profiles(n-n1+1)%latitude           = iv%instid(inst)%info%lat(1,n)
 266             profiles(n-n1+1)%longitude          = iv%instid(inst)%info%lon(1,n)
 267             profiles(n-n1+1)%zenangle           = iv%instid(inst)%satzen(n)
 268             profiles(n-n1+1)%skin%surftype      = iv%instid(inst)%surftype(n)
 269             profiles(n-n1+1)%skin%snow_fraction = iv%instid(inst)%snow_frac(n)
 270          end do
 271
 272          ! Retrieve values from atlas
 273          call rttov_get_emis(                   &
 274             errorstatus,                        & ! out
 275             opts(inst),                         & ! in
 276             chanprof,                           & ! in
 277             profiles,                           & ! in
 278             coefs(inst),                        & ! in
 279             atlas,                              &! in
 280             emissivity=emissivity(:)%emis_in )    ! out
 281          if ( errorstatus /= errorstatus_success ) then
 282             call da_error(__FILE__,__LINE__, &
 283                (/"failure in retrieving emissivity values"/))
 284          end if
 285          deallocate (profiles)
 286          deallocate (chanprof)
 287       end if
 288
 289       if ( use_mspps_emis(inst) ) then
 290          ! Only for AMSU-A over land
 291          if ( trim(rttov_inst_name(rtminit_sensor(inst))) == 'amsua' ) then
 292             allocate(em_mspps(nchan))
 293             do n = n1, n2
 294                if ( iv%instid(inst)%surftype(n) == 0 ) then
 295                   call da_mspps_emis(ob%instid(inst)%tb(1:nchan,n), nchan, em_mspps)
 296                   do k = 1, nchan
 297                      if ( emissivity((n-n1)*nchan+k)%emis_in  < 0.01 ) then
 298                         emissivity((n-n1)*nchan+k)%emis_in  = em_mspps(k)
 299                      end if
 300                   end do
 301                end if
 302             end do
 303             deallocate(em_mspps)
 304          end if
 305       end if
 306
 307       !$OMP PARALLEL DO &
 308       !$OMP PRIVATE ( n, temp, temp2, temp3 )
 309       do n=n1,n2
 310          con_vars(n) % nlevels = nlevels
 311          allocate (con_vars(n) % t(nlevels))
 312          allocate (con_vars(n) % q(nlevels))
 313          if ( use_rttov_kmatrix ) then
 314             allocate (con_vars(n) % t_jac(nchan,nlevels))
 315             allocate (con_vars(n) % q_jac(nchan,nlevels))
 316             allocate (con_vars(n) % ps_jac(nchan))
 317             con_vars(n) % t_jac(:,:) = 0.0
 318             con_vars(n) % q_jac(:,:) = 0.0
 319             con_vars(n) % ps_jac(:)  = 0.0
 320          end if
 321
 322          con_vars(n) % t(1:nlevels) = iv%instid(inst)%t(1:nlevels,n)
 323          con_vars(n) % q(1:nlevels) = iv%instid(inst)%mr(1:nlevels,n)
 324          con_vars(n) % ps           = iv%instid(inst)%ps(n)
 325
 326          aux_vars(n) % t2m      = iv%instid(inst)%t2m(n)
 327          aux_vars(n) % q2m      = iv%instid(inst)%mr2m(n)
 328          aux_vars(n) % u10      = iv%instid(inst)%u10(n)
 329          aux_vars(n) % v10      = iv%instid(inst)%v10(n)
 330          aux_vars(n) % surftype = iv%instid(inst)%surftype(n)
 331          aux_vars(n) % surft    = iv%instid(inst)%ts(n)
 332          aux_vars(n) % satzen   = iv%instid(inst)%satzen(n)
 333          aux_vars(n) % satazi   = iv%instid(inst)%satazi(n)
 334          aux_vars(n) % solzen   = iv%instid(inst)%solzen(n)
 335          aux_vars(n) % solazi   = iv%instid(inst)%solazi(n)
 336          aux_vars(n) % elevation = iv%instid(inst)%elevation(n) !iv%instid(inst)%info%elv(n)
 337          aux_vars(n) % rlat     = iv%instid(inst)%info%lat(1,n)
 338
 339          ! [1.3] Call RTM forward model
 340          ! da_rttov_direct nominally an array version, but doesn't handle arrays
 341          ! of surface flags properly
 342          ! da_rttov_k or da_rttov_direct is used one profile per call
 343          allocate(temp(nchan),temp2(nchan),temp3(nchan,nlevels-1))
 344          if ( use_rttov_kmatrix ) then
 345             call da_rttov_k (inst, nchan, 1, nlevels, &
 346                con_vars(n:n), aux_vars(n:n), &
 347                temp, temp2, temp3, emissivity((n-n1)*nchan+1:(n-n1+1)*nchan))
 348             iv%instid(inst)%emiss(:,n) = emissivity((n-n1)*nchan+1:(n-n1+1)*nchan)%emis_out
 349          else
 350             call da_rttov_direct (inst, nchan, 1, nlevels, &
 351                con_vars(n:n), aux_vars(n:n), &
 352                temp, temp2, temp3, emissivity((n-n1)*nchan+1:(n-n1+1)*nchan))
 353             iv%instid(inst)%emiss(:,n) = emissivity((n-n1)*nchan+1:(n-n1+1)*nchan)%emis_out
 354          end if
 355          iv%instid(inst)%tb_xb(:,n)=temp(:)
 356          iv%instid(inst)%rad_xb(:,n)=temp2(:)
 357          ! Overcast Radiances for AIRS Cloud Detection(MMR)
 358          iv%instid(inst)%rad_ovc(:,1:nlevels-1,n)=temp3(:,:) ! overcast nlayers=nlevels-1
 359          deallocate(temp,temp2,temp3)
 360
 361       end do
 362       !$OMP END PARALLEL DO
 363
 364       if ( use_rttov_kmatrix ) then
 365          do n = n1, n2
 366             do k = 1, nlevels
 367                iv%instid(inst)%t_jacobian(:,k,n) = con_vars(n)%t_jac(:,k)
 368                iv%instid(inst)%q_jacobian(:,k,n) = con_vars(n)%q_jac(:,k)
 369             end do
 370             iv%instid(inst)%ps_jacobian(:,n) = con_vars(n)%ps_jac(:)
 371          end do
 372          do n = n1, n2
 373             deallocate (con_vars(n) % t_jac)
 374             deallocate (con_vars(n) % q_jac)
 375             deallocate (con_vars(n) % ps_jac)
 376          end do
 377       end if
 378
 379       do n=n1,n2
 380          deallocate (con_vars(n) % t)
 381          deallocate (con_vars(n) % q)
 382       end do
 383
 384       !----------------------------------------------------------------
 385       ! [2.0] calculate components of innovation vector:
 386       !----------------------------------------------------------------
 387
 388       do n=n1,n2
 389          do k=1,nchan
 390             if (iv%instid(inst)%tb_inv(k,n) > missing_r) then
 391                iv%instid(inst)%tb_inv(k,n) = ob%instid(inst)%tb(k,n) - iv%instid(inst)%tb_xb(k,n)
 392             else
 393                iv%instid(inst)%tb_inv(k,n)   = missing_r
 394                iv%instid(inst)%tb_qc(k,n)    = qc_bad
 395             end if
 396          end do
 397       end do
 398
 399       deallocate (v_p)
 400       deallocate (clw)
 401       deallocate (dpf)
 402       deallocate (pf)
 403       deallocate (pres)
 404       deallocate (con_vars)
 405       deallocate (aux_vars)
 406       deallocate (emissivity)
 407
 408       if ( rttov_emis_atlas_ir > 0 .or. rttov_emis_atlas_mw > 0 ) then
 409          call rttov_deallocate_emis_atlas(atlas)
 410       end if
 411
 412    end do        ! end loop for sensor
 413
 414    if ( use_clddet_zz ) deallocate ( geoht_full )
 415
 416    if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
 417 #else
 418     call da_error(__FILE__,__LINE__, &
 419        (/"Must compile with $RTTOV option for radiances"/))
 420 #endif
 421
 422 end subroutine da_get_innov_vector_rttov
 423