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