Skip to content

Commit

Permalink
Add the capability for assimilating GOES-ABI radiance data (#1983)
Browse files Browse the repository at this point in the history 
TYPE: new feature

KEYWORDS: ABI, cloud detection, all-sky obs error model

SOURCE: JJ Guerrette (NCAR/MMM, now at tomorrow.io), Deqin Li (Liaoning Meteorological Bureau of CMA), Jake Liu (NCAR/MMM)

DESCRIPTION OF CHANGES:
This PR addes the assimilation of GOES-16/17 ABI's 3 water vapor channels' radiance data. This includes reading of ABI's full-disk, CONUS, and meso1&2 data files, superobbing and thinning of ABI data, IR-based cloud detection scheme as part of quality control, and all-sky obs error model. Cloud detection scheme should be the same as for AHI in principle, but the actual code implementation is not the same. No attempt made to make ABI's cloud detection code consistent with AHI's when bringing the code originally developed back in 2019-2020 into the latest develop branch. Some technical information is provided below for the use of this new capability.

1. Read ABI files:
Raw netcdf ABI data files (one file for one channel) need to be listed in 4 list files:
'file_list_GOES-16-ABI_C' for CONUS scan files
'file_list_GOES-16-ABI_F'  for full-disk scan files (e.g., OR_ABI-L1b-RadF-M6C08_G16_s20191211200263_e20191211209571_c20191211210021.nc)
'file_list_GOES-16-ABI_M1' for meso1 scan files
'file_list_GOES-16-ABI_M2' for meso2 scan files

ABI data reading code will automatically determine which file(s) to read in matching ABI file's time and analysis time.

2. Cloud detection scheme needs to read in a terrain file OR_ABI-TERR_G16.nc or OR_ABI-TERR_G17.nc for GOES-16 ABI or GOES-17 ABI.

3. Related namelist settings:
```
&wrfvar4
 use_goesabiobs  = true, ! read goes-16 and goes-17 ABI data
/
```
```
&wrfvar14
 rtminit_nsensor= 1,
 rtminit_platform=  4, ! goes
 rtminit_satid=  16,      ! goes-16
 rtminit_sensor=  44,  ! abi
 thinning= true,
 thinning_mesh= 30.0,
 qc_rad=true,
 write_iv_rad_ascii=true,
 write_oa_rad_ascii=true,
 rtm_option= 2,
 crtm_cloud= false,
 only_sea_rad=false,
 use_varbc=true,
 varbc_nobsmin=500,
 crtm_irland_coef= "IGBP.IRland.EmisCoeff.bin",
 use_clddet_zz=true, ! IR-based cloud detection
 abi_superob_halfwidth=3, ! this will do supperobbing with 7x7 pixels
 /
```

See also AHI DA related PRs:
#1139
#1173
#1774

LIST OF MODIFIED FILES: 41
M       Registry/registry.var
M       var/build/depend.txt
M       var/da/da_define_structures/da_define_structures.f90
M       var/da/da_monitor/da_rad_diags.f90
M       var/da/da_radiance/da_allocate_rad_iv.inc
M       var/da/da_radiance/da_deallocate_radiance.inc
M       var/da/da_radiance/da_get_innov_vector_crtm.inc
M       var/da/da_radiance/da_get_innov_vector_rttov.inc
A       var/da/da_radiance/da_get_sat_angles.inc
A       var/da/da_radiance/da_get_sat_angles_1d.inc
A       var/da/da_radiance/da_get_solar_angles.inc
A       var/da/da_radiance/da_get_solar_angles_1d.inc
M       var/da/da_radiance/da_initialize_rad_iv.inc
A       var/da/da_radiance/da_qc_goesabi.inc
M       var/da/da_radiance/da_qc_rad.inc
M       var/da/da_radiance/da_radiance.f90
M       var/da/da_radiance/da_radiance1.f90
M       var/da/da_radiance/da_radiance_init.inc
A       var/da/da_radiance/da_read_obs_ncgoesabi.inc
M       var/da/da_radiance/da_rttov.f90
M       var/da/da_radiance/da_setup_radiance_structures.inc
M       var/da/da_radiance/da_write_iv_rad_ascii.inc
M       var/da/da_radiance/da_write_oa_rad_ascii.inc
M       var/da/da_radiance/module_radiance.f90
M       var/da/da_setup_structures/da_setup_obs_structures.inc
M       var/da/da_setup_structures/da_setup_structures.f90
A       var/da/da_tools/da_llxy_1d.inc
A       var/da/da_tools/da_llxy_default_1d.inc
A       var/da/da_tools/da_llxy_global_1d.inc
A       var/da/da_tools/da_llxy_kma_global_1d.inc
A       var/da/da_tools/da_llxy_latlon_1d.inc
A       var/da/da_tools/da_llxy_lc_1d.inc
A       var/da/da_tools/da_llxy_merc_1d.inc
A       var/da/da_tools/da_llxy_ps_1d.inc
A       var/da/da_tools/da_llxy_rotated_latlon_1d.inc
A       var/da/da_tools/da_llxy_wrf_1d.inc
A       var/da/da_tools/da_togrid_1d.inc
M       var/da/da_tools/da_tools.f90
M       var/run/VARBC.in
A       var/run/radiance_info/goes-16-abi.info
A       var/run/radiance_info/goes-17-abi.info

TESTS CONDUCTED: 
1. WRFDA regression test passed on Derecho.
2. Clear-sky ABI DA is tested with a full-disk data file. 
3. the Jenkins tests all passing.

RELEASE NOTE: Add the capability for assimilating GOES-ABI radiance data
  • Loading branch information
@liujake
liujake authored Jan 17, 2024
1 parent 715df80 commit 32de5e4
Show file tree
Hide file tree
Showing 41 changed files with 5,026 additions and 49 deletions.
3 changes: 3 additions & 0 deletions Registry/registry.var
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,7 @@ rconfig logical use_amsr2obs namelist,wrfvar4 1 .false. - "use
rconfig logical use_ahiobs namelist,wrfvar4 1 .false. - "use_ahiobs" "" ""
rconfig logical use_gmiobs namelist,wrfvar4 1 .false. - "use_gmiobs" "" ""
rconfig logical use_goesimgobs namelist,wrfvar4 1 .false. - "use_goesimgobs" "" ""
rconfig logical use_goesabiobs namelist,wrfvar4 1 .false. - "use_goesabiobs" "" ""
rconfig logical use_kma1dvar namelist,wrfvar4 1 .false. - "use_kma1dvar" "" ""
rconfig logical use_filtered_rad namelist,wrfvar4 1 .false. - "use_filtered_rad" "" ""
rconfig logical use_obs_errfac namelist,wrfvar4 1 .false. - "use_obs_errfac" "" ""
Expand Down Expand Up @@ -468,6 +469,7 @@ rconfig integer varbc_nobsmin namelist,wrfvar14 1 10 - "va
rconfig integer use_clddet namelist,wrfvar14 1 2 - "use_clddet" "0: off, 1: mmr, 2: pf, 3: ecmwf" ""
rconfig logical use_clddet_zz namelist,wrfvar14 1 .false. - "use_clddet_zz" "cloud detection scheme from Zhuge X. and Zou X. JAMC, 2016." ""
rconfig integer ahi_superob_halfwidth namelist,wrfvar14 1 0 - "ahi_superob_halfwidth" "" ""
rconfig integer abi_superob_halfwidth namelist,wrfvar14 1 0 - "abi_superob_halfwidth" "" ""
rconfig logical airs_warmest_fov namelist,wrfvar14 1 .false. - "airs_warmest_fov" "" ""
rconfig logical use_satcv namelist,wrfvar14 2 .false. - "use_satcv" "" ""
rconfig logical use_blacklist_rad namelist,wrfvar14 1 .true. - "use_blacklist_rad" "" ""
Expand All @@ -477,6 +479,7 @@ 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 abi_use_symm_obs_err namelist,wrfvar14 1 .false. - "abi_use_symm_obs_err" "" ""
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" "" ""
Expand Down
12 changes: 6 additions & 6 deletions var/build/depend.txt
View file

Large diffs are not rendered by default.

20 changes: 12 additions & 8 deletions var/da/da_define_structures/da_define_structures.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -574,10 +574,12 @@ module da_define_structures
real, pointer :: vtox(:,:)
end type varbc_type
type clddet_geoir_type
real :: RTCT, RFMFT, TEMPIR, terr_hgt
real :: tb_stddev_10, tb_stddev_13,tb_stddev_14
real :: CIRH2O
!real, allocatable :: CIRH2O(:,:,:)
real :: RTCT, RFMFT, TEMPIR, terr_hgt ! for both ABI and AHI
real :: tb_stddev_10, tb_stddev_13,tb_stddev_14 ! only for AHI
real :: CIRH2O ! for both ABI and AHI
real, allocatable :: CIRH2O_abi(:,:,:) ! only for ABI
real, allocatable :: tb_stddev_3x3(:) ! only for ABI
integer :: RFMFT_ij(2) ! only for ABI
end type clddet_geoir_type
type superob_type
real, allocatable :: tb_obs(:,:)
Expand Down Expand Up @@ -618,6 +620,8 @@ module da_define_structures
integer, pointer :: cloud_flag(:,:)
integer, pointer :: cloudflag(:)
integer, pointer :: rain_flag(:)
real, pointer :: cloud_mod(:,:) ! only for ABI
real, pointer :: cloud_obs(:,:) ! only for ABI
real, allocatable :: cloud_frac(:)
real, pointer :: satzen(:)
real, pointer :: satazi(:)
Expand All @@ -632,10 +636,10 @@ module da_define_structures
real, pointer :: lod(:,:,:) ! layer_optical_depth
real, pointer :: trans(:,:,:) ! layer transmittance
real, pointer :: der_trans(:,:,:) ! d(transmittance)/dp
real, pointer :: kmin_t(:)
real, pointer :: kmax_p(:)
real, pointer :: sensitivity_ratio(:,:,:)
real, pointer :: p_chan_level(:,:)
real, pointer :: kmin_t(:)
real, pointer :: kmax_p(:)
real, pointer :: sensitivity_ratio(:,:,:)
real, pointer :: p_chan_level(:,:)
real, pointer :: qrn(:,:)
real, pointer :: qcw(:,:)
real, pointer :: qci(:,:)
Expand Down
48 changes: 44 additions & 4 deletions var/da/da_monitor/da_rad_diags.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ program da_rad_diags
integer :: ncid, dimid, varid
integer, dimension(3) :: ishape, istart, icount
!
logical :: amsr2
logical :: amsr2, abi
logical :: isfile, prf_found, jac_found
integer, parameter :: datelen1 = 10
integer, parameter :: datelen2 = 19
Expand All @@ -62,9 +62,9 @@ program da_rad_diags
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
integer, dimension(:,:), allocatable :: tb_qc, cloud_flag
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 :: cloud_mod, cloud_obs, 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
Expand Down Expand Up @@ -139,6 +139,7 @@ program da_rad_diags
write(0,*) trim(instid(iinst))

amsr2 = index(instid(iinst),'amsr2') > 0
abi = index(instid(iinst),'abi') > 0

nerr = 0
total_npixel = 0
Expand Down Expand Up @@ -263,6 +264,12 @@ program da_rad_diags
allocate ( tb_oma(1:nchan,1:total_npixel) )
allocate ( tb_err(1:nchan,1:total_npixel) )
allocate ( tb_qc(1:nchan,1:total_npixel) )
if ( abi ) then
allocate ( cloud_mod(1:nchan,1:total_npixel) )
allocate ( cloud_obs(1:nchan,1:total_npixel) )
allocate ( cloud_flag(1:nchan,1:total_npixel))
cloud_flag = 0
end if
allocate ( ems(1:nchan,1:total_npixel) )
if ( jac_found ) then
allocate ( ems_jac(1:nchan,1:total_npixel) )
Expand Down Expand Up @@ -333,6 +340,11 @@ program da_rad_diags
tb_inv = missing_r
tb_oma = missing_r
tb_err = missing_r
if ( abi ) then
cloud_mod = missing_r
cloud_obs = missing_r
end if

ncname = 'diags_'//trim(instid(iinst))//"_"//datestr1(itime)//'.nc'
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
Expand Down Expand Up @@ -392,7 +404,15 @@ program da_rad_diags
read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_err(:,ipixel)
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! QC
read(unit=iunit(iproc),fmt='(10i11)',iostat=ios ) tb_qc(:,ipixel)
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf
if ( abi .and. buf(1:4) == "CMOD" ) then ! read cloud_mod, cloud_obs, cloud_flag for abi
read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) cloud_mod(:,ipixel)
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! CMOD
read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) cloud_obs(:,ipixel)
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! COBS
read(unit=iunit(iproc),fmt='(10i11)',iostat=ios ) cloud_flag(:,ipixel)
read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! cloud_flag
end if
if ( buf(1:4) == "INFO" ) then
backspace(iunit(iproc))
cycle npixel_loop
Expand Down Expand Up @@ -523,6 +543,13 @@ program da_rad_diags
end if
ios = NF_DEF_VAR(ncid, 'tb_err', NF_FLOAT, 2, ishape(1:2), varid)
ios = NF_DEF_VAR(ncid, 'tb_qc', NF_INT, 2, ishape(1:2), varid)
if ( abi ) then
ios = NF_DEF_VAR(ncid, 'cloud_mod', 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, 'cloud_obs', 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, 'cloud_flag', NF_INT, 2, ishape(1:2), varid)
end if
!
! define 2-D array with dimensions nlev * total_npixel
!
Expand Down Expand Up @@ -669,6 +696,14 @@ program da_rad_diags
ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_err)
ios = NF_INQ_VARID (ncid, 'tb_qc', varid)
ios = NF_PUT_VARA_INT(ncid, varid, istart(1:2), icount(1:2), tb_qc)
if ( abi ) then
ios = NF_INQ_VARID (ncid, 'cloud_mod', varid)
ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), cloud_mod)
ios = NF_INQ_VARID (ncid, 'cloud_obs', varid)
ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), cloud_obs)
ios = NF_INQ_VARID (ncid, 'cloud_flag', varid)
ios = NF_PUT_VARA_INT(ncid, varid, istart(1:2), icount(1:2), cloud_flag)
end if
!
! output 2-D array with dimensions nlev * total_npixel
!
Expand Down Expand Up @@ -890,6 +925,11 @@ program da_rad_diags
deallocate ( tb_bak_clr )
deallocate ( weightfunc_peak )
deallocate ( tb_inv )
if ( abi ) then
deallocate ( cloud_mod )
deallocate ( cloud_obs )
deallocate ( cloud_flag )
end if
deallocate ( tb_oma )
deallocate ( ems )
if ( jac_found ) deallocate ( ems_jac )
Expand Down
18 changes: 16 additions & 2 deletions var/da/da_radiance/da_allocate_rad_iv.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,10 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
end if
if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then
allocate (iv%instid(i)%cloudflag(iv%instid(i)%num_rad))
end if
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
allocate (iv%instid(i)%cloud_mod(nchan,iv%instid(i)%num_rad))
allocate (iv%instid(i)%cloud_obs(nchan,iv%instid(i)%num_rad))
end if
if ( index(iv%instid(i)%rttovid_string, 'gmi') > 0 ) then
allocate (iv%instid(i)%clw(iv%instid(i)%num_rad))
Expand Down Expand Up @@ -112,16 +116,26 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
allocate (iv%instid(i)%gamma_jacobian(nchan,iv%instid(i)%num_rad))
allocate (iv%instid(i)%cloud_frac(iv%instid(i)%num_rad))
if ( use_clddet_zz ) then
iv%instid(i)%superob_width = 2*ahi_superob_halfwidth+1
! here we assume AHI and ABI (they cover different regions) are not used simultaneously
if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) &
iv%instid(i)%superob_width = 2*ahi_superob_halfwidth+1
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) &
iv%instid(i)%superob_width = 2*abi_superob_halfwidth+1

