Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Enhancements for AHI radiance DA #1774

Merged
merged 16 commits into from
Dec 8, 2022
Merged

Enhancements for AHI radiance DA #1774

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
9f2fa47
Add 2 new namelist entries for AHI assimilation in WRFDA
weather4evr Oct 3, 2022
5512227
Change comment in ahi_info to add JAXA
weather4evr Oct 3, 2022
9867a08
Add new AHI namelist variables so they can be used by subroutines
weather4evr Oct 3, 2022
bd3809f
Remove unneeded comments
weather4evr Oct 4, 2022
2834cba
Expand satinfo structure with new variables for all-sky radiances
weather4evr Oct 4, 2022
bc57756
Allocate/deallocate/initialize stuff, mostly cloud-ice path (CIP), ne…
weather4evr Oct 4, 2022
8af5097
Read extra variables for AHI from radiance_info and allocate new sati…
weather4evr Oct 4, 2022
4a61dcb
Add clear-sky brightness temp calculation for AHI and cloud-ice path …
weather4evr Oct 4, 2022
a15075e
Upgrades to reading AHI data in JAXA format and implementation of bia…
weather4evr Oct 4, 2022
ec9e953
Add symmetric error model for obs error, highlight/fix some questiona…
weather4evr Oct 4, 2022
e5559d7
Output extra fields into ascii file if AHI for background field.
weather4evr Oct 4, 2022
d3ba3bc
Add new variables for AHI
weather4evr Oct 4, 2022
79a4d92
Change AHI-8 channel numbers to be more accurate. Useful for output.
weather4evr Oct 4, 2022
9523624
Merge branch the latest develop code into ahiUpgrades
liujake Nov 23, 2022
158c7a2
Minor changes to da_qc_ahi.inc based on review comments
weather4evr Dec 6, 2022
0f74466
Merge branch 'ahiUpgrades' of https://github.com/weather4evr/WRF into…
weather4evr Dec 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Registry/registry.var
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" "" ""
Expand Down
1 change: 1 addition & 0 deletions var/da/da_define_structures/da_define_structures.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
45 changes: 41 additions & 4 deletions var/da/da_monitor/da_rad_diags.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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) )
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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 )
Expand Down
1 change: 1 addition & 0 deletions var/da/da_radiance/da_allocate_rad_iv.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
13 changes: 13 additions & 0 deletions var/da/da_radiance/da_deallocate_radiance.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down
22 changes: 18 additions & 4 deletions var/da/da_radiance/da_get_innov_vector_crtm.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions var/da/da_radiance/da_initialize_rad_iv.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
51 changes: 41 additions & 10 deletions var/da/da_radiance/da_qc_ahi.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -163,7 +164,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv)

if ( isflg > 0 ) then
do k = 1, nchan
if ( k /= 2 .and. k /= 3 .and. k /= 4 ) then
if ( k /= 2 .and. k /= 3 .and. k /= 4 ) then ! This is questionable...
liujake marked this conversation as resolved.
Show resolved Hide resolved
if (only_sea_rad) then
tb_qc(k) = qc_bad
nrej_land = nrej_land + 1
Expand Down Expand Up @@ -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)
tb_xb_clr(:,n) ! CSS could be a bug; variable not initialized.
liujake marked this conversation as resolved.
Show resolved Hide resolved
end if

if (print_cld_debug) write(stdout,'(A,I8,8F12.4,2x,A)') 'PIXEL_DEBUG5: ', n, &
Expand Down Expand Up @@ -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)
Expand All @@ -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
liujake marked this conversation as resolved.
Show resolved Hide resolved
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
Expand Down Expand Up @@ -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 ! commented out; doesn't seem to be needed
liujake marked this conversation as resolved.
Show resolved Hide resolved
if (.not. crtm_cloud ) then
do k = 1, nchan
if (use_error_factor_rad) then
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion var/da/da_radiance/da_radiance.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion var/da/da_radiance/da_radiance1.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading