Enhancements for AHI radiance DA (#1774)

TYPE: Enhancement 

KEYWORDS: WRFDA, AHI, Himawari-8

SOURCE: Craig Schwartz (NCAR/MMM)

DESCRIPTION OF CHANGES:
This PR makes several enhancements for assimilating Himawari-8 radiance data, including
(1) Introduction of an all-sky obs error model (Harnisch et al., 2016) for all-sky AHI DA;
(2) Optional read and use of AHI level-2 product (e.g., cloud mask); 
(3) More efficient read of a sub-area of the full disk data;
(4) Allow the use of offline statistics of constant bias correction values;
(5) More diagnostic output (peak of weighting functions, cloud ice water path, cloud flag) in omb/oma files.

LIST OF MODIFIED FILES:
var/da/da_radiance/module_radiance.f90 Registry/registry.var
var/da/da_define_structures/da_define_structures.f90
var/da/da_monitor/da_rad_diags.f90
var/da/da_radiance/da_allocate_rad_iv.inc
var/da/da_radiance/da_deallocate_radiance.inc
var/da/da_radiance/da_get_innov_vector_crtm.inc
var/da/da_radiance/da_initialize_rad_iv.inc
var/da/da_radiance/da_qc_ahi.inc
var/da/da_radiance/da_radiance.f90
var/da/da_radiance/da_radiance1.f90
var/da/da_radiance/da_radiance_init.inc
var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc
var/da/da_radiance/da_write_iv_rad_ascii.inc
var/da/da_radiance/module_radiance.f90
var/run/ahi_info
var/run/radiance_info/himawari-8-ahi.info

TESTING:
(1) WRFDA regression test passed;
(2) AHI all-sky DA runs also Ok.

Release Note: Enhancements for AHI radiance DA, including all-sky observation error model, Level-2 AHI product read, and more diagnostic output.
Xu, D. M., Z. Q. Liu, S. Y. Fan, M. Chen, and F. F. Shen, 2021: Assimilating all-sky infrared radiances from Himawari-8 using the 3DVar method for the prediction of a severe storm over North China. Adv. Atmos. Sci., 38(4), 661-676.

Registry/registry.var

@@ -467,6 +467,8 @@ rconfig   character crtm_irwater_coef       namelist,wrfvar14  1  "Nalli.IRwater
 rconfig   character crtm_mwwater_coef       namelist,wrfvar14  1  "FASTEM5.MWwater.EmisCoeff.bin" - "crtm_mwwater_coef" "" ""
 rconfig   character crtm_irland_coef        namelist,wrfvar14  1  "USGS.IRland.EmisCoeff.bin"     - "crtm_irland_coef"  "" ""
 rconfig   character crtm_visland_coef       namelist,wrfvar14  1  "USGS.VISland.EmisCoeff.bin"    - "crtm_visland_coef" "" ""
+rconfig   logical   ahi_use_symm_obs_err    namelist,wrfvar14  1  .false.  - "ahi_use_symm_obs_err"         ""  ""
+rconfig   logical   ahi_apply_clrsky_bias   namelist,wrfvar14  1  .false.  - "ahi_apply_clrsky_bias"         ""  ""
 rconfig   integer   num_pseudo              namelist,wrfvar15  1  0        - "num_pseudo"              ""  ""   
 rconfig   real      pseudo_x                namelist,wrfvar15  1  1.0      - "pseudo_x"                ""  ""   
 rconfig   real      pseudo_y                namelist,wrfvar15  1  1.0      - "pseudo_y"                ""  ""   

var/da/da_define_structures/da_define_structures.f90

@@ -636,6 +636,7 @@ module da_define_structures
       real,    pointer     :: vegtyp(:)
       real,    pointer     :: vegfra(:)
       real,    pointer     :: clwp(:) ! model/guess clwp
+      real,    pointer     :: cip(:)  ! model/guess cloud-ice path
       real,    pointer     :: clw(:)  ! currently AMSR2 only
       real,    pointer     :: ps_jacobian(:,:) ! only RTTOV
       real,    pointer     :: ts_jacobian(:,:) ! only over water CRTM

var/da/da_monitor/da_rad_diags.f90

@@ -60,8 +60,12 @@ program da_rad_diags
    real*4,  dimension(:), allocatable     :: ret_clw
    real*4,  dimension(:), allocatable     :: satzen, satazi, solzen, solazi, t2m, mr2m, u10, v10, ps, ts
    real*4,  dimension(:), allocatable     :: smois, tslb, snowh, vegfra, clwp, cloud_frac
+   real*4,  dimension(:), allocatable     :: cip ! cloud-ice path
+   integer, dimension(:), allocatable     :: cloudflag ! cloudflag from L2 AHI
    integer, dimension(:,:), allocatable   :: tb_qc
    real*4,  dimension(:,:), allocatable   :: tb_obs, tb_bak, tb_inv, tb_oma, tb_err, ems, ems_jac
+   real*4,  dimension(:,:), allocatable   :: tb_bak_clr ! clear-sky brightness temp
+   real*4,  dimension(:,:), allocatable   :: weightfunc_peak ! peak of weighting function
    real*4,  dimension(:,:), allocatable   :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water
    real*4,  dimension(:,:), allocatable   :: prf_ice, prf_rain, prf_snow, prf_grau, prf_hail
    real*4,  dimension(:,:), allocatable   :: prf_water_reff, prf_ice_reff, prf_rain_reff
@@ -248,9 +252,13 @@ program da_rad_diags
             allocate (       vegfra(1:total_npixel) )
             allocate (         elev(1:total_npixel) )
             allocate (         clwp(1:total_npixel) ) !model/guess clwp
+            allocate (         cip(1:total_npixel) ) ! model/guess cloud-ice path
+            allocate (   cloudflag(1:total_npixel) ) ! AHI Level-2 cloud flag
             allocate (   cloud_frac(1:total_npixel) )
             allocate ( tb_obs(1:nchan,1:total_npixel) )
             allocate ( tb_bak(1:nchan,1:total_npixel) )
+            allocate ( tb_bak_clr(1:nchan,1:total_npixel) )
+            allocate ( weightfunc_peak(1:nchan,1:total_npixel) )
             allocate ( tb_inv(1:nchan,1:total_npixel) )
             allocate ( tb_oma(1:nchan,1:total_npixel) )
             allocate ( tb_err(1:nchan,1:total_npixel) )
@@ -320,11 +328,14 @@ program da_rad_diags
             ! initialize
             tb_obs = missing_r
             tb_bak = missing_r
+            tb_bak_clr = missing_r
+            weightfunc_peak = missing_r
             tb_inv = missing_r
             tb_oma = missing_r
             tb_err = missing_r
             ncname = 'diags_'//trim(instid(iinst))//"_"//datestr1(itime)//'.nc'
-            ios = NF_CREATE(trim(ncname), NF_CLOBBER, ncid)  ! NF_CLOBBER specifies the default behavior of
+            ios = NF_CREATE(trim(ncname), NF_NETCDF4, ncid)  ! Change to output netcdf4 files
+           !ios = NF_CREATE(trim(ncname), NF_CLOBBER, ncid)  ! NF_CLOBBER specifies the default behavior of
                                                              ! overwritting any existing dataset with the 
                                                              ! same file name
             if ( ios /= 0 ) then
@@ -349,15 +360,23 @@ program da_rad_diags
                lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel), solzen(ipixel), &
                solazi(ipixel), ret_clw(ipixel)
 
