New DA capability of QC and obs error modeling for AMSR-2 cloudy…

var/da/da_define_structures/da_define_structures.f90

@@ -497,6 +497,7 @@ module da_define_structures
       integer, pointer     :: tb_qc(:,:)
       real,    pointer     :: tb_error(:,:)
       real,    pointer     :: tb_xb(:,:) 
+      real,    pointer     :: tb_xb_clr(:,:) 
       real,    pointer     :: tb_sens(:,:)
       real,    pointer     :: tb_imp(:,:)
       real,    pointer     :: rad_xb(:,:)

var/da/da_radiance/da_allocate_rad_iv.inc

@@ -80,6 +80,9 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
       end if
       allocate (iv%instid(i)%ps(iv%instid(i)%num_rad))
       allocate (iv%instid(i)%tb_xb(nchan,iv%instid(i)%num_rad))
+      if ( crtm_cloud ) then
+         allocate (iv%instid(i)%tb_xb_clr(nchan,iv%instid(i)%num_rad))
+      end if
       allocate (iv%instid(i)%tb_qc(nchan,iv%instid(i)%num_rad))
       allocate (iv%instid(i)%tb_inv(nchan,iv%instid(i)%num_rad))
       allocate (iv%instid(i)%tb_error(nchan,iv%instid(i)%num_rad))

var/da/da_radiance/da_deallocate_radiance.inc

@@ -21,6 +21,7 @@
          deallocate (satinfo(i) % ichan)
          deallocate (satinfo(i) % iuse)
          deallocate (satinfo(i) % error)
+         deallocate (satinfo(i) % error_cld)
          deallocate (satinfo(i) % polar)
 
          deallocate (satinfo(i) % scanbias) 
@@ -102,6 +103,9 @@
             end if
             deallocate (iv%instid(i)%ps)
             deallocate (iv%instid(i)%tb_xb)
+            if ( crtm_cloud ) then
+               deallocate (iv%instid(i)%tb_xb_clr)
+            end if
             deallocate (iv%instid(i)%tb_qc)
             deallocate (iv%instid(i)%tb_inv)
             deallocate (iv%instid(i)%tb_error)

var/da/da_radiance/da_get_innov_vector_crtm.inc

@@ -83,7 +83,9 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
 
 !! for crtm cloud
    real :: qcw(1),qrn(1), qci(1),qsn(1),qgr(1)
-
+   logical :: calc_tb_clr
+   type (CRTM_RTSolution_type), allocatable :: RTSolution_clr(:,:)
+
    ! Initializations for AIRS (MMR) Cloud Detection
    integer           :: ilev, jlev, nclouds
    real, allocatable :: hessian(:,:)
@@ -222,7 +224,28 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
          call da_error(__FILE__,__LINE__, &
             (/"Error in allocating CRTM RTSolution Stucture"/))
       end if
-
+
+      calc_tb_clr = .false.
+      if ( crtm_cloud .and. &
+           trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' ) then
+         !Tb_clear_sky is only needed for symmetric obs error model
+         !symmetric obs error model only implemented for amsr2 for now
+         calc_tb_clr = .true.
+      end if
+
+      if ( calc_tb_clr ) then
+         allocate (RTSolution_clr(ChannelInfo(inst)%n_Channels,1), STAT = Allocate_Status )
+         if ( Allocate_Status /= 0 ) then
+            call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution_clr"/))
+         end if
+         call CRTM_RTSolution_Create( RTSolution_clr , &              !Output
+                                      Atmosphere(1)%n_layers )        !Input
+         if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution_clr))) ) then
+            call da_error(__FILE__,__LINE__, &
+               (/"Error in allocating CRTM RTSolution_clr Stucture"/))
+         end if
+      end if !calc_tb_clr
+
       if (use_crtm_kmatrix) then
          allocate (RTSolution_K(ChannelInfo(inst)%n_Channels,1), &
                    Atmosphere_K(ChannelInfo(inst)%n_Channels,1), &
@@ -820,6 +843,22 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
          iv%instid(inst)%land_coverage(n)  = Surface(1)%land_coverage
          iv%instid(inst)%ice_coverage(n)   = Surface(1)%ice_coverage                              
          iv%instid(inst)%snow_coverage(n)  = Surface(1)%snow_coverage
+
+         if ( calc_tb_clr ) then
+            atmosphere(1)%n_clouds = 0
+            call da_crtm_direct(1, nchanl, 1, Atmosphere,   &
+               Surface,                                     &
+               GeometryInfo,                                &
+               ChannelInfo(inst:inst),                      &
+               RTSolution_clr,                              &
+               Options)
+            do k = 1, nchanl
+               iv%instid(inst)%tb_xb_clr(k,n) = RTSolution_clr(k,1)%Brightness_Temperature
+            end do
+            ! resetting n_clouds
+            atmosphere(1)%n_clouds = n_clouds
+         end if !calc_tb_clr
+
       end do pixel_loop      !  end loop for pixels
 
       if ( use_mspps_emis(inst) ) deallocate(em_mspps)
@@ -841,6 +880,14 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
          call da_error(__FILE__,__LINE__, &
             (/"Error in deallocating RTSolution"/))
 
+      if ( calc_tb_clr ) then
+         call CRTM_RTSolution_Destroy(RTSolution_clr)
+         deallocate( RTSolution_clr, STAT = Allocate_Status )
+         if (Allocate_Status /= 0) &
+            call da_error(__FILE__,__LINE__, &
+               (/"Error in deallocating RTSolution_clr"/))
+      end if !calc_tb_clr
+
       if (use_crtm_kmatrix) then
          call CRTM_RTSolution_Destroy(RTSolution_K)
          deallocate( RTSolution_K, STAT = Allocate_Status )

var/da/da_radiance/da_initialize_rad_iv.inc

@@ -74,6 +74,9 @@ subroutine da_initialize_rad_iv (i, n, iv, p)
    end if
    iv%instid(i)%ps(n)           = 0.0
    iv%instid(i)%tb_xb(:,n)      = 0.0
+   if ( crtm_cloud ) then
+      iv%instid(i)%tb_xb_clr(:,n) = 0.0
+   end if
    iv%instid(i)%tb_inv(:,n)     = p%tb_inv(:)
    iv%instid(i)%tb_qc(:,n)      = 0
    iv%instid(i)%tb_error(:,n)   = 500.0

var/da/da_radiance/da_qc_amsr2.inc

@@ -25,6 +25,7 @@ subroutine da_qc_amsr2 (it, i, nchan, ob, iv)
                 nrej_land, nrej_glint
 
    character(len=30)  :: filename
+   real      :: c37_mean
 
    if (trace_use) call da_trace_entry("da_qc_amsr2")
 
@@ -44,6 +45,13 @@ subroutine da_qc_amsr2 (it, i, nchan, ob, iv)
       if (iv%instid(i)%info%proc_domain(1,n)) &
             num_proc_domain = num_proc_domain + 1
 
+      if ( crtm_cloud ) then
+         ! calculate c37_mean
+         c37_mean = 1.0-(ob%instid(i)%tb(11,n)-ob%instid(i)%tb(12,n)+  &
+                    iv%instid(i)%tb_xb(11,n)-iv%instid(i)%tb_xb(12,n))/ &
+                    (2.0*(iv%instid(i)%tb_xb_clr(11,n)-iv%instid(i)%tb_xb_clr(12,n)))
+      end if
+
       !  0.0  initialise QC by flags assuming good obs
       !-----------------------------------------------------------------
       iv%instid(i)%tb_qc(:,n) = qc_good
@@ -96,38 +104,61 @@ subroutine da_qc_amsr2 (it, i, nchan, ob, iv)
          end do
       end if
 
+      ! assigning obs errors
+      if (.not. crtm_cloud ) then
+         do k = 1, nchan
+            if (use_error_factor_rad) then
+               iv%instid(i)%tb_error(k,n) = &
+                   satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
+            else
+               iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
+            end if
+         end do ! nchan
+
+      else !crtm_cloud
+         ! symmetric error model, Geer and Bauer (2011)
+         do k = 1, nchan
+            if (c37_mean.lt.0.05) then
+               iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)
+            else if (c37_mean.ge.0.05.and.c37_mean.lt.0.5) then
+               iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)+ &
+                  (c37_mean-0.05)*(satinfo(i)%error_cld(k)-satinfo(i)%error_std(k))/(0.5-0.05)
+            else
+               iv%instid(i)%tb_error(k,n)= satinfo(i)%error_cld(k)
+            end if
+         end do ! nchan
+
+      end if
+
       !  5.0 check innovation
       !-----------------------------------------------------------------
-      do k = 1, nchan
-
+      if (.not. crtm_cloud ) then
          ! absolute departure check