allocate (iv%instid(i)%superob(iv%instid(i)%superob_width, &
iv%instid(i)%superob_width))
do iy = 1, iv%instid(i)%superob_width
do ix = 1, iv%instid(i)%superob_width
allocate (iv%instid(i)%superob(ix,iy)%cld_qc(iv%instid(i)%num_rad))
allocate (iv%instid(i)%superob(ix,iy)%tb_obs(nchan,iv%instid(i)%num_rad))
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
do n = 1, iv%instid(i)%num_rad
allocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3(nchan))
end do
end if
end do
end do
end if
end if
if ( use_rttov_kmatrix .or. use_crtm_kmatrix ) then
allocate(iv%instid(i)%ts_jacobian(nchan,iv%instid(i)%num_rad))
allocate(iv%instid(i)%ps_jacobian(nchan,iv%instid(i)%num_rad))
Expand Down
23 changes: 21 additions & 2 deletions var/da/da_radiance/da_deallocate_radiance.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,13 @@
deallocate ( satinfo(i) % clearSkyBias)
endif

! Deallocate extra variables for ABI
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
deallocate (satinfo(i) % error_cld_y)
deallocate (satinfo(i) % error_cld_x)
endif


if (use_error_factor_rad) then
deallocate (satinfo(i) % error_factor)
endif
Expand Down Expand Up @@ -115,6 +122,10 @@
end if
if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then
deallocate (iv%instid(i)%cloudflag)
end if
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
deallocate (iv%instid(i)%cloud_mod)
deallocate (iv%instid(i)%cloud_obs)
end if
if ( index(iv%instid(i)%rttovid_string,'gmi') > 0 ) then
deallocate (iv%instid(i)%clw)
Expand Down Expand Up @@ -149,8 +160,16 @@
if ( use_clddet_zz ) then
do iy = 1, iv%instid(i)%superob_width
do ix = 1, iv%instid(i)%superob_width
deallocate (iv%instid(i)%superob(ix,iy)%cld_qc)
deallocate (iv%instid(i)%superob(ix,iy)%tb_obs)
if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
do n = 1,iv%instid(i)%num_rad
if ( allocated (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3) ) &
deallocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3)
if ( allocated (iv%instid(i)%superob(ix,iy)%cld_qc(n)%CIRH2O_abi) ) &
deallocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%CIRH2O_abi)
end do
end if
deallocate (iv%instid(i)%superob(ix,iy)%cld_qc)
deallocate (iv%instid(i)%superob(ix,iy)%tb_obs)
end do
end do
deallocate (iv%instid(i)%superob)
Expand Down
14 changes: 11 additions & 3 deletions var/da/da_radiance/da_get_innov_vector_crtm.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
real, allocatable :: hessian(:,:)
real*8, allocatable :: eignvec(:,:), eignval(:)
real :: rad_clr, rad_ovc_ilev, rad_ovc_jlev

