Skip to content

A new cloud detection scheme based on Particle Filter #887

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

Merged
liujake merged 7 commits into from
May 3, 2019
Merged

A new cloud detection scheme based on Particle Filter #887

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4f9df08
use_clddet_pf
Apr 22, 2019
850c77f
fix
Apr 22, 2019
3b2a213
Delete Makefile
xdmjolly Apr 23, 2019
fa2f62a
change use_clddet flag
Apr 30, 2019
31816c7
Merge branch 'pf-cloud-detect' of https://github.com/xdmjolly/WRF int…
Apr 30, 2019
a1f5a43
Update registry.var
liujake May 3, 2019
1a5922f
Update da_cloud_detect.inc
liujake May 3, 2019
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
3 changes: 1 addition & 2 deletions Registry/registry.var
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -439,8 +439,7 @@ rconfig logical freeze_varbc namelist,wrfvar14 1 .false. - "fr
rconfig real varbc_factor namelist,wrfvar14 1 1.0 - "varbc_factor" "" ""
rconfig integer varbc_nbgerr namelist,wrfvar14 1 5000 - "varbc_nbgerr" "" ""
rconfig integer varbc_nobsmin namelist,wrfvar14 1 10 - "varbc_nobsmin" "" ""
rconfig logical use_clddet_mmr namelist,wrfvar14 1 .false. - "use_clddet_mmr" "" ""
rconfig logical use_clddet_ecmwf namelist,wrfvar14 1 .false. - "use_clddet_ecmwf" "" ""
rconfig integer use_clddet namelist,wrfvar14 1 2 - "use_clddet" "0: off, 1: mmr, 2: pf, 3: ecmwf" ""
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 Down
2 changes: 1 addition & 1 deletion var/build/depend.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ da_qscat.o : da_qscat.f90 da_calculate_grady_qscat.inc da_transform_xtoy_qscat_a
da_rad_diags.o : da_rad_diags.f90
da_radar.o : da_radar.f90 da_write_oa_radar_ascii.inc da_max_error_qc_radar.inc da_calculate_grady_radar.inc da_radial_velocity_adj.inc da_radial_velocity_lin.inc da_radial_velocity.inc da_radar_rf.inc da_get_innov_vector_radar.inc da_check_max_iv_radar.inc da_transform_xtoy_radar_adj.inc da_transform_xtoy_radar.inc da_print_stats_radar.inc da_oi_stats_radar.inc da_residual_radar.inc da_jo_and_grady_radar.inc da_ao_stats_radar.inc da_tools_serial.o da_reporting.o da_tracing.o da_tools.o da_statistics.o da_par_util1.o da_par_util.o da_interpolation.o da_define_structures.o da_control.o module_domain.o
da_radiance.o : da_radiance.f90 da_blacklist_rad.inc da_read_pseudo_rad.inc da_get_innov_vector_radiance.inc da_radiance_init.inc da_setup_radiance_structures.inc da_sort_rad.inc da_read_kma1dvar.inc da_initialize_rad_iv.inc da_allocate_rad_iv.inc da_read_obs_bufrssmis.inc da_read_obs_bufrairs.inc da_read_obs_bufriasi.inc da_read_obs_bufrseviri.inc da_read_obs_bufrtovs.inc da_write_filtered_rad.inc da_read_simulated_rad.inc da_read_filtered_rad.inc da_calculate_grady_rad.inc gsi_thinning.o da_wrf_interfaces.o da_varbc.o da_tracing.o da_tools.o da_statistics.o da_rttov.o da_reporting.o da_radiance1.o da_physics.o da_par_util.o da_par_util1.o da_tools_serial.o da_interpolation.o da_define_structures.o da_crtm.o da_control.o module_radiance.o module_domain.o amsr2time_.c da_read_obs_hdf5amsr2.inc da_deallocate_radiance.inc da_read_obs_ncgoesimg.inc da_get_satzen.inc da_read_obs_hdf5ahi.inc da_read_obs_netcdf4ahi_jaxa.inc da_read_obs_netcdf4ahi_geocat.inc
da_radiance1.o : da_radiance1.f90 da_mspps_ts.inc da_mspps_emis.inc da_setup_satcv.inc da_qc_rad.inc da_print_stats_rad.inc da_oi_stats_rad.inc da_ao_stats_rad.inc da_cld_eff_radius.inc da_detsurtyp.inc da_write_oa_rad_ascii.inc da_write_iv_rad_ascii.inc da_qc_mhs.inc da_qc_ssmis.inc da_qc_hirs.inc da_qc_amsub.inc da_qc_amsua.inc da_qc_airs.inc da_cloud_detect_airs.inc da_cloud_sim.inc da_qc_seviri.inc da_qc_iasi.inc da_cloud_detect_iasi.inc da_qc_crtm.inc da_predictor_crtm.inc da_predictor_rttov.inc da_write_biasprep.inc da_biasprep.inc da_read_biascoef.inc da_biascorr.inc da_residual_rad.inc da_jo_and_grady_rad.inc gsi_constants.o da_tracing.o da_tools_serial.o da_tools.o da_statistics.o da_reporting.o da_par_util1.o da_par_util.o module_dm.o da_define_structures.o da_control.o module_radiance.o da_qc_amsr2.inc da_qc_goesimg.inc da_qc_ahi.inc
da_radiance1.o : da_radiance1.f90 da_mspps_ts.inc da_mspps_emis.inc da_setup_satcv.inc da_qc_rad.inc da_print_stats_rad.inc da_oi_stats_rad.inc da_ao_stats_rad.inc da_cld_eff_radius.inc da_detsurtyp.inc da_write_oa_rad_ascii.inc da_write_iv_rad_ascii.inc da_qc_mhs.inc da_qc_ssmis.inc da_qc_hirs.inc da_qc_amsub.inc da_qc_amsua.inc da_qc_airs.inc da_cloud_detect.inc da_cloud_sim.inc da_qc_seviri.inc da_qc_iasi.inc da_qc_crtm.inc da_predictor_crtm.inc da_predictor_rttov.inc da_write_biasprep.inc da_biasprep.inc da_read_biascoef.inc da_biascorr.inc da_residual_rad.inc da_jo_and_grady_rad.inc gsi_constants.o da_tracing.o da_tools_serial.o da_tools.o da_statistics.o da_reporting.o da_par_util1.o da_par_util.o module_dm.o da_define_structures.o da_control.o module_radiance.o da_qc_amsr2.inc da_qc_goesimg.inc da_qc_ahi.inc
da_rain.o : da_rain.f90 da_calculate_grady_rain.inc da_get_innov_vector_rain.inc da_get_hr_rain.inc da_check_max_iv_rain.inc da_transform_xtoy_rain_adj.inc da_transform_xtoy_rain.inc da_print_stats_rain.inc da_oi_stats_rain.inc da_residual_rain.inc da_jo_and_grady_rain.inc da_ao_stats_rain.inc da_tracing.o da_tools.o da_statistics.o da_par_util.o da_par_util1.o da_interpolation.o da_define_structures.o da_control.o module_comm_dm.o module_dm.o module_domain.o
da_recursive_filter.o : da_recursive_filter.f90 da_apply_rf_adj.inc da_apply_rf.inc da_apply_rf_1v_adj.inc da_apply_rf_1v.inc da_transform_through_rf_adj.inc da_transform_through_rf.inc da_recursive_filter_1d_adj.inc da_recursive_filter_1d.inc da_calculate_rf_factors.inc da_transform_through_rf_dual_res.inc da_transform_through_rf_adj_dual_res.inc da_perform_2drf.inc da_rf_cv3.o da_rfz_cv3.o da_tracing.o da_par_util.o da_define_structures.o da_control.o module_domain.o
da_reporting.o : da_reporting.f90 da_message2.inc da_message.inc da_warning.inc da_error.inc da_control.o
Expand Down
2 changes: 1 addition & 1 deletion 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 @@ -142,7 +142,7 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
allocate(iv%instid(i)%der_trans(nchan,iv%instid(i)%nlevels,iv%instid(i)%num_rad))
end if

