diff --git a/Registry/registry.var b/Registry/registry.var index 15013f9238..5524d7446f 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -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" "" "" diff --git a/var/build/depend.txt b/var/build/depend.txt index b40827b2e9..438b62951e 100644 --- a/var/build/depend.txt +++ b/var/build/depend.txt @@ -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 diff --git a/var/da/da_radiance/da_allocate_rad_iv.inc b/var/da/da_radiance/da_allocate_rad_iv.inc index 1c01c421d6..9dc6c8e3a9 100644 --- a/var/da/da_radiance/da_allocate_rad_iv.inc +++ b/var/da/da_radiance/da_allocate_rad_iv.inc @@ -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)) diff --git a/var/da/da_radiance/da_cloud_detect_iasi.inc b/var/da/da_radiance/da_cloud_detect.inc similarity index 70% rename from var/da/da_radiance/da_cloud_detect_iasi.inc rename to var/da/da_radiance/da_cloud_detect.inc index 05cc2fe57c..7cebadd00c 100644 --- a/var/da/da_radiance/da_cloud_detect_iasi.inc +++ b/var/da/da_radiance/da_cloud_detect.inc @@ -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. ! --------- @@ -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 @@ -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) = & @@ -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 @@ -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 ! !--------------------! @@ -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 ! !-----------------! @@ -268,4 +353,4 @@ REAL :: ZHOOK_HANDLE #endif end if -end subroutine da_cloud_detect_iasi +end subroutine da_cloud_detect diff --git a/var/da/da_radiance/da_cloud_detect_airs.inc b/var/da/da_radiance/da_cloud_detect_airs.inc deleted file mode 100644 index 8162bb07f3..0000000000 --- a/var/da/da_radiance/da_cloud_detect_airs.inc +++ /dev/null @@ -1,267 +0,0 @@ -subroutine da_cloud_detect_airs(isensor,nchannels,ndim,kts,kte,n,iv) - -!** *CLOUD_DETECT_AIRS* - CLOUD FLAGGING FOR AIRS AND IASI - -! AUTHOR: THOMAS AULIGNE DATE : 01/08/2005 -! -! PURPOSE. -! ------- -! FLAG THE PRESENCE OF CLOUD CONTAMINATION IN AIRS AND IASI CHANNELS -! -!** INTERFACE. -! --------- -! WHERE nchannels : Number of channels -! kts : model level corresponding to 100hPa (top of initial cloud search) -! kte : model level corresponding to surface (lower extent of cloud) -! rad_obs : Potentially cloudy observations -! rad_clr : Clear radiance from Model -! rad_ovc : Model overcast radiance estimates -! cloud_flag : Cloud flag by channel; 1=clear, -1=cloudy -! -!** EXTERNALS -! --------- -! N2QN1 - Minimization algorithm (double-precision constrained version of M1QN3) -! -! MODIFICATIONS -! ------------- -! NONE -!** ----------------------------------------- - -IMPLICIT NONE - -!* 0.1 Global arrays -INTEGER,INTENT(IN) :: isensor ! sensor index. -INTEGER,INTENT(IN) :: nchannels ! number of channels -INTEGER,INTENT(IN) :: ndim ! model levels between surface (lower extent of cloud) and 100hPa (top of cloud search) -INTEGER,INTENT(IN) :: kts ! model level corresponding to 100hPa (top of initial cloud search) -INTEGER,INTENT(IN) :: kte ! model level corresponding to surface (lower extent of cloud) -INTEGER,INTENT(IN) :: n ! pixel index -type (iv_type), intent(inout) :: iv ! O-B structure. - -INTEGER,PARAMETER :: NITER = 100 -INTEGER,PARAMETER :: NBAND = 1 -LOGICAL,PARAMETER :: LPRECON = .false. -INTEGER,PARAMETER :: NEIGNVEC = 4 -INTEGER,PARAMETER :: AIRS_Max_Channels = 2378 - -!! local declarations -INTEGER :: ichan(nchannels) ! AIRS and IASI channel IDs -REAL :: rad_obs(nchannels) ! Observed radiance -REAL :: rad_clr(nchannels) ! Model clear radiance estimates -REAL :: rad_ovc(nchannels,ndim-1) ! RT overcast radiance estimates -double precision :: px(ndim) !neignvec) ! Cloud fractions -REAL :: rad_cld(nchannels) -INTEGER :: ich,ilev,jlev,i,j,JBAND -double precision :: ZF, ZF_CLR -double precision :: ZG(ndim) -double precision :: binf(ndim), bsup(ndim) -REAL :: AMAT(nchannels,ndim) -INTEGER :: NCHAN -LOGICAL :: LMATCH -INTEGER :: Band_Size(5) -INTEGER :: Bands(AIRS_Max_Channels,5) -integer :: cldtoplevel - -! Hessian evaluation - REAL :: hessian(ndim,ndim), eignvec(ndim,ndim), eignval(ndim) - -!! local declarations for N2QN1 !! -INTEGER :: NRZ, impres, io, IMODE, NSIM, nit, izs(2) -double precision :: ZDF1, ZDXMIN, ZEPSG -double precision ,ALLOCATABLE :: ZRZ(:) -real, allocatable :: RZS(:) -INTEGER, ALLOCATABLE :: IZ(:) -DOUBLE PRECISION, ALLOCATABLE :: DZS(:) - -REAL :: ZHOOK_HANDLE - -! Initializations - Band_Size(:) = 0 - Bands(:,:) = 0 - 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(2),2) = & -!& (/ 1003, 1012, 1019, 1024, 1030, 1038, 1048, 1069, 1079, 1082, & -!& 1083, 1088, 1090, 1092, 1095, 1104, 1111, 1115, 1116, 1119, & -!& 1120, 1123, 1130, 1138, 1142, 1178, 1199, 1206, 1221, 1237, & -!& 1252, 1260, 1263, 1266, 1278, 1285 /) - -! Bands(1:Band_Size(3),3) = & -!& (/ 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, & !& 1290, 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, & -!& 1466, 1477, 1500, 1519, 1538, 1545, & !& 1466, 1471, 1477, 1479, 1488, 1500, 1519, 1520, 1538, 1545, & -!& 1565, 1574, 1583, 1593, 1627, 1636, 1652, 1669, & !& 1565, 1574, 1583, 1593, 1614, 1627, 1636, 1644, 1652, 1669, & -!& 1694, 1708, 1723, 1740, 1748, 1756, & !& 1674, 1681, 1694, 1708, 1717, 1723, 1740, 1748, 1751, 1756, & -!& 1766, 1771, 1777, 1783, 1794, 1800, 1806, & !& 1763, 1766, 1771, 1777, 1780, 1783, 1794, 1800, 1803, 1806, & -!& 1826, 1843 /) !& 1812, 1826, 1843 /) - - 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 /) - -! Bands(1:Band_Size(5),5) = & -!& (/ 1938, 1939, 1941, 1946, 1947, 1948, 1958, 1971, 1973, 1988, & -!& 1995, 2084, 2085, 2097, 2098, 2099, 2100, 2101, 2103, 2104, & -!& 2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, 2115, & -!& 2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123, 2128, 2134, & -!& 2141, 2145, 2149, 2153, 2164, 2189, 2197, 2209, 2226, 2234, & -!& 2280, 2318, 2321, 2325, 2328, 2333, 2339, 2348, 2353, 2355, & -!& 2363, 2370, 2371, 2377 /) - - 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 - rad_ovc = iv%instid(isensor)%rad_ovc(1:nchannels,kts+1:kte,n) - - nchan = 0 - AMAT(:,:) = 0.0 - px(1:ndim-1) = 0.0 - px(ndim) = 1.0 - ZF_CLR = 0.0 - nit = niter - -!do ich=1,nchannels -! CALL CRTM_Planck_Radiance(11,ichan(ich),tb_obs(ich),rad_obs(ich)) -! CALL CRTM_Planck_Radiance(11,ichan(ich),tb_clr(ich),rad_clr(ich)) -!end do - - !--------------------! - ! Loop over band ! - !--------------------! - BAND_LOOP: DO JBAND = 1, NBAND - DO i = 1, Band_Size(JBAND) - LMATCH = .FALSE. - DO ich=1,nchannels - IF (ichan(ich)/= Bands(i,JBAND)) CYCLE - IF ((rad_obs(ich)<=0.0).OR.(rad_obs(ich)>1000.0)) CYCLE - IF ((rad_clr(ich)<=0.0).OR.(rad_clr(ich)>1000.0)) CYCLE - IF (ANY(rad_ovc(ich,1:NDIM-1)<=0.0)) CYCLE - IF (ANY(rad_ovc(ich,1:NDIM-1)>1000.0)) CYCLE - - LMATCH = .TRUE. !! Found match for channel - 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 - ENDDO - IF (.NOT. LMATCH) WRITE(*,*) & - 'CLOUD_DETECT_AIRS: No matching for channel:',i,Bands(i,JBAND) - ENDDO - ENDDO BAND_LOOP ! Loop over band - - !--------------------! - ! Hessian evaluation ! - !--------------------! - IF (LPRECON) THEN -! write(76,*) '**** HESSIAN ****' - hessian(:,:)= 0.0 - DO ilev=1, NDIM - DO jlev=ilev, NDIM - DO J=1,NCHAN - hessian(ilev,jlev) = hessian(ilev,jlev) + & - (AMAT(J,ilev)-AMAT(J,NDIM)) * & - (AMAT(J,jlev)-AMAT(J,NDIM)) - ENDDO - hessian(jlev,ilev) = hessian(ilev,jlev) - ENDDO -! write(76,*) hessian(ilev,1:NDIM) - ENDDO -!!! call da_eof_decomposition(ndim, hessian, eignvec, eignval) - ENDIF - - !-----------------! - ! n2qn1 minimizer ! - !-----------------! - impres = 2 - io = 66 - NSIM = NITER+5 - ZDXMIN = 1.e-6 - ZEPSG = 1.e-3 !e-9 - IMODE = 1 - NRZ = NDIM*(NDIM+9)/2 ! N2QN1 - ALLOCATE(IZ(2*NDIM +1)) - ALLOCATE(ZRZ(NRZ)) - ALLOCATE(DZS(NCHAN*NDIM)) - allocate(rzs(ndim*neignvec)) - binf = -1000.0 - bsup = 1000.0 - izs(1) = nchan - izs(2) = neignvec - rzs = 0.0 - ZRZ = 0.0 - dzs(1:NCHAN*NDIM)=RESHAPE(AMAT(1:NCHAN,1:NDIM),(/NCHAN*NDIM/)) - - IF (LPRECON) THEN - IMODE = 2 - i = 0 - DO ilev=1, NDIM - DO jlev=ilev, NDIM - i = i + 1 - ZRZ(i) = hessian(jlev,ilev) - ENDDO - ENDDO - ENDIF -! rzs(1:ndim*neignvec) = RESHAPE(eignvec(1:ndim,1:neignvec),(/ndim*neignvec/)) -! rzs(ndim*neignvec+1:(ndim+1)*neignvec)= eignval(1:neignvec) - - call da_cloud_sim(0,NDIM,px,ZF,ZG,izs,RZS,DZS) - ZDF1 = 1.e-1*ZF - - - call da_error(__FILE__,__LINE__, & - (/"inria_n2qn1 is not implemented here, please contact the author of this subroutine."/)) -! call inria_n2qn1(da_cloud_sim,NDIM,px,ZF,ZG,(/(ZDXMIN,jlev=1,NDIM)/),ZDF1, & -! ZEPSG,impres,io,IMODE,nit,NSIM,binf,bsup,IZ,ZRZ,izs,RZS,DZS) - - IF (ALLOCATED(IZ)) DEALLOCATE(IZ) - IF (ALLOCATED(ZRZ)) DEALLOCATE(ZRZ) - IF (ALLOCATED(DZS)) DEALLOCATE(DZS) - if (allocated(rzs)) deallocate(rzs) - - !-----------------! - ! Cloudy radiance ! - !-----------------! - DO ich=1,nchannels - rad_cld(ich) = SUM(px(1:ndim-1) * rad_ovc(ich,1:ndim-1)) + px(ndim) * rad_clr(ich) - - if (ABS(rad_cld(ich)-rad_clr(ich)) < 0.01*rad_clr(ich)) then - iv%instid(isensor)%cloud_flag(ich,n) = qc_good - else - iv%instid(isensor)%cloud_flag(ich,n) = qc_bad - end if - ENDDO - - ! Dump cloud top pressure - do ilev = kte, kts+2, -1 - if (px(ilev-kts+1) > 0.01) cldtoplevel = ilev - end do - - if (rtm_option == rtm_option_rttov) then -#ifdef RTTOV - iv%instid(isensor)%clwp(n) = coefs(isensor)%coef%ref_prfl_p(cldtoplevel) -#endif - elseif (rtm_option == rtm_option_crtm) then -#ifdef CRTM - iv%instid(isensor)%clwp(n) = iv%instid(isensor)%pm(cldtoplevel,n) -#endif - end if - - -end subroutine da_cloud_detect_airs diff --git a/var/da/da_radiance/da_crtm.f90 b/var/da/da_radiance/da_crtm.f90 index e73a1d548d..d183ced5d9 100644 --- a/var/da/da_radiance/da_crtm.f90 +++ b/var/da/da_radiance/da_crtm.f90 @@ -32,11 +32,11 @@ module da_crtm biasprep, qc_rad,missing_r,rtminit_sensor,rtminit_nsensor, filename_len, & use_error_factor_rad,read_biascoef, analysis_date,time_window_max, & time_window_min,num_fgat_time,rtminit_platform, print_detail_rad, & - rtminit_satid, global,kms,kme,ims,ime,jms,jme,kts,kte,use_clddet_mmr, & + rtminit_satid, global,kms,kme,ims,ime,jms,jme,kts,kte,use_clddet, & use_crtm_kmatrix, use_varbc, freeze_varbc, use_pseudo_rad, & use_antcorr, time_slots, use_satcv, use_simulated_rad, simulated_rad_io, & simulated_rad_ngrid, interp_option, use_mspps_emis, use_mspps_ts, calc_weightfunc, & - use_clddet_ecmwf,its,ite,jts,jte, & + its,ite,jts,jte, & crtm_coef_path, crtm_irwater_coef, crtm_mwwater_coef, crtm_irland_coef, crtm_visland_coef, & cloud_cv_options use da_interpolation, only : da_interp_lin_2d_partial,da_interp_lin_2d_adj_partial, & diff --git a/var/da/da_radiance/da_get_innov_vector_crtm.inc b/var/da/da_radiance/da_get_innov_vector_crtm.inc index 403337a44a..6bc85db10b 100644 --- a/var/da/da_radiance/da_get_innov_vector_crtm.inc +++ b/var/da/da_radiance/da_get_innov_vector_crtm.inc @@ -86,7 +86,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) logical :: calc_tb_clr type (CRTM_RTSolution_type), allocatable :: RTSolution_clr(:,:) - ! Initializations for AIRS (MMR) Cloud Detection + ! Initializations for Cloud Detection integer :: ilev, jlev, nclouds real, allocatable :: hessian(:,:) real*8, allocatable :: eignvec(:,:), eignval(:) @@ -616,9 +616,9 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) end if - ! Compute Overcast Radiances for AIRS Cloud Detection (MMR) +! Compute Overcast Radiances for MMR in Auligné (2014).or. PF in Xu et al., (2016).or. MW03 ecmwf method in McNally AP and Watts PD (2003) !---------------------------------------------------------------- - if (use_clddet_mmr .or. use_clddet_ecmwf ) then + if (use_clddet==1 .or. use_clddet==2.or. use_clddet==3 ) then do i = 1, nchanl !crtm_2.2.3 contains Upwelling_Overcast_Radiance iv%instid(inst)%rad_ovc(i,kts:kte,n)=RTSolution(i,1)%Upwelling_Overcast_Radiance(:) @@ -773,7 +773,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) end do end if - if (use_clddet_ecmwf) then !begin ecmwf cloud detection + if (use_clddet==3) then !begin ecmwf cloud detection !kmin_t,kmax_p corresponding tropopause and KBPL ! find the k index of tropopause t_tropp = 999.0 ! initialize diff --git a/var/da/da_radiance/da_initialize_rad_iv.inc b/var/da/da_radiance/da_initialize_rad_iv.inc index e3efab0098..a6d6a0892a 100644 --- a/var/da/da_radiance/da_initialize_rad_iv.inc +++ b/var/da/da_radiance/da_initialize_rad_iv.inc @@ -147,7 +147,7 @@ subroutine da_initialize_rad_iv (i, n, iv, p) iv%instid(i)%der_trans(:,:,n) = 0.0 end if - if (use_clddet_ecmwf) then + if (use_clddet==3) then iv%instid(i)%kmin_t(n) = 0.0 iv%instid(i)%kmax_p(n) = 0.0 iv%instid(i)%sensitivity_ratio(:,:,n) = 0.0 diff --git a/var/da/da_radiance/da_qc_ahi.inc b/var/da/da_radiance/da_qc_ahi.inc index 91a0328c5a..f9411eaedf 100644 --- a/var/da/da_radiance/da_qc_ahi.inc +++ b/var/da/da_radiance/da_qc_ahi.inc @@ -3,7 +3,6 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) !--------------------------------------------------------------------------- ! Purpose: perform quality control for ahi data. ! Method: Assume cloud flag coming from GEOCAT processing - ! To be developed: built in cloud_detection method !--------------------------------------------------------------------------- implicit none @@ -27,6 +26,13 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) character(len=30) :: filename real :: c37_mean + ! mmr or pf Cloud Detection Variables + integer :: kts_100hPa(1), kte_surf, ndim + integer :: numrad_local(nchan), numrad_global(nchan) + real :: tstore + real :: bias_local(nchan), bias_global(nchan) + integer :: kmin, kmax + integer, allocatable :: k_cloud_flag(:) ! cloud flags if (trace_use) call da_trace_entry("da_qc_ahi") ngood(:) = 0 @@ -160,6 +166,45 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) nrej_omb_abs(k) = nrej_omb_abs(k) + 1 end if end do ! nchan + + ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016) + !--------------------------------------------- + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then + iv%instid(i)%cloud_flag(:,n) = qc_good + + if (rtm_option == rtm_option_rttov) then +#ifdef RTTOV + kte_surf = iv%instid(i)%nlevels + kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), & + MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0) + do k=1,nchan + tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * & + (ob%instid(i)%tb(k,n) - bias_global(k)) + iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / & + (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0) + end do +#endif + elseif (rtm_option == rtm_option_crtm) then + kte_surf = kte + kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), & + MASK = iv%instid(i)%pm(kts:kte,n) < 100.0) + + do k = 1, nchan + CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), & + iv%instid(i)%rad_obs(k,n)) + end do + + end if + + ndim = kte_surf - kts_100hPa(1) + 1 + + call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) + end if + + do k = 1, nchan + if (iv%instid(i)%cloud_flag(k,n) == qc_bad) iv%instid(i)%tb_qc(k,n) = qc_bad + end do + end if do k = 1, nchan diff --git a/var/da/da_radiance/da_qc_airs.inc b/var/da/da_radiance/da_qc_airs.inc index 024567ee85..50a682508f 100644 --- a/var/da/da_radiance/da_qc_airs.inc +++ b/var/da/da_radiance/da_qc_airs.inc @@ -266,9 +266,9 @@ subroutine da_qc_airs (it, i, nchan, ob, iv) end do ! chan - ! k. ecmwf cloud_detection + ! k. MW03 ecmwf method in McNally AP and Watts PD (2003) !----------------------------------------------------------- - if (use_clddet_ecmwf) then + if (use_clddet==3) then iv%instid(i)%cloud_flag(:,n) = qc_good allocate ( k_cloud_flag(1:nchan) ) call da_error(__FILE__,__LINE__, & @@ -310,7 +310,7 @@ subroutine da_qc_airs (it, i, nchan, ob, iv) ! Do inter-processor communication to gather statistics. - if (use_clddet_mmr .and. (.not.use_satcv(2))) then + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then do k = 1, nchan bias_global(k) = wrf_dm_sum_real(bias_local(k)) numrad_global(k) = wrf_dm_sum_integer(numrad_local(k)) @@ -320,9 +320,9 @@ subroutine da_qc_airs (it, i, nchan, ob, iv) do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2 - ! 1. Cloud detection scheme (NESDIS or MMR) + ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016) !--------------------------------------------- - if (use_clddet_mmr .and. (.not.use_satcv(2))) then + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then iv%instid(i)%cloud_flag(:,n) = qc_good if (rtm_option == rtm_option_rttov) then @@ -351,7 +351,7 @@ subroutine da_qc_airs (it, i, nchan, ob, iv) ndim = kte_surf - kts_100hPa(1) + 1 - call da_cloud_detect_airs(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) + call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) end if do k = 1, nchan diff --git a/var/da/da_radiance/da_qc_goesimg.inc b/var/da/da_radiance/da_qc_goesimg.inc index 91809d3836..eff20e0bf6 100644 --- a/var/da/da_radiance/da_qc_goesimg.inc +++ b/var/da/da_radiance/da_qc_goesimg.inc @@ -31,6 +31,13 @@ subroutine da_qc_goesimg(it, i, nchan, ob, iv) character(len=30) :: filename + ! mmr or pf Cloud Detection Variables + integer :: kts_100hPa(1), kte_surf, ndim + integer :: numrad_local(nchan), numrad_global(nchan) + real :: tstore + real :: bias_local(nchan), bias_global(nchan) + integer :: kmin, kmax + integer, allocatable :: k_cloud_flag(:) ! cloud flags if (trace_use_dull) call da_trace_entry("da_qc_goesimg.inc") ngood(:) = 0 @@ -101,7 +108,45 @@ subroutine da_qc_goesimg(it, i, nchan, ob, iv) ! nrej_eccloud = nrej_eccloud + 1 ! end if !end if - end if + + ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016) + !--------------------------------------------- + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then + iv%instid(i)%cloud_flag(:,n) = qc_good + + if (rtm_option == rtm_option_rttov) then +#ifdef RTTOV + kte_surf = iv%instid(i)%nlevels + kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), & + MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0) + do k=1,nchan + tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * & + (ob%instid(i)%tb(k,n) - bias_global(k)) + iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / & + (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0) + end do +#endif + else if (rtm_option == rtm_option_crtm) then + kte_surf = kte + kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), & + MASK = iv%instid(i)%pm(kts:kte,n) < 100.0) + + do k = 1, nchan + CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), & + iv%instid(i)%rad_obs(k,n)) + end do + + end if + + ndim = kte_surf - kts_100hPa(1) + 1 + + call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) + end if + + do k = 1, nchan + if (iv%instid(i)%cloud_flag(k,n) == qc_bad) iv%instid(i)%tb_qc(k,n) = qc_bad + end do + end if !not crtm ! c. check innovation !----------------------------------------------------------- diff --git a/var/da/da_radiance/da_qc_iasi.inc b/var/da/da_radiance/da_qc_iasi.inc index 6d5c31b723..44f22e853f 100644 --- a/var/da/da_radiance/da_qc_iasi.inc +++ b/var/da/da_radiance/da_qc_iasi.inc @@ -134,9 +134,9 @@ subroutine da_qc_iasi (it, i, nchan, ob, iv) end if end do ! chan - ! e. ecmwf detection - !---------------------------------------- - if (use_clddet_ecmwf) then + ! k. MW03 ecmwf method in McNally AP and Watts PD (2003) + !----------------------------------------------------------- + if (use_clddet==3) then iv%instid(i)%cloud_flag(:,n) = qc_good allocate ( k_cloud_flag(1:nchan) ) call da_error(__FILE__,__LINE__, & @@ -180,7 +180,7 @@ subroutine da_qc_iasi (it, i, nchan, ob, iv) ! Do inter-processor communication to gather statistics. - if (use_clddet_mmr .and. (.not.use_satcv(2))) then + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then do k = 1, nchan bias_global(k) = wrf_dm_sum_real(bias_local(k)) numrad_global(k) = wrf_dm_sum_integer(numrad_local(k)) @@ -190,9 +190,9 @@ subroutine da_qc_iasi (it, i, nchan, ob, iv) do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2 - ! 1. Cloud detection scheme (NESDIS or MMR) + ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016) !--------------------------------------------- - if (use_clddet_mmr .and. (.not.use_satcv(2))) then + if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then iv%instid(i)%cloud_flag(:,n) = qc_good if (rtm_option == rtm_option_rttov) then @@ -221,7 +221,7 @@ subroutine da_qc_iasi (it, i, nchan, ob, iv) ndim = kte_surf - kts_100hPa(1) + 1 - call da_cloud_detect_iasi(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) + call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv) end if do k = 1, nchan diff --git a/var/da/da_radiance/da_radiance1.f90 b/var/da/da_radiance/da_radiance1.f90 index 1b1f29a86d..ba7ca04c72 100644 --- a/var/da/da_radiance/da_radiance1.f90 +++ b/var/da/da_radiance/da_radiance1.f90 @@ -20,8 +20,8 @@ module da_radiance1 rtm_option_rttov,rtm_option_crtm, radiance, only_sea_rad, & global, gas_constant, gravity, monitor_on,kts,kte,use_rttov_kmatrix, & use_pseudo_rad, pi, t_triple, crtm_cloud, DT_cloud_model,write_jacobian, & - use_crtm_kmatrix,use_clddet_mmr, use_satcv, cv_size_domain, & - cv_size_domain_js, calc_weightfunc, use_clddet_ecmwf, deg_to_rad, rad_to_deg + use_crtm_kmatrix,use_clddet, use_satcv, cv_size_domain, & + cv_size_domain_js, calc_weightfunc, deg_to_rad, rad_to_deg 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 @@ -223,8 +223,7 @@ module da_radiance1 #include "da_qc_crtm.inc" #endif #include "da_cloud_sim.inc" -#include "da_cloud_detect_airs.inc" -#include "da_cloud_detect_iasi.inc" +#include "da_cloud_detect.inc" #include "da_qc_airs.inc" #include "da_qc_amsua.inc" #include "da_qc_amsub.inc" diff --git a/var/da/da_radiance/da_rttov.f90 b/var/da/da_radiance/da_rttov.f90 index 9150a87a4e..46e71c55b5 100644 --- a/var/da/da_radiance/da_rttov.f90 +++ b/var/da/da_radiance/da_rttov.f90 @@ -30,7 +30,7 @@ module da_rttov rtminit_print, rttov_scatt,comm,ierr,biasprep, qc_rad, & num_fgat_time,stdout,trace_use, use_error_factor_rad, & qc_good, qc_bad,myproc,biascorr, global,ims,ime,jms,jme, & - use_clddet_mmr, time_slots, rttov_emis_atlas_ir, rttov_emis_atlas_mw, & + use_clddet, time_slots, rttov_emis_atlas_ir, rttov_emis_atlas_mw, & use_mspps_emis, use_mspps_ts use da_interpolation, only : da_to_zk_new, & da_interp_lin_2d, da_interp_lin_3d, da_interp_lin_3d_adj, da_interp_lin_2d_adj