integer :: Band_Size(5), Bands(AIRS_Max_Channels,5)
!For Zhuge and Zou cloud detection
real, allocatable :: geoht_full(:,:,:)
Expand Down Expand Up @@ -243,9 +243,10 @@ 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' .or. &
trim( crtm_sensor_name(rtminit_sensor(inst))) == 'abi' .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
!symmetric obs error model only implemented for amsr2 & abi/ahi for now
calc_tb_clr = .true.
end if

Expand Down Expand Up @@ -443,7 +444,6 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
end if


call da_interp_2d_partial (grid%xb%u10, iv%instid(inst)%info, 1, n, n, model_u10(n:n))
call da_interp_2d_partial (grid%xb%v10, iv%instid(inst)%info, 1, n, n, model_v10(n:n))
call da_interp_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
Expand Down Expand Up @@ -476,6 +476,14 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
cycle pixel_loop
end if
end do
!if ( all(ob%instid(inst)%tb(1:nchanl,n) < 0.) ) then
! write(message(1),'(a,2i5.0,a)') ' Skipping the pixel at loc ', i, j, &
! ' where all observed BTs are < 0'
! call da_warning(__FILE__,__LINE__,message(1:1))
! iv%instid(inst)%tb_inv(:,n) = missing_r
! iv%instid(inst)%info%proc_domain(:,n) = .false.
! cycle pixel_loop
!end if