if (use_clddet_ecmwf) then
if (use_clddet==3) then
allocate (iv%instid(i)%kmin_t (iv%instid(i)%num_rad))
allocate (iv%instid(i)%kmax_p (iv%instid(i)%num_rad))
allocate (iv%instid(i)%sensitivity_ratio (nchan,iv%instid(i)%nlevels,iv%instid(i)%num_rad))
Expand Down
109 changes: 97 additions & 12 deletions var/da/da_radiance/da_cloud_detect_iasi.inc → var/da/da_radiance/da_cloud_detect.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
subroutine da_cloud_detect_iasi(isensor,nchannels,ndim,kts,kte,n,iv)
subroutine da_cloud_detect(isensor,nchannels,ndim,kts,kte,n,iv)

!** *CLOUD_DETECT_AIRS* - CLOUD FLAGGING FOR AIRS AND IASI
!** *CLOUD_DETECT* - CLOUD FLAGGING FOR IR Channels

! AUTHOR: THOMAS AULIGNE DATE : 01/08/2005
!
! PURPOSE.
! -------
! FLAG THE PRESENCE OF CLOUD CONTAMINATION IN AIRS AND IASI CHANNELS
! FLAG THE PRESENCE OF CLOUD CONTAMINATION IN IR CHANNELS
!
!** INTERFACE.
! ---------
Expand All @@ -21,10 +21,14 @@ subroutine da_cloud_detect_iasi(isensor,nchannels,ndim,kts,kte,n,iv)
!** EXTERNALS
! ---------
! N2QN1 - Minimization algorithm (double-precision constrained version of M1QN3)
!
!** METHODS
! CLOUD DETECTION SCHEME MMR FROM AULIGNÉ T. MWR (2014).OR. PF FROM XU ET AL., GMD (2016)
! MODIFICATIONS
! -------------
! NONE
! BY DONGMEI XU 201904
! PURPOSE, ADD CLOUD DETECTION METHOD BASED ON PARTICLE FILTER
! METHOD, Xu et al., 2016: A method for retrieving clouds with satellite infrared radiances using the particle filter. Geosci. Model Dev., 9, 3919–3932
!
!** -----------------------------------------