-         if ( k <= 7 .or. k == 11 .or. k == 12) then
-            if (abs(iv%instid(i)%tb_inv(k,n)) > 6.0) then
-               iv%instid(i)%tb_qc(k,n)  = qc_bad
-               if (iv%instid(i)%info%proc_domain(1,n)) &
-                     nrej_omb_abs(k) = nrej_omb_abs(k) + 1
-            end if
-         else if ( k == 8 .or. k == 9 ) then
-            if (abs(iv%instid(i)%tb_inv(k,n)) > 8.0) then
-               iv%instid(i)%tb_qc(k,n)  = qc_bad
-               if (iv%instid(i)%info%proc_domain(1,n)) &
-                     nrej_omb_abs(k) = nrej_omb_abs(k) + 1
-            end if
-         else
-            if (abs(iv%instid(i)%tb_inv(k,n)) > 10.0) then
-               iv%instid(i)%tb_qc(k,n)  = qc_bad
-               if (iv%instid(i)%info%proc_domain(1,n)) &
-                     nrej_omb_abs(k) = nrej_omb_abs(k) + 1
+         do k = 1, nchan
+            if ( k <= 7 .or. k == 11 .or. k == 12) then
+               if (abs(iv%instid(i)%tb_inv(k,n)) > 6.0) then
+                  iv%instid(i)%tb_qc(k,n)  = qc_bad
+                  if (iv%instid(i)%info%proc_domain(1,n)) &
+                        nrej_omb_abs(k) = nrej_omb_abs(k) + 1
+               end if
+            else if ( k == 8 .or. k == 9 ) then
+               if (abs(iv%instid(i)%tb_inv(k,n)) > 8.0) then
+                  iv%instid(i)%tb_qc(k,n)  = qc_bad
+                  if (iv%instid(i)%info%proc_domain(1,n)) &
+                        nrej_omb_abs(k) = nrej_omb_abs(k) + 1
+               end if
+            else
+               if (abs(iv%instid(i)%tb_inv(k,n)) > 10.0) then
+                  iv%instid(i)%tb_qc(k,n)  = qc_bad
+                  if (iv%instid(i)%info%proc_domain(1,n)) &
+                        nrej_omb_abs(k) = nrej_omb_abs(k) + 1
+               end if
             end if
-         end if
+         end do ! nchan
+      end if
 
+      do k = 1, nchan
          ! relative departure check
-         if (use_error_factor_rad) then
-            iv%instid(i)%tb_error(k,n) = &
-                satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
-         else
-            iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
-         end if
          if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
             iv%instid(i)%tb_qc(k,n)  = qc_bad
             if (iv%instid(i)%info%proc_domain(1,n)) &

var/da/da_radiance/da_radiance_init.inc

@@ -33,6 +33,8 @@ subroutine da_radiance_init(iv,ob)
  integer             :: idum, wmo_sensor_id, sensor_type, iost
  integer             :: iunit
  character(len=filename_len)  :: filename
+ character(len=20)   :: cdum
+ real                :: error_cld
 
 ! local variables for tuning error factor
 !----------------------------------------
@@ -161,18 +163,40 @@ subroutine da_radiance_init(iv,ob)
    allocate ( satinfo(n) % ichan(nchanl(n)) )
    allocate ( satinfo(n) % iuse (nchanl(n)) )
    allocate ( satinfo(n) % error(nchanl(n)) )
+   allocate ( satinfo(n) % error_cld(nchanl(n)) )
    allocate ( satinfo(n) % polar(nchanl(n)) )
 
+   satinfo(n) % error_cld(:) = 500.0 !initialize
+
    read(iunit,*)
    do j = 1, nchanl(n)
-     read(iunit,'(1x,5i5,2e18.10)')    &
+      read(iunit,'(1x,5i5,2e18.10,a20)') &
+                     wmo_sensor_id, &
+               satinfo(n)%ichan(j), &
+                       sensor_type, &
+               satinfo(n)%iuse(j) , &
+                              idum, &
+               satinfo(n)%error(j), &
+               satinfo(n)%polar(j), &
+               cdum
+      !in the current radiance info files, the last column
+      !can be either sensor_id_string or blank
+      if ( len_trim(cdum) > 0 .and. index(cdum,'-') == 0 ) then
+         ! read the line again to get error_cld when it is available
+         backspace(iunit)
+         read(iunit,'(1x,5i5,2e18.10,f10.5)')  &
                      wmo_sensor_id, &
                satinfo(n)%ichan(j), &
                        sensor_type, &
                satinfo(n)%iuse(j) , &
                               idum, &
                satinfo(n)%error(j), &
-               satinfo(n)%polar(j)
+               satinfo(n)%polar(j), &
+               error_cld
+         if ( error_cld > 0.0 ) then
+            satinfo(n)%error_cld(j) = error_cld
+         end if
+      end if
      iv%instid(n)%ichan(j) = satinfo(n)%ichan(j)
      ob%instid(n)%ichan(j) = satinfo(n)%ichan(j)
    end do

var/da/da_radiance/module_radiance.f90

@@ -143,9 +143,10 @@ module module_radiance
 
    type satinfo_type
       integer, pointer   :: ichan(:)      ! channel index
-      integer, pointer   :: iuse (:)      ! usage flag (-1: not use) from GSI info file
-      real   , pointer   :: error(:)      ! error Standard Deviation from GSI info file
-      real   , pointer   :: polar(:)      ! polarisation (0:ver; 1:hori) from GSI info file
+      integer, pointer   :: iuse (:)      ! usage flag (-1: not use) from radiance info file
+      real   , pointer   :: error(:)      ! error Standard Deviation from radiance info file
+      real   , pointer   :: error_cld(:)  ! error Standard Deviation for cloudy radiance from radiance info file
+      real   , pointer   :: polar(:)      ! polarisation (0:ver; 1:hori) from radiance info file
       real   , pointer   :: error_factor(:) ! error tuning factor ! from error tuning file
      ! new air mass bias correction coefs.
       real   , pointer   :: scanbias(:,:) ! scan bias without latitude band variation

var/run/radiance_info/gcom-w-1-amsr2.info

@@ -1,15 +1,15 @@
 sensor channel IR/MW use idum  varch  polarisation(0:vertical;1:horizontal)
-   478    1    1   -1    0  0.7260000000E+00  0.0000000000E+00  gcom-w1-amsr2
-   478    2    1   -1    0  0.9560000000E+00  1.0000000000E+00
-   478    3    1   -1    0  0.7750000000E+00  0.0000000000E+00
-   478    4    1   -1    0  0.9910000000E+00  1.0000000000E+00
-   478    5    1    1    0  0.8660000000E+00  0.0000000000E+00
-   478    6    1    1    0  1.1290000000E+00  1.0000000000E+00
-   478    7    1    1    0  1.2270000000E+00  0.0000000000E+00
-   478    8    1    1    0  1.7470000000E+00  1.0000000000E+00
-   478    9    1    1    0  1.6000000000E+00  0.0000000000E+00
-   478   10    1    1    0  2.6790000000E+00  1.0000000000E+00
-   478   11    1    1    0  1.1790000000E+00  0.0000000000E+00
-   478   12    1    1    0  2.2680000000E+00  1.0000000000E+00
-   478   13    1   -1    0  2.1310000000E+00  0.0000000000E+00
-   478   14    1   -1    0  4.0750000000E+00  1.0000000000E+00
+   478    1    1   -1    0  0.7260000000E+00  0.0000000000E+00  10.14561
+   478    2    1   -1    0  0.9560000000E+00  1.0000000000E+00  18.24548
+   478    3    1   -1    0  0.7750000000E+00  0.0000000000E+00  11.14696
+   478    4    1   -1    0  0.9910000000E+00  1.0000000000E+00  20.18668
+   478    5    1    1    0  0.8660000000E+00  0.0000000000E+00  21.93555
+   478    6    1    1    0  1.1290000000E+00  1.0000000000E+00  40.92418
+   478    7    1    1    0  1.2270000000E+00  0.0000000000E+00  28.30175
+   478    8    1    1    0  1.7470000000E+00  1.0000000000E+00  57.58830
+   478    9    1    1    0  1.6000000000E+00  0.0000000000E+00  12.69287
+   478   10    1    1    0  2.6790000000E+00  1.0000000000E+00  27.33099
+   478   11    1    1    0  1.1790000000E+00  0.0000000000E+00  23.24269
+   478   12    1    1    0  2.2680000000E+00  1.0000000000E+00  53.35099
+   478   13    1   -1    0  2.1310000000E+00  0.0000000000E+00  36.07700
+   478   14    1   -1    0  4.0750000000E+00  1.0000000000E+00  33.61592