! convert cloud content unit from kg/kg to kg/m^2
if (crtm_cloud) then
Expand Down
34 changes: 33 additions & 1 deletion var/da/da_radiance/da_get_innov_vector_rttov.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,12 +49,30 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
real, allocatable :: em_mspps(:) ! emissivity caluclated using MSPPS algorithm
real :: ts_mspps ! surface temperature calcualted using MSPPS algorithm

!For Zhuge and Zou cloud detection
real, allocatable :: geoht_full(:,:,:)
real :: geoht_pixel(kts:min(kte,kme-1))
real :: tt_pixel(kts:min(kte,kme-1))
real :: pp_pixel(kts:min(kte,kme-1))

if (trace_use) call da_trace_entry("da_get_innov_vector_rttov")

!------------------------------------------------------
! [1.0] calculate the background bright temperature
!-------------------------------------------------------

if ( use_clddet_zz ) then
allocate ( geoht_full(ims:ime,jms:jme,kms:kme-1) )
do k = kms, kme-1
do j = jms, jme
do i = ims, ime
geoht_full(i,j,k) = 0.5 * ( grid%ph_2(i,j,k) + grid%phb(i,j,k) + &
grid%ph_2(i,j,k+1) + grid%phb(i,j,k+1) ) / gravity
end do
end do
end do
end if