IMPLICIT NONE
Expand Down Expand Up @@ -76,9 +80,22 @@ DOUBLE PRECISION, ALLOCATABLE :: DZS(:)
INTEGER :: gn
REAL :: ZHOOK_HANDLE


logical :: iasi,airs, modis,imager,sounder,cris,giirs,ahi
double precision, allocatable :: ppx(:,:),wx(:),jo(:)
integer :: p1,ii,jj
double precision :: tmp
! Initializations
Band_Size(:) = 0
Bands(:,:) = 0


iasi = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'iasi'
airs = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'airs'
imager = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'imager'
ahi = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'ahi'
Band_Size(:) = 0
Bands(:,:) = 0

if ( iasi ) then
Band_Size(1:5) = (/ 193, 15, 116, 4, 15 /)

Bands(1:Band_Size(1),1) = &
Expand Down Expand Up @@ -129,6 +146,47 @@ REAL :: ZHOOK_HANDLE
& (/ 6982, 6985, 6987, 6989, 6991, 6993, 6995, 6997, 7267, 7269, &
& 7424, 7426, 7428, 7885, 8007 /)

else if (airs) then
Band_Size(1:5) = (/86, 0, 0, 16, 0 /)

Bands(1:Band_Size(1),1) = &
& (/ & !& 1, 6, 7, 10, 11, 15, 16, 17, 20, 21, &
& & !& 22, 24, 27, 28, 30, 36, 39, 40, 42, 51, &
& & !& 52, 54, 55, 56, 59, 62, 63, 68, 69, 71, &
& & !& 72, 73, 74, 75, 76, 77, 78, 79, 80, 82, &
& 92, 93, 98, 99, 101, 104, 105, & !& 83, 84, 86, 92, 93, 98, 99, 101, 104, 105, &
& 108, 110, 111, 113, 116, 117, 123, 124, 128, 129, &
& 138, 139, 144, 145, 150, 151, 156, 157, 159, 162, &
& 165, 168, 169, 170, 172, 173, 174, 175, 177, 179, &
& 180, 182, 185, 186, 190, 192, 198, 201, 204, & !& 180, 182, 185, 186, 190, 192, 193, 198, 201, 204, &
& 207, 210, 215, 216, 221, 226, 227, & !& 207, 210, 213, 215, 216, 218, 221, 224, 226, 227, &
& 232, 252, 253, 256, 257, 261, & !& 232, 239, 248, 250, 251, 252, 253, 256, 257, 261, &
& 262, 267, 272, 295, 299, 305, 310, & !& 262, 267, 272, 295, 299, 300, 305, 308, 309, 310, &
& 321, 325, 333, 338, 355, 362, 375, 453, 475, & !& 318, 321, 325, 333, 338, 355, 362, 375, 453, 475, &
& 484, 497, 528, 587, 672, 787, 791, 843, 870, 914, &
& 950 /)