-            read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2,f15.5)',iostat=ios)  &
+            read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2,f15.5,f15.5,i7)',iostat=ios)  &
                t2m(ipixel), mr2m(ipixel), u10(ipixel), v10(ipixel), ps(ipixel), ts(ipixel), &
                smois(ipixel), tslb(ipixel), snowh(ipixel), isflg(ipixel), soiltyp(ipixel),  &
-               vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel), cloud_frac(ipixel)
+               vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel), cloud_frac(ipixel), cip(ipixel), cloudflag(ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! OBS
             read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_obs(:,ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! BAK
             read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_bak(:,ipixel)
-            read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! IVBC
+            read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! IVBC or BAK_clr or WEIGHTFUNC_PEAK
+            if ( buf(1:7) == "BAK_clr" ) then
+               read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_bak_clr(:,ipixel)
+               read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf        ! IVBC or WEIGHTFUNC_PEAK
+            end if
+            if ( buf(1:15) == "WEIGHTFUNC_PEAK" ) then
+               read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) weightfunc_peak(:,ipixel)
+               read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf        ! IVBC
+            end if
             read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_inv(:,ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! OMA or EMS
             if ( buf(1:3) == "OMA" ) then
@@ -490,6 +509,10 @@ program da_rad_diags
       ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
       ios = NF_DEF_VAR(ncid, 'tb_bak', NF_FLOAT, 2, ishape(1:2), varid)
       ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
+      ios = NF_DEF_VAR(ncid, 'tb_bak_clr', NF_FLOAT, 2, ishape(1:2), varid)
+      ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
+      ios = NF_DEF_VAR(ncid, 'weightfunc_peak', NF_FLOAT, 2, ishape(1:2), varid)
+      ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
       ios = NF_DEF_VAR(ncid, 'tb_inv', NF_FLOAT, 2, ishape(1:2), varid)
       ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
       ios = NF_DEF_VAR(ncid, 'tb_oma', NF_FLOAT, 2, ishape(1:2), varid)
@@ -601,6 +624,8 @@ program da_rad_diags
       ios = NF_DEF_VAR(ncid, 'vegfra',       NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'elev',         NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'clwp',         NF_FLOAT, 1, ishape(1), varid)
+      ios = NF_DEF_VAR(ncid, 'cip',          NF_FLOAT, 1, ishape(1), varid)
+      ios = NF_DEF_VAR(ncid, 'cloudflag',    NF_INT,   1, ishape(1), varid)
       if ( amsr2 ) then
          ios = NF_DEF_VAR(ncid, 'ret_clw',   NF_FLOAT, 1, ishape(1), varid)
       end if
@@ -626,6 +651,10 @@ program da_rad_diags
       ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_obs)
       ios = NF_INQ_VARID (ncid, 'tb_bak', varid)
       ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_bak)