do inst = 1, iv%num_inst ! loop for sensor
if ( iv%instid(inst)%num_rad < 1 ) cycle
nlevels = iv%instid(inst)%nlevels
Expand Down Expand Up @@ -99,7 +117,6 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
call da_interp_lin_3d (grid%xb%t, iv%instid(inst)%info, iv%instid(inst)%t (:,n1:n2))
call da_interp_lin_3d (grid%xb%q, iv%instid(inst)%info, iv%instid(inst)%mr(:,n1:n2))


do n= n1,n2
do k=1, nlevels
if (iv%instid(inst)%info%zk(k,n) <= 0.0) then
Expand Down Expand Up @@ -132,6 +149,19 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
iv%instid(inst)%surftype(n) = 0
end if

if ( use_clddet_zz ) then
! Find tropopause temperature for Zhuge and Zou Cloud Detection
do k = kts, min(kte,kme-1)
call da_interp_2d_partial ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
call da_interp_2d_partial ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
call da_interp_2d_partial ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )

! call da_interp_lin_2d ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
! call da_interp_lin_2d ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
! call da_interp_lin_2d ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
end do
call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
end if
end do

call da_interp_lin_2d (grid%xb % u10, iv%instid(inst)%info, 1, iv%instid(inst)%u10(n1:n2))
Expand Down Expand Up @@ -381,6 +411,8 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)

end do ! end loop for sensor

if ( use_clddet_zz ) deallocate ( geoht_full )

if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
#else
call da_error(__FILE__,__LINE__, &
Expand Down
Loading

0 comments on commit 32de5e4

Please sign in to comment.