Bands(1:Band_Size(4),4) = &
& (/ 1852, 1865, 1866, 1868, 1869, 1872, 1873, 1876, & !& 1852, 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876,
& 1881, 1882, 1883, 1911, 1917, 1918, & !& 1877, 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, &
& 1924, 1928 /) !& 1921, 1923, 1924, 1928, 1937 /)

else if (imager) then
Band_Size(1) = 2
Bands(1:Band_Size(1),1) = &
& (/4,5/)
else if (ahi) then
Band_Size(1) = 4
Bands(1:Band_Size(1),1) = &
& (/7,8,9,10/)
end if
jo=0
allocate(ppx(ndim*11,ndim))
allocate(wx(ndim*11))
allocate(jo(ndim*11))
wx=1.0
px(1:ndim-1) = 1.0/ndim
px(ndim) = 1.0 - SUM(px(1:ndim-1))
ichan = iv%instid(isensor)%ichan(1:nchannels)
rad_clr = iv%instid(isensor)%rad_xb(1:nchannels,n) !iv%instid(isensor)%tb_xb(1:nchan,n)
rad_obs = iv%instid(isensor)%rad_obs(1:nchannels,n) !iv%instid(isensor)%tb_inv(1:nchan,n) + rad_clr
Expand Down Expand Up @@ -163,17 +221,43 @@ REAL :: ZHOOK_HANDLE
nchan = nchan +1
AMAT(nchan,1:ndim-1) = rad_ovc(ich,1:NDIM-1) / rad_obs(ich)
AMAT(nchan,ndim) = rad_clr(ich) / rad_obs(ich)
ZF_CLR = ZF_CLR + 0.5*(AMAT(nchan,ndim)-1.0)**2
ZF_CLR = ZF_CLR + 0.5*(AMAT(nchan,ndim)-1.0)**2

if (use_clddet==2) then
p1=0
do ii=0,10
do jj=1,ndim-1
p1=p1+1
ppx(p1,1:ndim)=0 !initialization
ppx(p1,jj)=real(ii)/10
ppx(p1,ndim)=1-ppx(p1,jj) !initialization
jo(p1)=iv%instid(isensor)%tb_error(ich,n)*(1-SUM(ppx(p1,1:ndim)*amat(nchan,1:ndim)))**2
end do
end do
! step 2 calculate the weight
do ii=1,p1
! jo(ii)=jo(ii)-mv
tmp=exp(-jo(ii))/sum(exp(-jo(1:p1))) ! normalize the weight
wx(ii)=tmp*wx(ii) !1/jo !exp(-jo)
end do
wx(:)=wx(:)/sum(wx(1:p1))
px=0
do k = 1,ndim
px(k)= SUM(wx(1:p1)*ppx(1:p1,k))
end do
px(1:ndim) = px(1:ndim) / SUM(px(1:ndim)) ! Re-normalization
end if
ENDDO
IF (.NOT. LMATCH) then
if (print_detail_rad) then
write(unit=message(1),fmt='(A,2I8)') 'CLOUD_DETECT_IASI: No match for channel:',i,Bands(i,JBAND)
write(unit=message(1),fmt='(A,2I8)') 'CLOUD_DETECT: No match for channel:',i,Bands(i,JBAND)
call da_message(message(1:1))
endif
ENDIF
ENDDO
ENDDO BAND_LOOP ! Loop over band


if (use_clddet==1) then
!--------------------!
! Hessian evaluation !
!--------------------!
Expand Down Expand Up @@ -239,7 +323,8 @@ REAL :: ZHOOK_HANDLE
IF (ALLOCATED(ZRZ)) DEALLOCATE(ZRZ)
IF (ALLOCATED(DZS)) DEALLOCATE(DZS)
if (allocated(rzs)) deallocate(rzs)

end if !mmr

!-----------------!
! Cloudy radiance !
!-----------------!
Expand Down Expand Up @@ -268,4 +353,4 @@ REAL :: ZHOOK_HANDLE
#endif
end if

end subroutine da_cloud_detect_iasi
end subroutine da_cloud_detect
Loading