+      ios = NF_INQ_VARID (ncid, 'tb_bak_clr', varid)
+      ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_bak_clr)
+      ios = NF_INQ_VARID (ncid, 'weightfunc_peak', varid)
+      ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), weightfunc_peak)
       ios = NF_INQ_VARID (ncid, 'tb_inv', varid)
       ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_inv)
       ios = NF_INQ_VARID (ncid, 'tb_oma', varid)
@@ -809,6 +838,10 @@ program da_rad_diags
       ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), elev)
       ios = NF_INQ_VARID (ncid, 'clwp', varid)
       ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), clwp)
+      ios = NF_INQ_VARID (ncid, 'cip', varid)
+      ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), cip)
+      ios = NF_INQ_VARID (ncid, 'cloudflag', varid)
+      ios = NF_PUT_VARA_INT(ncid,  varid, istart(2), icount(2), cloudflag)
       if ( amsr2 ) then
          ios = NF_INQ_VARID (ncid, 'ret_clw', varid)
          ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), ret_clw)
@@ -848,10 +881,14 @@ program da_rad_diags
       deallocate ( vegfra )
       deallocate ( elev )
       deallocate ( clwp )
+      deallocate ( cip )
+      deallocate ( cloudflag )
       deallocate ( cloud_frac )
       deallocate ( ret_clw )
       deallocate ( tb_obs )
       deallocate ( tb_bak )
+      deallocate ( tb_bak_clr )
+      deallocate ( weightfunc_peak )
       deallocate ( tb_inv )
       deallocate ( tb_oma )
       deallocate ( ems )

var/da/da_radiance/da_allocate_rad_iv.inc

@@ -88,6 +88,7 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
       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))
+         allocate (iv%instid(i)%cip(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))

var/da/da_radiance/da_deallocate_radiance.inc

@@ -30,6 +30,18 @@
          deallocate (satinfo(i) % bcoef0)
          deallocate (satinfo(i) % error_std)
 
+	 ! Deallocate extra variables for AHI
+         if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then
+	    deallocate ( satinfo(i) % BTLim)
+	    deallocate ( satinfo(i) % ca1)
+	    deallocate ( satinfo(i) % ca2)
+	    deallocate ( satinfo(i) % clearSkyBias)
+         endif
+
+         if (use_error_factor_rad) then
+            deallocate (satinfo(i) % error_factor)
+         endif
+
          deallocate (ob%instid(i) % ichan)
          deallocate (iv%instid(i) % ichan)
 
@@ -111,6 +123,7 @@
             deallocate (iv%instid(i)%tb_xb)
             if ( crtm_cloud ) then
                deallocate (iv%instid(i)%tb_xb_clr)
+               deallocate (iv%instid(i)%cip)
             end if
             deallocate (iv%instid(i)%tb_qc)
             deallocate (iv%instid(i)%tb_inv)

var/da/da_radiance/da_get_innov_vector_crtm.inc

@@ -41,6 +41,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
    integer :: model_isflg
 
    real    :: model_qcw(kms:kme)
+   real    :: model_qci(kms:kme)
    real    :: model_rho(kms:kme)
    real, allocatable :: model_snow(:)  ! snow water equivalent, different from model_snowh,
                                        ! used in calculating reff_water
@@ -55,9 +56,9 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
    integer :: inst, nchanl, n1,n2
    integer :: ipred, npred, gammapred
 
-   ! variables for computing clwp
-   real    :: clw(kms:kme), dpf(kms:kme)
-   real    :: clwp
+   ! variables for computing clwp (clw/clwp) and cloud ice (ci/cip)
+   real    :: clw(kms:kme), dpf(kms:kme), ci(kms:kme)
+   real    :: clwp, cip
 
    real    :: total_od
 
@@ -241,7 +242,8 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
 
       calc_tb_clr = .false.
       if ( crtm_cloud .and. &
-           trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' ) then
+           ( trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' .or. &
+             trim( crtm_sensor_name(rtminit_sensor(inst))) == 'ahi') ) 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.
@@ -389,6 +391,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
             if (crtm_cloud) then
 
                call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,qci)
+               call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,model_qci(kte-k+1:kte-k+1))
 
                call da_interp_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n,qrn)
 
@@ -511,6 +514,16 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
             clwp  = clwp + clw(k)
          end do
 
+         ! computing ice path from guess
+         if (crtm_cloud) then
+            cip = 0.0
+            do k = kts,kte ! from top to bottom
+               dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
+               ci(k) = model_qci(k)*dpf(k)/gravity ! kg/m2 or mm
+               cip  = cip + ci(k)
+            end do
+	 endif
+
          ! CRTM GeometryInfo Structure
          GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
          GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
@@ -853,6 +866,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
          iv%instid(inst)%vegtyp(n)    = model_ivgtyp
          iv%instid(inst)%vegfra(n)    = model_vegfra(n)
          iv%instid(inst)%clwp(n)      = clwp
+         if (crtm_cloud) iv%instid(inst)%cip(n) = cip
          iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage
          iv%instid(inst)%land_coverage(n)  = Surface(1)%land_coverage
          iv%instid(inst)%ice_coverage(n)   = Surface(1)%ice_coverage                              

var/da/da_radiance/da_initialize_rad_iv.inc

@@ -84,6 +84,7 @@ subroutine da_initialize_rad_iv (i, n, iv, p)
    iv%instid(i)%tb_xb(:,n)      = 0.0
    if ( crtm_cloud ) then
       iv%instid(i)%tb_xb_clr(:,n) = 0.0
+      iv%instid(i)%cip(n)         = 0.0
    end if
    iv%instid(i)%tb_inv(:,n)     = p%tb_inv(:)
    iv%instid(i)%tb_qc(:,n)      = 0

var/da/da_radiance/da_qc_ahi.inc

@@ -58,6 +58,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
                                       num_clddet_tests)
    integer            :: superob_center
    integer            :: isuper, jsuper
+   real               :: cm, co, ca ! variables for all-sky obs error
 
    real, pointer      :: tb_ob(:,:), tb_xb(:,:), tb_inv(:,:), tb_xb_clr(:,:), ca_mean(:,:)
    integer            :: tb_qc(nchan), tb_qc_clddet(nchan)
@@ -215,7 +216,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
             tb_ob(:,n)
          if (crtm_cloud ) then
             if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG4: ', n, &
-               tb_xb_clr(:,n)
+               iv%instid(i)%tb_xb_clr(:,n)
          end if
 
          if (print_cld_debug) write(stdout,'(A,I8,8F12.4,2x,A)') 'PIXEL_DEBUG5: ', n, &
@@ -410,7 +411,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
                              qual_clddet )
 
             if (reject_clddet) then
-               tb_qc_clddet = qc_bad
+               tb_qc_clddet = qc_bad ! CSS do we want to set it bad for a given pixel within a superob?
                if (iv%instid(i)%info%proc_domain(1,n)) then
                   nrej_clddet(:,itest) = nrej_clddet(:,itest) + 1
         !          write(stdout,"(A,F14.6,A,I4,2D12.4)") trim(crit_names_clddet(itest)), crit_clddet, " isflg", isflg, iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n)
@@ -422,12 +423,10 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
          end do ! isuper
          end do ! jsuper
 
-         if ( iv%instid(i)%superob_width > 1 ) then
-            iv%instid(i)%cloud_frac(n) = &
-               real( count(sum(clddet_tests,3) > 0),8) / real( iv%instid(i)%superob_width**2,8)
-         end if
+         iv%instid(i)%cloud_frac(n) = &
+            real( count(sum(clddet_tests,3) > 0),8) / real( iv%instid(i)%superob_width**2,8)
 
-         if (.not. crtm_cloud ) tb_qc = tb_qc_clddet
+         if (.not. crtm_cloud ) tb_qc = tb_qc_clddet ! CSS logic here isn't quite right if superobbing
 
 !JJGDEBUG
          if (print_cld_debug) write(stdout,'(A,I8,*(2x,I1:))') 'PIXEL_DEBUG6: ', n, clddet_tests
@@ -476,7 +475,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
 
       ! ---------------------------calculate and save ca_mean for crtm_cloud and crtm_clr
       ! 5.0 assigning obs errors
-	     tb_xb_clr => iv%instid(i)%tb_xb_clr  	 
+	   ! tb_xb_clr => iv%instid(i)%tb_xb_clr ! currently not used
       if (.not. crtm_cloud ) then
          do k = 1, nchan
             if (use_error_factor_rad) then
@@ -486,6 +485,35 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
                iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
             end if
          end do ! nchan
+
+      else ! Added this else block...until now obs error = 500 if crtm_cloud = T...not good!
+
+         if ( ahi_use_symm_obs_err ) then
+	    do k = 1, nchan
+	      ! Okamato et al. (2014)
+	      !cm = iv%instid(i)%tb_xb(k,n) - iv%instid(i)%tb_xb_clr(k,n)
+	      !co = ob%instid(i)%tb(k,n)    - iv%instid(i)%tb_xb_clr(k,n)
+
+	      ! Harnisch et al. (2016)
+	       cm = max(0.0, satinfo(i)%BTLim(k) - iv%instid(i)%tb_xb(k,n))
+	       co = max(0.0, satinfo(i)%BTLim(k) - ob%instid(i)%tb(k,n))
+
+	       ! Symmetric cloud amount
+	       ca = 0.5*( abs(cm) + abs(co) )
+
+	       ! Figure out observation error as a function of ca
+	       if (ca.lt.satinfo(i)%ca1(k)) then
+		  iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)
+	       else if (ca.ge.satinfo(i)%ca1(k) .and. ca.lt.satinfo(i)%ca2(k)) then
+		  iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)+ &
+		     (ca-satinfo(i)%ca1(k))*(satinfo(i)%error_cld(k)-satinfo(i)%error_std(k))/(satinfo(i)%ca2(k)-satinfo(i)%ca1(k))
+	       else
+		  iv%instid(i)%tb_error(k,n)= satinfo(i)%error_cld(k)
+	       end if
+	    end do ! nchan
+	 else
+	    iv%instid(i)%tb_error(:,n)= 500.0 ! this is the default
+	 endif
       end if
 
       !  6.0 gross 8k check -clr,sdobs-clr
@@ -528,7 +556,10 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)
       end do	 ! nchan
 
       !final QC decsion	 	
-      if (crtm_cloud ) tb_qc(2:4) = qc_good  ! no qc  for crtm_cloud
+      ! CSS comment out below. Dangerous, especially if satinfo(i)%iuse(k) = -1
+      ! CSS also okay to fail 3std check if using symmetric error model for obs errors
+      !if (crtm_cloud ) tb_qc(2:4) = qc_good  ! no qc  for crtm_cloud
+
       iv%instid(i)%tb_qc(:,n) = tb_qc
       do k = 1, nchan     		 
            if (iv%instid(i)%tb_qc(k,n) == qc_bad) then

var/da/da_radiance/da_radiance.f90

@@ -58,7 +58,7 @@ module da_radiance
       use_rad,crtm_cloud, DT_cloud_model, global, use_varbc, freeze_varbc, &
       airs_warmest_fov, time_slots, interp_option, ids, ide, jds, jde, &
       ips, ipe, jps, jpe, simulated_rad_ngrid, obs_qc_pointer, use_blacklist_rad, use_satcv, &
-      use_goesimgobs, pi, earth_radius, satellite_height,use_clddet_zz, ahi_superob_halfwidth
+      use_goesimgobs, pi, earth_radius, satellite_height,use_clddet_zz, ahi_superob_halfwidth, ahi_apply_clrsky_bias
 
 #ifdef CRTM
    use da_crtm, only : da_crtm_init, da_get_innov_vector_crtm

var/da/da_radiance/da_radiance1.f90

@@ -22,7 +22,7 @@ module da_radiance1
       use_pseudo_rad, pi, t_triple, crtm_cloud, DT_cloud_model,write_jacobian, &
       use_crtm_kmatrix,use_clddet, use_satcv, cv_size_domain, &
       cv_size_domain_js, calc_weightfunc, deg_to_rad, rad_to_deg,use_clddet_zz, &
-      ahi_superob_halfwidth
+      ahi_superob_halfwidth, ahi_use_symm_obs_err
    use da_define_structures, only : info_type,model_loc_type,maxmin_type, &
       iv_type, y_type, jo_type,bad_data_type,bad_data_type,number_type, &
       be_type, clddet_geoir_type, superob_type