diff --git a/Registry/registry.var b/Registry/registry.var index 3405aa694e..b7147ffbe7 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -451,6 +451,7 @@ rconfig real varbc_factor namelist,wrfvar14 1 1.0 - "va rconfig integer varbc_nbgerr namelist,wrfvar14 1 5000 - "varbc_nbgerr" "" "" rconfig integer varbc_nobsmin namelist,wrfvar14 1 10 - "varbc_nobsmin" "" "" 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 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/da.make b/var/build/da.make index 70eebd2c71..2725792246 100644 --- a/var/build/da.make +++ b/var/build/da.make @@ -69,6 +69,7 @@ WRFVAR_OBJS = \ gsi_kinds.o \ gsi_constants.o \ gsi_thinning.o \ + mod_clddet_geoir.o \ da_wrfvar_io.o \ da_airsr.o \ da_wrfvar_top.o \ diff --git a/var/build/depend.txt b/var/build/depend.txt index 6f27f8f1f3..1b733744e2 100644 --- a/var/build/depend.txt +++ b/var/build/depend.txt @@ -110,7 +110,7 @@ da_blas.o : da_blas.f90 dgemv.inc dgemm.inc xerbla.inc lsame.inc dswap.inc dtrmm da_bogus.o : da_bogus.f90 da_calculate_grady_bogus.inc da_get_innov_vector_bogus.inc da_check_max_iv_bogus.inc da_transform_xtoy_bogus_adj.inc da_transform_xtoy_bogus.inc da_print_stats_bogus.inc da_oi_stats_bogus.inc da_residual_bogus.inc da_jo_and_grady_bogus.inc da_ao_stats_bogus.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_par_util1.o da_par_util.o da_interpolation.o da_define_structures.o da_control.o module_domain.o da_buoy.o : da_buoy.f90 da_calculate_grady_buoy.inc da_get_innov_vector_buoy.inc da_check_max_iv_buoy.inc da_transform_xtoy_buoy_adj.inc da_transform_xtoy_buoy.inc da_print_stats_buoy.inc da_oi_stats_buoy.inc da_residual_buoy.inc da_jo_and_grady_buoy.inc da_ao_stats_buoy.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_grid_definitions.o da_par_util.o da_par_util1.o da_interpolation.o da_define_structures.o da_control.o module_domain.o da_control.o : da_control.f90 module_driver_constants.o -da_crtm.o : da_crtm.f90 da_det_crtm_climat.inc da_crtm_sensor_descriptor.inc da_crtm_init.inc da_crtm_ad.inc da_crtm_direct.inc da_crtm_k.inc da_crtm_tl.inc da_get_innov_vector_crtm.inc da_transform_xtoy_crtm_adj.inc da_transform_xtoy_crtm.inc da_tracing.o da_tools.o da_tools_serial.o da_reporting.o da_radiance1.o module_dm.o da_interpolation.o da_control.o module_radiance.o da_define_structures.o module_domain.o +da_crtm.o : da_crtm.f90 da_det_crtm_climat.inc da_crtm_sensor_descriptor.inc da_crtm_init.inc da_crtm_ad.inc da_crtm_direct.inc da_crtm_k.inc da_crtm_tl.inc da_get_innov_vector_crtm.inc da_transform_xtoy_crtm_adj.inc da_transform_xtoy_crtm.inc da_tracing.o da_tools.o da_tools_serial.o da_reporting.o da_radiance1.o module_dm.o da_physics.o da_interpolation.o da_control.o module_radiance.o da_define_structures.o module_domain.o da_define_structures.o : da_define_structures.f90 da_gauss_noise.inc da_random_seed.inc da_initialize_cv.inc da_zero_vp_type.inc da_zero_y.inc da_zero_x.inc da_deallocate_y.inc da_deallocate_observations.inc da_deallocate_background_errors.inc da_allocate_y.inc da_allocate_observations.inc da_allocate_background_errors.inc da_wavelet.o da_reporting.o da_tools_serial.o da_tracing.o da_control.o module_domain.o da_allocate_y_rain.inc da_allocate_y_radar.inc da_allocate_observations_rain.inc da_allocate_obs_info.inc da_dynamics.o : da_dynamics.f90 da_wz_base.inc da_uv_to_vorticity.inc da_w_adjustment_adj.inc da_w_adjustment_lin.inc da_uv_to_divergence_adj.inc da_uv_to_divergence.inc da_psichi_to_uv_adj.inc da_psichi_to_uv.inc da_hydrostaticp_to_rho_lin.inc da_hydrostaticp_to_rho_adj.inc da_balance_geoterm_lin.inc da_balance_geoterm_adj.inc da_balance_equation_lin.inc da_balance_equation_adj.inc da_balance_cycloterm_lin.inc da_balance_cycloterm_adj.inc da_balance_cycloterm.inc da_wpec_constraint.inc da_wpec_constraint_adj.inc da_wpec_constraint_cycloterm.inc da_wpec_constraint_geoterm.inc da_wpec_constraint_lin.inc da_tools.o da_tracing.o da_ffts.o da_reporting.o da_define_structures.o module_comm_dm.o module_dm.o module_domain.o da_control.o da_divergence_constraint.inc da_divergence_constraint_adj.inc da_etkf.o : da_etkf.f90 da_solve_etkf.inc da_matmultiover.inc da_matmulti.inc da_innerprod.inc da_lapack.o da_gen_be.o da_control.o @@ -230,6 +230,7 @@ gen_spectra.o : gen_spectra.f90 da_tracing.o da_spectral.o da_control.o gsi_constants.o : gsi_constants.f90 init_constants_derived.inc gsi_kinds.o gsi_kinds.o : gsi_kinds.f90 gsi_thinning.o : gsi_thinning.f90 gsi_constants.o gsi_kinds.o +mod_clddet_geoir.o : mod_clddet_geoir.f90 da_control.o kma2netcdf.o : kma2netcdf.f90 module_kma_wave2grid.o module_kma2netcdf_interface.o da_wrfvar_io.o module_dm.o module_wrf_error.o module_timing.o module_driver_constants.o module_machine.o da_control.o module_domain.o module_configure.o module_ffts.o : module_ffts.f90 qpassm.inc fft661.inc fft551.inc module_kma2netcdf_interface.o : module_kma2netcdf_interface.f90 module_tiles.o module_configure.o module_driver_constants.o module_timing.o module_domain.o module_kma_wave2grid.o diff --git a/var/da/da_define_structures/da_define_structures.f90 b/var/da/da_define_structures/da_define_structures.f90 index a27a100220..5c117803f8 100644 --- a/var/da/da_define_structures/da_define_structures.f90 +++ b/var/da/da_define_structures/da_define_structures.f90 @@ -508,7 +508,12 @@ module da_define_structures real, pointer :: bgerr(:) 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(:,:,:) + end type clddet_geoir_type type cv_index_type integer :: ts integer :: nclouds @@ -548,6 +553,7 @@ module da_define_structures real, pointer :: satazi(:) real, pointer :: solzen(:) real, pointer :: solazi(:) + real, pointer :: tropt(:) !! Tropopause temperature, K. real, pointer :: t(:,:) real, pointer :: q(:,:) real, pointer :: mr(:,:) @@ -623,7 +629,7 @@ module da_define_structures real, pointer :: ice_coverage(:) real, pointer :: snow_coverage(:) integer, pointer :: crtm_climat(:) ! CRTM only - + type (clddet_geoir_type), pointer :: cld_qc(:) type (varbc_info_type) :: varbc_info type (varbc_type),pointer :: varbc(:) type (cv_index_type), pointer :: cv_index(:) diff --git a/var/da/da_physics/da_physics.f90 b/var/da/da_physics/da_physics.f90 index f138dafc4c..7f1391e4c0 100644 --- a/var/da/da_physics/da_physics.f90 +++ b/var/da/da_physics/da_physics.f90 @@ -107,6 +107,6 @@ module da_physics #include "da_integrat_dz.inc" #include "da_uv_to_sd_lin.inc" #include "da_uv_to_sd_adj.inc" - +#include "da_trop_wmo.inc" end module da_physics diff --git a/var/da/da_physics/da_trop_wmo.inc b/var/da/da_physics/da_trop_wmo.inc new file mode 100644 index 0000000000..a9c5936e76 --- /dev/null +++ b/var/da/da_physics/da_trop_wmo.inc @@ -0,0 +1,114 @@ +subroutine da_trop_wmo (t, z, p, nlev, tropt, tropp, tropk) + + !---------------------------------------------------------------------------- + ! * Computes tropopause T, P, and/or level based on code from Cameron Homeyer + ! and WMO definition + ! + ! * WMO tropopause definition: + ! The boundary between the troposphere and the stratosphere, where an + ! abrupt change in lapse rate usually occurs. It is defined as the lowest + ! level at which the lapse rate decreases to 2 °C/km or less, provided + ! that the average lapse rate between this level and all higher levels + ! within 2 km does not exceed 2 °C/km. + !---------------------------------------------------------------------------- + + implicit none + + ! Assumed shape inputs for single column (size=nlev) + ! ordered from bottom to top of model + real, intent(in) :: t(:) ! Temperature, K. (3D) + +!JJG: "z" is supposed to be height, not geopotential height. Does it matter? + real, intent(in) :: z(:) ! Geopotential height above m.s.l., m. + real, intent(in) :: p(:) ! Pressure, mb. + real, optional, intent(out) :: tropt ! Tropopause temperature, K. + real, optional, intent(out) :: tropp ! Tropopause pressure, mb. + integer, optional, intent(out) :: tropk + integer, intent(in) :: nlev + + real :: dtdz,laps !LAPS + integer :: dtdztest(nlev), ztest(nlev) + integer :: i, j, k, kk, ktrop + +! real, parameter :: tropz_min = 5000.0 +! real, parameter :: tropz_max = 19000.0 + + if (.not.present(tropt) .and. & + .not.present(tropp) .and. & + .not.present(tropk)) return + + if (present(tropt)) tropt = missing_r + if (present(tropp)) tropp = missing_r + + !Loop over levels to find tropopause (single column) + ktrop = nlev-1 + trop_loop: do k = 1, nlev-1 + if ( p(k) .le. 500.0 ) then + ! Compute lapse rate (-dT/dz) + dtdz = ( t(k+1) - t(k) ) / & + ( z(k) - z(k+1) ) + else + ! Set lapse rate for p > 500 hPa + dtdz = 999.9 + endif + !Check if local lapse rate <= 2 K/km + if (dtdz .le. 0.002) then + ! Initialize lapse rate and altitude test arrays + dtdztest = 0 + ztest = 0 + + ! Compute average lapse rate across levels above current candidate + do kk = k+1, nlev-1 + dtdz = ( t(kk+1) - t(k) ) / & + ( z(k) - z(kk+1) ) + + !If avg. lapse rate <= 2 K/km and z <= trop + 2 km, set pass flag + if ( ( dtdz .le. 0.002 ) .and. & + ( (z(k) - z(kk)) .le. 2000. ) ) THEN + dtdztest(kk) = 1 + endif + + ! If z <= trop + 2 km, set pass flag + IF ( (z(k) - z(kk)) .le. 2000.0 ) THEN + ztest(kk) = 1 + endif + enddo !kk loop + laps=dtdz !LAPS + IF (SUM(dtdztest) .eq. SUM(ztest)) THEN + ! If qualified as tropopause, set altitude index and return value + ktrop = k + exit trop_loop + ENDIF + ENDIF + end do trop_loop + +! ! Filter ktrop using tpause height thresholds +! ztest = 0 +! if ( z(ktrop) .gt. tropz_max ) then +! where ( z.le.tropz_max ) +! ztest = 1 +! end where +! do k = nlev, 1, -1 +! if (ztest(k) .eq. 1) then +! ktrop = k +! exit +! end if +! end do +! else if ( z(ktrop) .lt. tropz_min ) then +! where ( z.ge.tropz_min ) +! ztest = 1 +! end where +! do k = 1, nlev +! if (ztest(k) .eq. 1) then +! ktrop = k +! exit +! end if +! end do +! end if + + if (present(tropt)) tropt = t(ktrop) + if (present(tropp)) tropp = p(ktrop) + if (present(tropk)) tropk = ktrop + +end subroutine da_trop_wmo + diff --git a/var/da/da_radiance/da_allocate_rad_iv.inc b/var/da/da_radiance/da_allocate_rad_iv.inc index 9dc6c8e3a9..a34bbb0c65 100644 --- a/var/da/da_radiance/da_allocate_rad_iv.inc +++ b/var/da/da_radiance/da_allocate_rad_iv.inc @@ -10,7 +10,7 @@ subroutine da_allocate_rad_iv (i, nchan, iv) integer , intent (in) :: i, nchan type (iv_type) , intent (inout) :: iv - + integer :: n call da_trace_entry("da_allocate_rad_iv") allocate (iv%instid(i)%info%date_char(iv%instid(i)%num_rad)) @@ -104,7 +104,11 @@ subroutine da_allocate_rad_iv (i, nchan, iv) allocate (iv%instid(i)%satazi(iv%instid(i)%num_rad)) allocate (iv%instid(i)%solzen(iv%instid(i)%num_rad)) allocate (iv%instid(i)%solazi(iv%instid(i)%num_rad)) + allocate (iv%instid(i)%tropt(iv%instid(i)%num_rad)) allocate (iv%instid(i)%gamma_jacobian(nchan,iv%instid(i)%num_rad)) + if ( use_clddet_zz ) then + allocate (iv%instid(i)%cld_qc(iv%instid(i)%num_rad)) + 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)) diff --git a/var/da/da_radiance/da_crtm.f90 b/var/da/da_radiance/da_crtm.f90 index d183ced5d9..dd1a538842 100644 --- a/var/da/da_radiance/da_crtm.f90 +++ b/var/da/da_radiance/da_crtm.f90 @@ -32,7 +32,7 @@ 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, & + rtminit_satid, global,kms,kme,ims,ime,jms,jme,kts,kte,use_clddet, use_clddet_zz,& 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, & @@ -41,6 +41,7 @@ module da_crtm cloud_cv_options use da_interpolation, only : da_interp_lin_2d_partial,da_interp_lin_2d_adj_partial, & da_interp_2d_partial + use da_physics, only: da_trop_wmo use module_dm, only : wrf_dm_sum_real, wrf_dm_sum_reals use da_radiance1, only : da_biasprep,da_detsurtyp,da_biascorr, & da_biasprep,da_cld_eff_radius, da_mspps_emis, da_mspps_ts diff --git a/var/da/da_radiance/da_deallocate_radiance.inc b/var/da/da_radiance/da_deallocate_radiance.inc index fa5dec69eb..76d8249887 100644 --- a/var/da/da_radiance/da_deallocate_radiance.inc +++ b/var/da/da_radiance/da_deallocate_radiance.inc @@ -127,8 +127,11 @@ deallocate (iv%instid(i)%satazi) deallocate (iv%instid(i)%solzen) deallocate (iv%instid(i)%solazi) + deallocate (iv%instid(i)%tropt) deallocate (iv%instid(i)%gamma_jacobian) - + if ( use_clddet_zz ) then + deallocate (iv%instid(i)%cld_qc) + end if if (ANY(use_satcv)) then if (use_satcv(2)) then do n = 1,iv%instid(i)%num_rad 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 6bc85db10b..31b18b2ef0 100644 --- a/var/da/da_radiance/da_get_innov_vector_crtm.inc +++ b/var/da/da_radiance/da_get_innov_vector_crtm.inc @@ -93,7 +93,11 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) 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(:,:,:) + real :: geoht_pixel(kts:min(kte,kme-1)) + real :: tt_pixel(kts:min(kte,kme-1)) + real :: pp_pixel(kts:min(kte,kme-1)) Band_Size(1:5) = (/86, 0, 0, 16, 0 /) Bands(:,:) = 0 Bands(1:Band_Size(1),1) = & @@ -195,7 +199,17 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) else if ( n_vegtype == IGBP_n_type ) then wrf_to_crtm_mw = igbp_to_crtm_mw end if - + 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 !------------------------------------------------------ ! [1.0] calculate the background bright temperature !------------------------------------------------------- @@ -416,6 +430,16 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) end if end do + if (use_clddet_zz) then + ! Find tropopause temperature for Zhuge and Zou Cloud Detection + do k = kts, min(kte,kme-1) + tt_pixel(k) = Atmosphere(1)%Temperature(kte-k+1) + pp_pixel(k) = Atmosphere(1)%Pressure(kte-k+1) + call da_interp_2d_partial ( 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 + 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)) @@ -514,7 +538,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) end if if (Surface(1)%Land_Coverage > 0.0) then - Surface(1)%Land_Type = nint(model_ivgtyp) + Surface(1)%Land_Type = min(11,nint(model_ivgtyp)) Surface(1)%Vegetation_Type = max(1,wrf_to_crtm_mw(nint(model_ivgtyp))) Surface(1)%Soil_Type = max(1,wrf_to_crtm_soil(nint(model_isltyp))) ! glacial land ice soil type and glacial vegetation type are not valid land types @@ -920,7 +944,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) end do ! end loop for sensor deallocate (wrf_to_crtm_mw) - + if ( use_clddet_zz ) deallocate ( geoht_full ) call CRTM_Atmosphere_Destroy (Atmosphere) IF ( CRTM_Atmosphere_Associated( Atmosphere(1) ) ) & call da_error(__FILE__,__LINE__, & diff --git a/var/da/da_radiance/da_initialize_rad_iv.inc b/var/da/da_radiance/da_initialize_rad_iv.inc index a6d6a0892a..68aa250be6 100644 --- a/var/da/da_radiance/da_initialize_rad_iv.inc +++ b/var/da/da_radiance/da_initialize_rad_iv.inc @@ -100,8 +100,18 @@ subroutine da_initialize_rad_iv (i, n, iv, p) iv%instid(i)%satazi(n) = p%satazi iv%instid(i)%solzen(n) = p%solzen iv%instid(i)%solazi(n) = p%solazi + iv%instid(i)%tropt(n) = 0.0 ! iv%instid(i)%solazi(n) = 0.0 - + if ( use_clddet_zz ) then + iv%instid(i)%cld_qc(n)%RTCT = p % cld_qc % RTCT + iv%instid(i)%cld_qc(n)%RFMFT = p % cld_qc % RFMFT + iv%instid(i)%cld_qc(n)%TEMPIR = p % cld_qc % TEMPIR + iv%instid(i)%cld_qc(n)%tb_stddev_10 = p % cld_qc % tb_stddev_10 + iv%instid(i)%cld_qc(n)%tb_stddev_13 = p % cld_qc % tb_stddev_13 + iv%instid(i)%cld_qc(n)%tb_stddev_14 = p % cld_qc % tb_stddev_14 + iv%instid(i)%cld_qc(n)%terr_hgt = p % cld_qc % terr_hgt + iv%instid(i)%cld_qc(n)%CIRH2O = p % cld_qc % CIRH2O + end if if ( rtm_option == rtm_option_rttov ) then iv%instid(i)%surftype(n) = 0 iv%instid(i)%snow_frac(n) = 0.0 diff --git a/var/da/da_radiance/da_qc_ahi.inc b/var/da/da_radiance/da_qc_ahi.inc index f9411eaedf..0790c1c6e3 100644 --- a/var/da/da_radiance/da_qc_ahi.inc +++ b/var/da/da_radiance/da_qc_ahi.inc @@ -3,6 +3,8 @@ 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 + ! HISTORY: 2020/03/01 - Add clear sky cloud detection procedures Dongmei Xu, NUIST, NCAR/MMM !--------------------------------------------------------------------------- implicit none @@ -12,164 +14,407 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) integer, intent(in) :: nchan ! number of channel type (y_type), intent(in) :: ob ! Observation structure. type (iv_type), intent(inout) :: iv ! O-B structure. - ! local variables logical :: lmix, cloud_detection integer :: n,k,isflg,ios,fgat_rad_unit integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), & - nrej_omb_std(nchan),nrej_eccloud(nchan), & + nrej_omb_std(nchan), & nrej_clw(nchan),num_proc_domain, & nrej_mixsurface,nrej_land - real :: inv_grosscheck - + real :: inv_grosscheck, lowBTcheck + ! isflg: SEA(0),ICE(1),LAND(2),SNOW(3),MSEA(4),MICE(5),MLND(6),MSNO(7) + integer, parameter :: sea_flag = 0 + integer, parameter :: ice_flag = 1 + integer, parameter :: land_flag = 2 + integer, parameter :: snow_flag = 3 + integer, parameter :: msea_flag = 4 + integer, parameter :: mice_flag = 5 + integer, parameter :: mland_flag = 6 + integer, parameter :: msnow_flag = 7 character(len=30) :: filename - real :: c37_mean - + logical :: print_cld_debug + + !! Additional variables used by Zhuge and Zou (2017) + integer :: itest + logical :: reject_clddet + real :: crit_clddet + real :: rad_O14, rad_M14, rad_tropt + real :: rad_o_ch7, rad_b_ch7, rad_o_ch14, rad_b_ch14 + real :: Relaz, Glintzen, tb_temp1 + real :: wave_num(10) + real :: plbc1(10), plbc2(10) + real :: plfk1(10), plfk2(10) + integer, parameter :: num_clddet_tests = 10 + integer, parameter :: num_clddet_cats = 4 + real :: eps_clddet(num_clddet_tests+2,num_clddet_cats) + integer :: index_clddet(num_clddet_tests), offset_clddet + integer :: isflgs_clddet(num_clddet_cats) + logical :: qual_clddet(num_clddet_cats) + character(len=10) :: crit_names_clddet(num_clddet_tests) + integer :: nrej_clddet(nchan,num_clddet_tests+1) + integer*2 :: clddet_tests(num_clddet_tests) + + real, pointer :: tb_ob(:,:), tb_xb(:,:), tb_inv(:,:), tb_xb_clr(:,:), ca_mean(:,:) + integer :: tb_qc(nchan), tb_qc_clddet(nchan) + + real :: big_num + real, parameter :: C1=1.19104276e-5 ! = 2 * h * c**2 mWm-2sr-1(cm-1)-4 + real, parameter :: C2=1.43877516 ! = h * c / b = 1.43877 K(cm-1)-1 + + ! h = Planck's constant + ! b = Boltzmann constant + ! c = velocity of light + + integer, parameter :: ch7 = 1 + integer, parameter :: ch10 = 4 + integer, parameter :: ch13 = 7 + integer, parameter :: ch14 = 8 + integer, parameter :: ch15 = 9 ! 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 + integer, allocatable :: k_cloud_flag(:) ! cloud flags if (trace_use) call da_trace_entry("da_qc_ahi") - + ! These values can change as SRF (spectral response function) is updated + ! It is recommended to acquire these from L1B files, not copy them from GOES R PUG L1b Vol. 3 + wave_num(1:10) = (/2570.373, 1620.528, 1443.554, 1363.228, 1184.220, & + 1040.891, 968.001, 894.000, 815.294, 753.790/) + plbc1(1:10) = (/0.43361, 1.55228, 0.34427, 0.05651, 0.18733, & + 0.09102, 0.07550, 0.22516, 0.21702, 0.06266/) + plbc2(1:10) = (/0.99939, 0.99667, 0.99918, 0.99986, 0.99948, & + 0.99971, 0.99975, 0.99920, 0.99916, 0.99974/) + + plfk1 = C1 * wave_num**3 + plfk2 = C2 * wave_num + + crit_names_clddet(1) = "rtct" !here + crit_names_clddet(2) = "etrop" + crit_names_clddet(3) = "pfmft" + crit_names_clddet(4) = "nfmft" + crit_names_clddet(5) = "rfmft" !here + crit_names_clddet(6) = "cirh2o" !here + crit_names_clddet(7) = "emiss4" + crit_names_clddet(8) = "ulst" !here + crit_names_clddet(9) = "notc" + crit_names_clddet(10) = "tempir" !here + + big_num = huge(big_num) + !! Table 4 from Zhuge X. and Zou X. JAMC, 2016. [modified from ABI Cloud Mask Algorithm] + !ocean land snow ice (assume same as snow) + eps_clddet = transpose( reshape( (/ & + 3.2, 4.1, big_num, big_num & + , 0.1, 0.3, 0.4, 0.4 & + , 0.8, 2.5, big_num, big_num & + , 1.0, 2.0, 5.0, 5.0 & + , 0.7, 1.0, big_num, big_num & + , 0.7, 0.7, 0.7, 0.7 & + , 0.1, 0.2, 0.3, 0.3 & ! Land values: 0.46 in ABI CM; 0.2 in ZZ16 + , 2.86, big_num, big_num, big_num & + , 0.05, 0.1, 0.12, 0.12 & + , 15., 21., 10., 10. & + , 11., 15., 4.5, 4.5 & + , 2.0, 2.0, 2.0, 2.0 & + /), (/ size(eps_clddet, 2), size(eps_clddet, 1) /)) ) + index_clddet = (/1, 2, 3, 4, 5, 6, 7, 9, 10, 12/) + isflgs_clddet = (/sea_flag, land_flag, snow_flag, ice_flag/) + print_cld_debug = .false. ngood(:) = 0 nrej(:) = 0 nrej_omb_abs(:) = 0 nrej_omb_std(:) = 0 - nrej_eccloud(:) = 0 nrej_clw(:) = 0 nrej_mixsurface = 0 nrej_land = 0 + nrej_clddet(:,:)= 0 num_proc_domain = 0 - + tb_ob => ob%instid(i)%tb + tb_xb => iv%instid(i)%tb_xb + tb_inv => iv%instid(i)%tb_inv do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2 if (iv%instid(i)%info%proc_domain(1,n)) & num_proc_domain = num_proc_domain + 1 - if ( crtm_cloud ) then - ! calculate c37_mean - c37_mean = 1.0-(ob%instid(i)%tb(11,n)-ob%instid(i)%tb(12,n)+ & - iv%instid(i)%tb_xb(11,n)-iv%instid(i)%tb_xb(12,n))/ & - (2.0*(iv%instid(i)%tb_xb_clr(11,n)-iv%instid(i)%tb_xb_clr(12,n))) - end if - ! 0.0 initialise QC by flags assuming good obs !----------------------------------------------------------------- - iv%instid(i)%tb_qc(:,n) = qc_good + tb_qc = qc_good ! 1.0 reject all channels over mixture surface type !------------------------------------------------------ isflg = iv%instid(i)%isflg(n) lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7) if (lmix) then - iv%instid(i)%tb_qc(:,n) = qc_bad + tb_qc = qc_bad if (iv%instid(i)%info%proc_domain(1,n)) & nrej_mixsurface = nrej_mixsurface + 1 end if 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 if (only_sea_rad) then - iv%instid(i)%tb_qc(k,n) = qc_bad - nrej_land = nrej_land + 1 + tb_qc(k) = qc_bad + nrej_land = nrej_land + 1 + end if end if - !end if end do end if - - ! 3.0 check iuse + ! 2.0 check iuse !----------------------------------------------------------------- do k = 1, nchan if (satinfo(i)%iuse(k) .eq. -1) & - iv%instid(i)%tb_qc(k,n) = qc_bad + tb_qc(k) = qc_bad end do - ! 4.0 check cloud + ! 3.0 check clw in fg !----------------------------------------------------------------- if (.not. crtm_cloud ) then - do k = 1, nchan - + do k = 1, nchan if (iv%instid(i)%clwp(n) >= 0.2) then - iv%instid(i)%tb_qc(k,n) = qc_bad + tb_qc = qc_bad if (iv%instid(i)%info%proc_domain(1,n)) & nrej_clw(k) = nrej_clw(k) + 1 end if + end do + end if + ! METHOD,Zhuge, X. and Zou, X., Test of a Modified Infrared-Only ABI Cloud Mask Algorithm for AHI Radiance Observations. J. Appl. Meteor. Climatol., 2016, 55: 2529–2546. + ahi_clddet: if ( use_clddet_zz & + .and. all(tb_inv( (/ch7,ch14,ch15/), n ) .gt. missing_r) & + .and. all(tb_ob( (/ch7,ch14,ch15/), n ) .gt. missing_r) & + .and. all(tb_xb( (/ch7,ch14,ch15/), n ) .gt. missing_r) & + ) then + !!=============================================================================== + !!=============================================================================== + !! + !! 4.0 AHI IR-only Cloud Mask Algorithm, combines: + !! (*) Heidinger A. and Straka W., ABI Cloud Mask, version 3.0, 11 JUN, 2013. + !! (*) Zhuge X. and Zou X. JAMC, 2016. + !! + !!=============================================================================== + !!=============================================================================== + +!JJGDEBUG + if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG1: ', n, & + tb_inv(:,n) + if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG2: ', n, & + tb_xb(:,n) + if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG3: ', n, & + 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) + end if - cloud_detection=.false. - if (cloud_detection) then - if (iv%instid(i)%landsea_mask(n) == 0 ) then - if (iv%instid(i)%tb_xb(3,n)-ob%instid(i)%tb(3,n)>3.5) then - iv%instid(i)%tb_qc(k,n) = qc_bad - if (iv%instid(i)%info%proc_domain(1,n)) & - nrej_eccloud(k) = nrej_eccloud(k) + 1 - end if - else - if (iv%instid(i)%tb_xb(3,n)-ob%instid(i)%tb(3,n)>2.5) then - iv%instid(i)%tb_qc(k,n) = qc_bad - if (iv%instid(i)%info%proc_domain(1,n)) & - nrej_eccloud(k) = nrej_eccloud(k) + 1 - end if - end if + if (print_cld_debug) write(stdout,'(A,I8,8F12.4,2x,A)') 'PIXEL_DEBUG5: ', n, & + iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n), & + iv%instid(i)%satzen(n), iv%instid(i)%satazi(n), & + iv%instid(i)%solzen(n), iv%instid(i)%solazi(n), & + iv%instid(i)%tropt(n), iv%instid(i)%cld_qc(n)%terr_hgt, & + iv%instid(i)%info%date_char(n) +!JJGDEBUG + + clddet_tests = 0 + + if ( tb_xb(ch14,n) .gt. 0. .and. iv%instid(i)%tropt(n) .gt. 0. ) then + tb_temp1 = tb_ob(ch14,n) + rad_O14 = plfk1(ch14) / & + ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1 ) ) -1 ) + tb_temp1 = tb_xb(ch14,n) + rad_M14 = plfk1(ch14) / & + ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1) ) -1 ) + tb_temp1 = iv%instid(i)%tropt(n) + rad_tropt = plfk1(ch14) / & + ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1) ) -1 ) else - if (iv%instid(i)%cloudflag(n) < 3) then ! only use confident clear pixels - iv%instid(i)%tb_qc(:,n) = qc_bad - if (iv%instid(i)%info%proc_domain(1,n)) & - nrej_eccloud(k) = nrej_eccloud(k) + 1 - end if + rad_O14 = missing_r + rad_M14 = missing_r + rad_tropt = missing_r end if - - end do - end if - ! assigning obs errors - if (.not. crtm_cloud ) then - do k = 1, nchan - if (use_error_factor_rad) then - iv%instid(i)%tb_error(k,n) = & - satinfo(i)%error_std(k)*satinfo(i)%error_factor(k) - else - iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k) - end if - end do ! nchan - - else !crtm_cloud - ! symmetric error model, Geer and Bauer (2011) - do k = 1, nchan - if (c37_mean.lt.0.05) then - iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k) - else if (c37_mean.ge.0.05.and.c37_mean.lt.0.5) then - iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)+ & - (c37_mean-0.05)*(satinfo(i)%error_cld(k)-satinfo(i)%error_std(k))/(0.5-0.05) - else - iv%instid(i)%tb_error(k,n)= satinfo(i)%error_cld(k) - end if - end do ! nchan + if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) then + tb_temp1 = tb_ob(ch7,n) + rad_o_ch7 = plfk1(ch7) / & + ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. ) + tb_temp1 = tb_xb(ch7,n) + rad_b_ch7 = plfk1(ch7) / & + ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. ) + tb_temp1 = tb_ob(ch14,n) + rad_o_ch14 = plfk1(ch7) / & + ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. ) + tb_temp1 = tb_xb(ch14,n) + rad_b_ch14 = plfk1(ch7) / & + ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. ) + else + rad_o_ch7 = missing_r + rad_b_ch7 = missing_r + rad_o_ch14 = missing_r + rad_b_ch14 = missing_r + end if - end if + tb_qc_clddet = tb_qc + + do itest = 1, num_clddet_tests + qual_clddet = .true. + offset_clddet = 0 + crit_clddet = missing_r + select case (itest) + case (1) + !-------------------------------------------------------------------------- + ! 4.1 Relative Thermal Contrast Test (RTCT) + !-------------------------------------------------------------------------- + crit_clddet = iv%instid(i)%cld_qc(n)%RTCT + qual_clddet(3:4) = .false. + case (2) + !-------------------------------------------------------------------------- + ! 4.2 Cloud check: step 1 + ! Emissivity at Tropopause Test (ETROP) + !-------------------------------------------------------------------------- + if ( tb_xb(ch14,n) .gt. 0. .and. iv%instid(i)%tropt(n) .gt. 0. ) & + crit_clddet = (rad_O14 - rad_M14) / (rad_tropt - rad_M14) + case (3) + !-------------------------------------------------------------------------- + ! 4.3 Cloud check: step 2 + ! Positive Fourteen Minus Fifteen Test (PFMFT) + !-------------------------------------------------------------------------- + ! See ABI Cloud Mask Description for qual_clddet + qual_clddet = (tb_xb(ch14,n).ge.tb_xb(ch15,n)) + if ( (tb_inv(ch14,n) + tb_xb(ch14,n)).le.310. .and. & + iv%instid(i)%cld_qc(n)%tb_stddev_14.ge.0.3 .and. & + tb_ob(ch14,n).gt.0. .and. tb_ob(ch15,n).gt.0. ) & + crit_clddet = ( tb_ob(ch14,n) - tb_ob(ch15,n) ) +! above using ob without VarBC +! ------------------------------- +! crit_clddet = (tb_inv(ch14,n) + tb_xb(ch14,n) - & +! (tb_inv(ch15,n) + tb_xb(ch15,n)) )-& +! (tb_xb(ch14,n) - tb_xb(ch15,n)) * & +! (tb_ob(ch14,n) - 260.) / (tb_xb(ch14,n) - 260.) +! above using ob with VarBC +! ------------------------------- + if ( crit_clddet.gt.missing_r .and. & + (tb_inv(ch14,n) + tb_xb(ch14,n)).gt.270. .and. & + tb_xb(ch14,n).gt.270. ) & + crit_clddet = crit_clddet - & + (tb_xb(ch14,n) - tb_xb(ch15,n)) * & + (tb_ob(ch14,n) - 260.) / (tb_xb(ch14,n) - 260.) +! above 1 line using ob without VarBC +! (tb_inv(ch14,n) + tb_xb(ch14,n) - 260.)/ & +! (tb_xb(ch14,n) - 260.) +! above 2 lines using ob with VarBC + + case (4) + !-------------------------------------------------------------------------- + ! 4.4 Negative Fourteen Minus Fifteen Test (NFMFT) + !-------------------------------------------------------------------------- + if (tb_ob(ch14,n) .gt. 0. .and. tb_ob(ch15,n) .gt. 0.) & + crit_clddet = tb_inv(ch15,n) - tb_inv(ch14,n) + case (5) + !-------------------------------------------------------------------------- + ! 4.5 Relative Fourteen Minus Fifteen Test (RFMFT) + !-------------------------------------------------------------------------- + ! See ABI Cloud Mask Description for qual_clddet + qual_clddet = ( tb_ob(ch14,n) - tb_ob(ch15,n) ) .lt. 1.0 + qual_clddet(2) = qual_clddet(2) .and. tb_ob(ch14,n) .le. 300. + qual_clddet(3:4) = .false. + + crit_clddet = iv%instid(i)%cld_qc(n)%RFMFT + case (6) + !-------------------------------------------------------------------------- + ! 4.6 Cirrus Water Vapor Test (CIRH2O) + !-------------------------------------------------------------------------- + ! See ABI Cloud Mask Description for qual_clddet + qual_clddet = & + iv%instid(i)%cld_qc(n)%terr_hgt .le. 2000. & + .and. iv%instid(i)%cld_qc(n)%terr_hgt .ge. 0. & + .and. iv%instid(i)%cld_qc(n)%tb_stddev_10 .gt. 0.5 & + .and. iv%instid(i)%cld_qc(n)%tb_stddev_14 .gt. 0.5 + crit_clddet = iv%instid(i)%cld_qc(n)%CIRH2O + case (7) + !-------------------------------------------------------------------------- + ! 4.7 Modified 4um Emissivity Test (M-4EMISS) + !-------------------------------------------------------------------------- + ! Modify EMISS for sun glint area may be not work, because we are at north land + ! - compute relative azimuth + Relaz = RELATIVE_AZIMUTH(iv%instid(i)%solazi(n),iv%instid(i)%satazi(n)) + + ! - compute glint angle + Glintzen = GLINT_ANGLE(iv%instid(i)%solzen(n),iv%instid(i)%satzen(n),Relaz ) + + if ( Glintzen.lt.40.0 .and. isflg==sea_flag) then + crit_clddet = - tb_inv(ch7,n) ! (B_ch7 - O_ch7) + offset_clddet = 1 + else + if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) & + crit_clddet = (rad_o_ch7/rad_o_ch14 - rad_b_ch7/rad_b_ch14)/ & + (rad_b_ch7 / rad_b_ch14) + end if + case (8) + !-------------------------------------------------------------------------- + ! 4.8 Uniform low stratus Test (ULST) + !-------------------------------------------------------------------------- +!JJG, AHI error: Changed this to solzen instead of solazi for night/day test + qual_clddet = iv%instid(i)%solzen(n) >= 85.0 + if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) & + crit_clddet = rad_b_ch7/rad_b_ch14 - rad_o_ch7/rad_o_ch14 + case (9) + !-------------------------------------------------------------------------- + ! 4.9 New Optically Thin Cloud Test (N-OTC) + !-------------------------------------------------------------------------- +!JJG, AHI error: Changed this to solzen instead of solazi for night/day test + if ( iv%instid(i)%solzen(n) .ge. 85.0 ) & + offset_clddet = 1 ! night time + + if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch15,n) .gt. 0.) & +! using ob without VarBC +! ------------------------------- + crit_clddet = tb_ob(ch7,n) - tb_ob(ch15,n) + +! using ob with VarBC +! ------------------------------- +! crit_clddet = tb_inv(ch7,n) + tb_xb(ch7,n) - & +! (tb_inv(ch15,n) + tb_xb(ch15,n)) + case (10) + !-------------------------------------------------------------------------- + ! 4.10 Temporal Infrared Test (TEMPIR) + !-------------------------------------------------------------------------- + crit_clddet = iv%instid(i)%cld_qc(n)%TEMPIR + case default + cycle + end select + +! call evaluate_clddet_test ( & +! isflg, isflgs_clddet, crit_clddet, eps_clddet(index_clddet(itest)+offset_clddet,:), qual_clddet, & +! iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n), & +! reject_clddet ) + + reject_clddet = crit_clddet.gt.missing_r .and. & + any( isflg.eq.isflgs_clddet .and. & + crit_clddet.gt.eps_clddet(index_clddet(itest)+offset_clddet,:) .and. & + qual_clddet ) + + if (reject_clddet) then + tb_qc_clddet = qc_bad + 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) + end if - ! 5.0 check innovation - !----------------------------------------------------------------- - if (.not. crtm_cloud ) then - ! absolute departure check - do k = 1, nchan - inv_grosscheck = 8.0 - if (use_satcv(2)) inv_grosscheck = 100.0 - if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then - iv%instid(i)%tb_qc(k,n) = qc_bad - if (iv%instid(i)%info%proc_domain(1,n)) & - nrej_omb_abs(k) = nrej_omb_abs(k) + 1 - end if - end do ! nchan + clddet_tests(itest) = 1 + end if + end do + if (.not. crtm_cloud ) tb_qc = tb_qc_clddet - ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016) +!JJGDEBUG + if (print_cld_debug) write(stdout,'(A,I8,*(2x,I1:))') 'PIXEL_DEBUG6: ', n, clddet_tests +!JJGDEBUG + else ! not clddet_zz + ! 4. 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 + 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 @@ -202,40 +447,83 @@ subroutine da_qc_ahi (it, i, nchan, ob, 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 + if (iv%instid(i)%cloud_flag(k,n) == qc_bad) tb_qc(k) = qc_bad end do - + end if ahi_clddet + ! ---------------------------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 + if (.not. crtm_cloud ) then + do k = 1, nchan + if (use_error_factor_rad) then + iv%instid(i)%tb_error(k,n) = & + satinfo(i)%error_std(k)*satinfo(i)%error_factor(k) + else + iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k) + end if + end do ! nchan + end if + + ! 6.0 gross 8k check -clr,sdobs-clr + !----------------------------------------------------------------- + if (.not. crtm_cloud ) then + ! absolute departure check + do k = 1, nchan + inv_grosscheck = 8.0 + if (use_satcv(2)) inv_grosscheck = 100.0 + if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then + tb_qc(k) = qc_bad + if (iv%instid(i)%info%proc_domain(1,n)) & + nrej_omb_abs(k) = nrej_omb_abs(k) + 1 + end if + end do ! nchan + if (use_clddet_zz) then + ! SDob cloud inhomogeneous check + if (iv%instid(i)%cld_qc(n)%tb_stddev_13 >= 2) then ! only use abs clear pixel + tb_qc = qc_bad + if (iv%instid(i)%info%proc_domain(1,n)) & + nrej_clddet(:,11)= nrej_clddet(:,11)+1 + end if + end if end if + ! 7.0 3std check + !----------------------------------------------------------------- do k = 1, nchan ! relative departure check if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then - iv%instid(i)%tb_qc(k,n) = qc_bad + tb_qc(k) = qc_bad if (iv%instid(i)%info%proc_domain(1,n)) & nrej_omb_std(k) = nrej_omb_std(k) + 1 + end if - - ! final QC decsion - if (iv%instid(i)%tb_qc(k,n) == qc_bad) then - iv%instid(i)%tb_error(k,n) = 500.0 - if (iv%instid(i)%info%proc_domain(1,n)) & + end do ! nchan + + !final QC decsion + 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 + iv%instid(i)%tb_error(k,n) = 500.0 + if (iv%instid(i)%info%proc_domain(1,n)) & nrej(k) = nrej(k) + 1 - else - if (iv%instid(i)%info%proc_domain(1,n)) & - ngood(k) = ngood(k) + 1 - end if - end do ! nchan - + else + if (iv%instid(i)%info%proc_domain(1,n)) & + ngood(k) = ngood(k) + 1 + end if + end do ! nchan end do ! end loop pixel ! Do inter-processor communication to gather statistics. call da_proc_sum_int (num_proc_domain) call da_proc_sum_int (nrej_mixsurface) call da_proc_sum_int (nrej_land) - call da_proc_sum_ints (nrej_eccloud) call da_proc_sum_ints (nrej_omb_abs) call da_proc_sum_ints (nrej_omb_std) call da_proc_sum_ints (nrej_clw) + do itest = 1, num_clddet_tests+1 + call da_proc_sum_ints (nrej_clddet(:,itest)) + end do call da_proc_sum_ints (nrej) call da_proc_sum_ints (ngood) @@ -256,11 +544,13 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string if(num_proc_domain > 0) write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface - write(fgat_rad_unit,'(a20,i7)') ' nrej_land = ', nrej_land - write(fgat_rad_unit,'(a20)') ' nrej_eccloud(:) = ' - write(fgat_rad_unit,'(10i7)') nrej_eccloud(:) + write(fgat_rad_unit,'(a20,i7)') ' nrej_land = ', nrej_land write(fgat_rad_unit,'(a20)') ' nrej_clw(:) = ' - write(fgat_rad_unit,'(10i7)') nrej_clw(:) + write(fgat_rad_unit,'(10i7)') nrej_clw(:) + do itest = 1, num_clddet_tests+1 + write(fgat_rad_unit,'(a20,i2,a2)') ' nrej_clddet',itest,"=" + write(fgat_rad_unit,'(10i7)') nrej_clddet(:,itest) + end do write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = ' write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:) write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = ' @@ -276,3 +566,40 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) if (trace_use) call da_trace_exit("da_qc_ahi") end subroutine da_qc_ahi + +function relative_azimuth ( sol_az ,sen_az ) + + implicit none + + real :: sol_az + real :: sen_az + real :: relative_azimuth + + relative_azimuth = abs(sol_az - sen_az) + if (relative_azimuth > 180.0) then + relative_azimuth = 360.0 - relative_azimuth + endif + relative_azimuth = 180.0 - relative_azimuth + +end function relative_azimuth + + +function glint_angle ( sol_zen , sat_zen , rel_az ) + !------------------------------------------------------------------------------------ + ! Glint angle (the angle difference between direct "specular" reflection off + ! the surface and actual reflection toward the satellite.) + !------------------------------------------------------------------------------------ + + implicit none + + real :: sol_zen + real :: sat_zen + real :: rel_az + real :: glint_angle + + glint_angle = cos(sol_zen * deg2rad) * cos(sat_zen * deg2rad) + & + sin(sol_zen * deg2rad) * sin(sat_zen * deg2rad) * cos(rel_az * deg2rad) + glint_angle = max(-1.0 , min( glint_angle ,1.0 )) + glint_angle = acos(glint_angle) / deg2rad + +end function glint_angle diff --git a/var/da/da_radiance/da_radiance.f90 b/var/da/da_radiance/da_radiance.f90 index 60eb490dc1..587a75348c 100644 --- a/var/da/da_radiance/da_radiance.f90 +++ b/var/da/da_radiance/da_radiance.f90 @@ -58,14 +58,14 @@ 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_goesimgobs, pi, earth_radius, satellite_height,use_clddet_zz #ifdef CRTM use da_crtm, only : da_crtm_init, da_get_innov_vector_crtm #endif use da_define_structures, only : maxmin_type, iv_type, y_type, jo_type, j_type, & bad_data_type, x_type, number_type, bad_data_type, & - airsr_type,info_type, model_loc_type, varbc_info_type, varbc_type + airsr_type,info_type, model_loc_type, varbc_info_type, varbc_type,clddet_geoir_type use da_interpolation, only : da_to_zk, da_to_zk_new use da_tools_serial, only : da_get_unit, da_free_unit use da_par_util1, only : da_proc_sum_int,da_proc_sum_ints @@ -98,6 +98,8 @@ module da_radiance dlat_grid,dlon_grid,thinning_grid, & makegrids,map2grids, & destroygrids + use mod_clddet_geoir, only: qc_SDob,qc_RTCT,qc_rtmt,qc_cwvt,qc_tit, & + find_std,calc_rtmt,calc_correlation implicit none diff --git a/var/da/da_radiance/da_radiance1.f90 b/var/da/da_radiance/da_radiance1.f90 index f74beb4df8..1327930964 100644 --- a/var/da/da_radiance/da_radiance1.f90 +++ b/var/da/da_radiance/da_radiance1.f90 @@ -5,7 +5,7 @@ module da_radiance1 !--------------------------------------------------------------------------- #if defined(RTTOV) || defined(CRTM) - use module_radiance, only : satinfo,q2ppmv,rttov_inst_name + use module_radiance, only : satinfo,q2ppmv,rttov_inst_name,deg2rad #ifdef CRTM use module_radiance, only : CRTM_Planck_Radiance, CRTM_Planck_Temperature #endif @@ -21,10 +21,10 @@ module da_radiance1 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, use_satcv, cv_size_domain, & - cv_size_domain_js, calc_weightfunc, deg_to_rad, rad_to_deg + cv_size_domain_js, calc_weightfunc, deg_to_rad, rad_to_deg,use_clddet_zz 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 + be_type, clddet_geoir_type use module_dm, only : wrf_dm_sum_real, wrf_dm_sum_integer use da_par_util, only : da_proc_stats_combine use da_par_util1, only : da_proc_sum_int,da_proc_sum_ints @@ -47,7 +47,7 @@ module da_radiance1 type (info_type) :: info type (model_loc_type) :: loc - + type (clddet_geoir_type) :: cld_qc integer :: ifgat, landsea_mask, rain_flag integer :: scanline, scanpos real :: satzen, satazi, solzen, solazi ! satellite and solar angles diff --git a/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc b/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc index c55b78d757..a4869e3d73 100644 --- a/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc +++ b/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc @@ -1,6 +1,6 @@ -subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) +subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) !-------------------------------------------------------- - ! Purpose: read in JAXA AHI Level-1 and Level-2 data in NETCDF4 format + ! Purpose: read in JAXA AHI Level-1 data in NETCDF4 format for different times ! and form innovation structure ! ! METHOD: use F90 sequantial data structure to avoid read the file twice @@ -10,6 +10,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) ! and deallocate sequential data structure ! ! HISTORY: 2016/10/22 - Creation Yuanbing Wang, NUIST/CAS, NCAR/NESL/MMM/DAS + ! HISTORY: 2020/03/01 - Add clear sky cloud detection procedures Dongmei Xu, NUIST, NCAR/MMM ! To be devoloped: 1.time information; 2.dimension sequence !------------------------------------------------------------------------------ @@ -17,7 +18,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) implicit none - character(len=*), intent(in) :: infile_tb, infile_clp + character(len=*), intent(in) :: infile_tb, infile_tbp type(iv_type), intent(inout) :: iv ! fixed parameter values @@ -34,8 +35,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) real(4), allocatable :: vlongitude(:) ! value for longitude real(4), allocatable :: tbb(:,:,:) ! tb for band 7-16 - real(4), allocatable :: sat_zenith(:,:) - integer(2), allocatable :: cloud_type(:,:) + real(4), allocatable :: sat_zenith(:,:),sat_azi(:,:),sol_azi(:,:),sol_zenith(:,:) real(r_kind),parameter :: tbmin = 50._r_kind real(r_kind),parameter :: tbmax = 550._r_kind @@ -54,6 +54,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) logical :: outside, outside_all, iuse integer :: i,j,k,l,m,n, ifile, landsea_mask + logical :: found, head_found, head_allocated ! Other work variables @@ -65,11 +66,24 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) character(80) :: filename1,filename2 integer :: nchan,nlongitude,nlatitude,ilongitude,ilatitude,ichannels integer :: LatDimID,LonDimID - integer :: latid,lonid,tbb_id,sazid,cltyid + integer :: latid,lonid,tbb_id,sazid,cltyid,ter_id integer :: nfile - character(80) :: fname_tb(nfile_max),fname_clp(nfile_max) + character(80) :: fname_tb(nfile_max),fname_tbp(nfile_max) logical :: fexist, got_clp_file - + +! for qc stuff + integer, parameter :: ch7 = 1 + integer, parameter :: ch10 = 4 + integer, parameter :: ch13 = 7 + integer, parameter :: ch14 = 8 + integer, parameter :: ch15 = 9 + real(4), allocatable :: tbbp(:,:,:) ! tb for band 7-16 + real(r_kind), allocatable :: tbb_temp1(:,:), tbb_temp2(:,:) + real(r_kind), allocatable :: rtct(:,:), rtct2(:,:),rtmt(:,:) + real(r_kind), allocatable :: tempir(:,:),cwvt1(:,:) + real(r_kind), allocatable :: ter_sat(:,:),SDob_10(:,:),SDob_13(:,:),SDob_14(:,:) + character(200) :: filenameStatic + type(clddet_geoir_type) ::cld_qc ! Allocatable arrays integer(i_kind),allocatable :: ptotal(:) real,allocatable :: in(:), out(:) @@ -123,24 +137,24 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) ! first check if ahi nc file is available filename1 = trim(infile_tb) - filename2 = trim(infile_clp) + filename2 = trim(infile_tbp) inquire (file=filename1, exist=fexist) if ( fexist ) then nfile = 1 fname_tb(nfile) = filename1 - fname_clp(nfile) = filename2 + fname_tbp(nfile) = filename2 else ! check if netcdf4 files are available for multiple input files ! here 0x is the input file sequence number ! do not confuse it with fgat time slot index do i = 1, nfile_max write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i - write(filename2,fmt='(A,A,I2.2,A)') trim(infile_clp),'-',i + write(filename2,fmt='(A,A,I2.2,A)') trim(infile_tbp),'-',i inquire (file=filename1, exist=fexist) if ( fexist ) then nfile = nfile + 1 fname_tb(nfile) = filename1 - fname_clp(nfile) = filename2 + fname_tbp(nfile) = filename2 else exit end if @@ -159,7 +173,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) num_ahi_local_local = 0 num_ahi_global_local = 0 - ! open NETCDF4 L1 file for read + ! open NETCDF4 L1 file for read for the current time iret = nf90_open(fname_tb(ifile), NF90_NOWRITE, ncid) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & @@ -221,6 +235,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) do j=1,nlatitude do i=1,nlongitude tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb + add_offset_tb + if( tbb(i,j,k) < tbmin .or. tbb(i,j,k) > tbmax ) tbb(i,j,k) = missing_r end do end do ! sample display @@ -245,32 +260,187 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) ! sample display write(unit=stdout,fmt=*)& 'satellite zenith angle(pixel=1,scan=1): ',sat_zenith(1,1) - + + + ! read array: satellite azimuth angle + ! read + iret = nf90_inq_varid(ncid, 'SAA', sazid) + allocate( sat_azi(nlongitude,nlatitude)) + iret = nf90_get_var(ncid, sazid, sat_azi) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__, & + (/"NETCDF4 read error for: satellite azimuth angle"/)) + endif + do j=1,nlatitude + do i=1,nlongitude + sat_azi(i,j)=sat_azi(i,j) * scale_factor_tb + end do + end do + ! sample display + write(unit=stdout,fmt=*)& + 'satellite azimuth angle(pixel=1,scan=1): ',sat_azi(1,1) + + ! read array: solar zenith angle + ! read + iret = nf90_inq_varid(ncid, 'SOZ', sazid) + allocate( sol_zenith(nlongitude,nlatitude)) + iret = nf90_get_var(ncid, sazid, sol_zenith) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__, & + (/"NETCDF4 read error for: satellite zenith angle"/)) + endif + do j=1,nlatitude + do i=1,nlongitude + sol_zenith(i,j)=sol_zenith(i,j) * scale_factor_tb + end do + end do + ! sample display + write(unit=stdout,fmt=*)& + 'solar zenith angle(pixel=1,scan=1): ',sol_zenith(1,1) + + + ! read array: solar azimuth angle + ! read + iret = nf90_inq_varid(ncid, 'SOA', sazid) + allocate( sol_azi(nlongitude,nlatitude)) + iret = nf90_get_var(ncid, sazid, sol_azi) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__, & + (/"NETCDF4 read error for: satellite azimuth angle"/)) + endif + do j=1,nlatitude + do i=1,nlongitude + sol_azi(i,j)=sol_azi(i,j) * scale_factor_tb + end do + end do + ! sample display + write(unit=stdout,fmt=*)& + 'solar azimuth angle(pixel=1,scan=1): ',sol_azi(1,1) ! close infile_tb file iret = nf90_close(ncid) - - !open infile_clp file - got_clp_file = .false. - iret = nf90_open(fname_clp(ifile), NF90_NOWRITE, ncid) - if ( iret == 0 ) then - got_clp_file = .true. + allocate( tbbp(nlongitude,nlatitude,nchan)) + allocate( ter_sat(nlongitude,nlatitude)) + allocate( SDob_10(nlongitude,nlatitude)) + allocate( SDob_13(nlongitude,nlatitude)) + allocate( SDob_14(nlongitude,nlatitude)) + allocate( rtct(nlongitude,nlatitude)) + allocate( rtct2(nlongitude,nlatitude)) + allocate( rtmt(nlongitude,nlatitude)) + allocate( tempir(nlongitude,nlatitude)) + allocate( cwvt1(nlongitude,nlatitude)) + allocate( tbb_temp1(nlongitude,nlatitude)) + allocate( tbb_temp2(nlongitude,nlatitude)) + + if (use_clddet_zz) then +! -------------------------- +! 0-call read static data: terrain and landmask ????? +! -------------------------- + ! open NETCDF4 L1 file for read for the previous time + iret = nf90_open(fname_tbp(ifile), NF90_NOWRITE, ncid) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__, & + (/"Cannot open NETCDF4 file "//trim(fname_tbp(ifile))/)) + cycle infile_loop + endif + + ! read array: tb for band 7-16 + ! read + + do k=1,nchan + iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id) + iret = nf90_get_var(ncid, tbb_id, tbbp(:,:,k)) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__, & + (/"NETCDF4 read error for: Brightness Temperature"/)) + endif + do j=1,nlatitude + do i=1,nlongitude + tbbp(i,j,k)=tbbp(i,j,k) * scale_factor_tb + add_offset_tb + if( tbbp(i,j,k) < tbmin .or. tbbp(i,j,k) > tbmax ) tbbp(i,j,k) = missing_r + end do + end do + ! sample display + write(unit=stdout,fmt=*)& + 'tbbp(pixel=1,scan=1,chan=',k,'): ', tbbp(1,1,k) + end do + + ! close infile_tbp file + iret = nf90_close(ncid) + +! -------------------------- +! 1-call read static data: terrain and landmask ????? +! -------------------------- + + + filenameStatic="static_ahi.nc" + iret = nf90_open(filenameStatic, nf90_NOWRITE, ncid) + if(iret /= 0)then + print*,"Cannot open NETCDF4 file "//trim(filenameStatic) endif + iret = nf90_inq_varid(ncid, "AhiTer", ter_id) + iret = nf90_get_var(ncid, ter_id, ter_sat(:,:)) + if(iret /= 0) print*,"NETCDF4 read error for: SatTer" + where (ter_sat == -999.) ter_sat = missing_r + iret = nf90_close(ncid) + ! print *, 'ter_sat',maxval(ter_sat), minval(ter_sat) + +!!2- ZZ test0-SDob + + tbb_temp1 = missing_r + tbb_temp1 = tbb(:,:,ch13) ! 10.4um + SDob_13=missing_r + call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_13) + !print *, 'SDob_13',maxval(SDob_13), minval(SDob_13) + + !!!!!!!!!! for ZZ test6 + tbb_temp1 = missing_r + tbb_temp1 = tbb(:,:,ch10) ! 7.3um + SDob_10=missing_r + call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_10) +!print *, 'SDob_10',maxval(SDob_10), minval(SDob_10) + tbb_temp1 = missing_r + tbb_temp1 = tbb(:,:,ch14) ! 11.2umum + SDob_14=missing_r + call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_14) +!print *, 'SDob_14',maxval(SDob_14), minval(SDob_14) + +!! ZZ test1- RTCT + tbb_temp1 = missing_r + tbb_temp1 = tbb(:,:,ch14) ! 11.2um + rtct=0.1 + call qc_RTCT(nlongitude,nlatitude,tbb_temp1,ter_sat,rtct,rtct2) + +!! ZZ test5-RFMFT + + tbb_temp1 = missing_r + tbb_temp2 = missing_r + tbb_temp1 = tbb(:,:,ch14); tbb_temp2 = tbb(:,:,ch15) ! 11.2um and 12.4um + rtmt=missing_r + call qc_rtmt(nlongitude,nlatitude,tbb_temp1,tbb_temp2,rtmt) +!print *, 'rtmt',maxval(rtmt), minval(rtmt) + +!! ZZ test6-CIRH2O + tbb_temp1 = missing_r + tbb_temp2 = missing_r + tbb_temp1 = tbb(:,:,ch10) + tbb_temp2 = tbb(:,:,ch14) ! 7.3um and 11.2um + cwvt1=missing_r + call qc_cwvt(nlongitude,nlatitude,tbb_temp1,tbb_temp2,cwvt1) +!print *, 'cwvt1',maxval(cwvt1), minval(cwvt1) + + +!! ZZ test10- TEMPIR + tbb_temp1 = missing_r + tbb_temp2 = missing_r + tbb_temp1 = tbbp(:,:,ch14); tbb_temp2 = tbb(:,:,ch14) ! 11.2um + where(tbb_temp1 /=missing_r .and. tbb_temp2 /= missing_r) tempir = tbb_temp1-tbb_temp2 + where(tbb_temp1 ==missing_r .or. tbb_temp2 == missing_r) tempir = missing_r +!print *, 'tempir',maxval(tempir), minval(tempir),nlongitude,nlatitude +end if - if ( got_clp_file ) then - ! read array: satellite zenith angle - rcode = nf90_inq_varid(ncid, 'CLTYPE', cltyid) - allocate( cloud_type(nlongitude,nlatitude)) - iret = nf90_get_var(ncid, cltyid, cloud_type) - if(iret /= 0)then - call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/)) - endif - ! sample display - write(unit=stdout,fmt=*)'cloud_type(pixel=1,scan=1): ',cloud_type(1,1) - - ! close infile_clp file - iret = nf90_close(ncid) - end if + + ! 2.0 Loop to read netcdf and assign information to a sequential structure !------------------------------------------------------------------------- ! Allocate arrays to hold data @@ -280,10 +450,9 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) p => head head_allocated = .true. end if - ! start scan_loop + scan_loop: do ilatitude=1, nlatitude - call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time) if ( obs_time < time_slots(0) .or. & obs_time >= time_slots(num_fgat_time) ) cycle scan_loop @@ -291,25 +460,30 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) if ( obs_time >= time_slots(ifgat-1) .and. & obs_time < time_slots(ifgat) ) exit end do - ! start fov_loop - fov_loop: do ilongitude=1, nlongitude - + fov_loop: do ilongitude=1, nlongitude if ( sat_zenith(ilongitude,ilatitude) > 65.0 ) cycle fov_loop num_ahi_file = num_ahi_file + 1 num_ahi_file_local = num_ahi_file_local + 1 info%lat = vlatitude(ilatitude) info%lon = vlongitude(ilongitude) - + if (use_clddet_zz) then + cld_qc%TEMPIR = tempir(ilongitude,ilatitude) + cld_qc%RTCT = rtct(ilongitude,ilatitude) + cld_qc%RFMFT = rtmt(ilongitude,ilatitude) + cld_qc%CIRH2O = cwvt1(ilongitude,ilatitude) + cld_qc%tb_stddev_10 = SDob_10(ilongitude,ilatitude) + cld_qc%tb_stddev_13 = SDob_13(ilongitude,ilatitude) + cld_qc%tb_stddev_14 = SDob_14(ilongitude,ilatitude) + cld_qc % terr_hgt=ter_sat(ilongitude,ilatitude) + end if call da_llxy (info, loc, outside, outside_all) if (outside_all) cycle fov_loop - num_ahi_global = num_ahi_global + 1 num_ahi_global_local = num_ahi_global_local + 1 ptotal(ifgat) = ptotal(ifgat) + 1 if (outside) cycle fov_loop ! No good for this PE - num_ahi_local = num_ahi_local + 1 num_ahi_local_local = num_ahi_local_local + 1 write(unit=info%date_char, & @@ -335,7 +509,6 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) cycle fov_loop end if end if - num_ahi_used = num_ahi_used + 1 data_all = missing_r @@ -344,40 +517,54 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) if( tb < tbmin .or. tb > tbmax ) tb = missing_r data_all(k)= tb end do - + ! 4.0 assign information to sequential radiance structure !-------------------------------------------------------------------------- - allocate ( p % tb_inv (1:nchan )) + allocate ( p % tb_inv (1:nchan )) p%info = info p%loc = loc p%landsea_mask = 1 p%scanpos = ilongitude !nint(sat_zenith(ilongitude,ilatitude))+1.001_r_kind ! p%satzen = sat_zenith(ilongitude,ilatitude) - p%satazi = 0 - p%solzen = 0 - p%solazi = 0 + p%satazi = sat_azi(ilongitude,ilatitude) + p%solzen = sol_zenith(ilongitude,ilatitude) + p%solazi = sol_azi(ilongitude,ilatitude) p%tb_inv(1:nchan) = data_all(1:nchan) p%sensor_index = inst p%ifgat = ifgat - p%cloudflag = cloud_type(ilongitude,ilatitude) - - allocate (p%next) ! add next data - p => p%next + p%cld_qc = cld_qc + allocate (p%next) ! add next data + p => p%next nullify (p%next) end do fov_loop end do scan_loop - write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_file : ',num_ahi_file_local write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_global : ',num_ahi_global_local write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_local : ',num_ahi_local_local - end do infile_loop + deallocate(data_all) ! Deallocate data arrays deallocate(vlatitude) deallocate(vlongitude) deallocate(tbb) deallocate(sat_zenith) - if( got_clp_file ) deallocate(cloud_type) + deallocate(sol_zenith) + deallocate(sat_azi) + deallocate(sol_azi) + if( use_clddet_zz ) deallocate(tbbp) + if( use_clddet_zz ) deallocate(ter_sat) + if( use_clddet_zz ) deallocate(SDob_10) + if( use_clddet_zz ) deallocate(SDob_13) + if( use_clddet_zz ) deallocate(SDob_14) + if( use_clddet_zz ) deallocate(rtct) + if( use_clddet_zz ) deallocate(rtct2) + if( use_clddet_zz ) deallocate(rtmt) + if( use_clddet_zz ) deallocate(tempir) + if( use_clddet_zz ) deallocate(cwvt1) + if( use_clddet_zz ) deallocate(tbb_temp1) + if( use_clddet_zz ) deallocate(tbb_temp2) + + end do infile_loop if (thinning .and. num_ahi_global > 0 ) then #ifdef DM_PARALLEL @@ -464,7 +651,6 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) end do end if ! End of thinning - iv%total_rad_pixel = iv%total_rad_pixel + num_ahi_used iv%total_rad_channel = iv%total_rad_channel + num_ahi_used*nchan @@ -519,3 +705,4 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_clp) if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_jaxa") end subroutine da_read_obs_netcdf4ahi_jaxa + diff --git a/var/da/da_radiance/da_setup_radiance_structures.inc b/var/da/da_radiance/da_setup_radiance_structures.inc index 2c19d2f883..4e92945fcb 100644 --- a/var/da/da_radiance/da_setup_radiance_structures.inc +++ b/var/da/da_radiance/da_setup_radiance_structures.inc @@ -234,7 +234,8 @@ subroutine da_setup_radiance_structures( grid, ob, iv ) end if if (data_format==3) then write(unit=stdout,fmt='(a)') 'Reading AHI data from JAXA NETCDF4' - call da_read_obs_netcdf4ahi_jaxa (iv, 'L1AHITBR', 'L2AHICLP') +! call da_read_obs_netcdf4ahi_jaxa (iv, 'L1AHITBR', 'L2AHICLP') + call da_read_obs_netcdf4ahi_jaxa (iv, 'L1AHITBR', 'L1AHITBRP') end if !if (data_format==4) then !filename = 'ahi' diff --git a/var/da/da_radiance/mod_clddet_geoir.f90 b/var/da/da_radiance/mod_clddet_geoir.f90 new file mode 100644 index 0000000000..e932d5dfd3 --- /dev/null +++ b/var/da/da_radiance/mod_clddet_geoir.f90 @@ -0,0 +1,338 @@ + module mod_clddet_geoir +! use netcdf +! use mod_para + use da_control, only: missing_r + implicit none + contains + + subroutine qc_SDob(nlongitude,nlatitude,tbb,SDob) + ! ------------QC: step 1 ---- + ! Kozo Okamoto, QJ, 2017 and Zhuge and Zou, JAMC, 2016 + ! 10.4wm(Okamoto, QJ, 2017),Chan13, 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14 + ! only 10.8 for AGRI, !!!cha12!!! + implicit none + integer, intent(in):: nlongitude,nlatitude + real, intent(in) :: tbb(nlongitude,nlatitude) + real, intent(inout) :: SDob(nlongitude,nlatitude) +! temp + integer, parameter :: npixs_sdob = 3*3 ! pixles + real :: temp_tb(npixs_sdob), temp_std + integer :: i, j + + SDob=missing_r +! call SD_ob(nlongitude,nlatitude,tbb, SDob) + do i=2, nlongitude-1 + do j=2, nlatitude-1 + temp_tb(1) = tbb(i-1,j-1) + temp_tb(2) = tbb(i, j-1) + temp_tb(3) = tbb(i+1,j-1) + temp_tb(4) = tbb(i-1,j) + temp_tb(5) = tbb(i, j) + temp_tb(6) = tbb(i+1,j) + temp_tb(7) = tbb(i-1,j+1) + temp_tb(8) = tbb(i, j+1) + temp_tb(9) = tbb(i+1,j+1) + call find_std(nlongitude,nlatitude,npixs_sdob, temp_tb, temp_std, missing_r, missing_r) + SDob(i,j) = temp_std + end do + end do + + end subroutine + + subroutine qc_RTCT(nlongitude,nlatitude,tbb,ter,rtct,rtct2) + ! ------------QC: step 2 ---- + ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14 + ! Relative thermal Contrast test + ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14, and terrain + ! only 10.8 for AGRI, !!!cha12!!! + implicit none + integer, intent(in):: nlongitude,nlatitude + real, intent(in) :: tbb(nlongitude,nlatitude) + real, intent(in) :: ter(nlongitude,nlatitude) + real, intent(out) :: rtct(nlongitude,nlatitude) + real, intent(inout) :: rtct2(nlongitude,nlatitude) +! real, intent(inout) :: rtct2(nlongitude,nlatitude) +! temp + integer, parameter :: npixs_rtct = 3*3 ! pixles + real :: temp_ter(npixs_rtct), temp_tb(npixs_rtct), temp_std + integer :: i, j, m, n + + rtct=missing_r ! all missing_r for variables are same +! call SD_ob(nlongitude,nlatitude,tbb, SDob) + do i=2, nlongitude-1 + do j=2, nlatitude-1 + temp_ter(1) = ter(i-1,j-1) + temp_ter(2) = ter(i, j-1) + temp_ter(3) = ter(i+1,j-1) + temp_ter(4) = ter(i-1,j ) + temp_ter(5) = ter(i, j ) + temp_ter(6) = ter(i+1,j ) + temp_ter(7) = ter(i-1,j+1) + temp_ter(8) = ter(i, j+1) + temp_ter(9) = ter(i+1,j+1) + temp_tb(1) = tbb(i-1,j-1 ) + temp_tb(2) = tbb(i, j-1 ) + temp_tb(3) = tbb(i+1,j-1 ) + temp_tb(4) = tbb(i-1,j ) + temp_tb(5) = tbb(i, j ) + temp_tb(6) = tbb(i+1,j ) + temp_tb(7) = tbb(i-1,j+1 ) + temp_tb(8) = tbb(i, j+1 ) + temp_tb(9) = tbb(i+1,j+1 ) + m = COUNT( temp_ter /= missing_r) + n = COUNT( temp_tb /= 0. ) + if ( n == npixs_rtct .and. m == npixs_rtct) then + call find_std(nlongitude,nlatitude,npixs_rtct, temp_ter, temp_std, missing_r, missing_r) + rtct(i,j) = ( maxval(temp_tb)-tbb(i,j) ) - 3*temp_std*7*0.001 + rtct2(i,j) = 3*temp_std*7*0.001 + else + rtct(i,j) = 1 !missing_r + end if + end do + end do + + end subroutine + + subroutine qc_rtmt(nlongitude,nlatitude,tb12,tb13,rtmt) + ! ------------QC: step 2 ---- + ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14 + ! Relative fourteen minus fifteen test for AHI + ! Relative thirteen minus twelve for AGRI + implicit none + integer, intent(in):: nlongitude,nlatitude + real, intent(in) :: tb12(nlongitude,nlatitude) + real, intent(in) :: tb13(nlongitude,nlatitude) + real, intent(out) :: rtmt(nlongitude,nlatitude) + integer, parameter :: npixs_rtmt = 11*11 ! pixles + real :: rtmt_tb12(npixs_rtmt),rtmt_tb13(npixs_rtmt) + real :: max_tb12,max_tb13 + integer :: nstar, nend + integer :: i, j, ij, m + RTMT = missing_r + do i=6, nlongitude-5 ! 11*11 + do j=6, nlatitude-5 ! 11*11 + ij=1 + do m=-5, 5 + nstar = (ij-1)*11+1 + nend = ij*11 + rtmt_tb12(nstar:nend) = tb12(i-5:i+5,j+m) + rtmt_tb13(nstar:nend) = tb13(i-5:i+5,j+m) + ij = ij + 1 + end do + call calc_rtmt(npixs_rtmt,rtmt_tb12,rtmt_tb13,max_tb12,max_tb13) + if (max_tb12/=missing_r .and. max_tb13/=missing_r) then + rtmt(i,j)=abs( (tb12(i,j)-tb13(i,j))- & + (max_tb12-max_tb13) ) + else + rtmt(i,j)=missing_r + end if + end do + end do + end subroutine + + subroutine qc_cwvt(nlongitude,nlatitude,tb10,tb12,cwvt1) + ! ------------QC: step 3 ---- + ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14 + ! Cirrus water vapor test for AHI using Ch10, Ch14 5*5pixel + ! Ch10, Ch12 for AGRI + implicit none + integer, intent(in):: nlongitude,nlatitude + real, intent(in) :: tb10(nlongitude,nlatitude) + real, intent(in) :: tb12(nlongitude,nlatitude) + real, intent(out) :: cwvt1(nlongitude,nlatitude) +! real, intent(out) :: cwvt2(nlongitude,nlatitude) + integer, parameter :: npixs_rtmt = 5*5 ! pixles + real :: cwvt_tb10(npixs_rtmt),cwvt_tb12(npixs_rtmt) + real :: cwvtmp + integer :: nstar, nend + integer :: i, j, ij, m + cwvt1 = missing_r + cwvtmp = 0 !missing_r + do i=3, nlongitude-2 ! 11*11 + do j=3, nlatitude-2 ! 11*11 + ij=1 + do m=-2, 2 + nstar = (ij-1)*5+1 + nend = ij*5 + cwvt_tb10(nstar:nend) = tb10(i-2:i+2,j+m) + cwvt_tb12(nstar:nend) = tb12(i-2:i+2,j+m) + ij = ij + 1 + end do + call calc_correlation (npixs_rtmt,cwvt_tb10,cwvt_tb12,cwvtmp) + cwvt1(i,j)=cwvtmp +! call calc_correlation2(npixs_rtmt,cwvt_tb10,cwvt_tb12,cwvt2(i,j)) + end do + end do + end subroutine + + subroutine qc_tit(nlongitude,nlatitude,tbb,tbb15m,tit,missing_r) + ! ------------QC: step 4 ---- + ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14 + ! Cirrus water vapor test for AHI using Ch14, Ch14min before + ! Ch12, Ch12 for AGRI + implicit none + integer, intent(in):: nlongitude,nlatitude + real, intent(in) :: tbb(nlongitude,nlatitude) + real, intent(in) :: tbb15m(nlongitude,nlatitude) + real, intent(in) :: missing_r + real, intent(out) :: tit(nlongitude,nlatitude) + integer :: nstar, nend + integer :: i, j, ij, m + + tit=tbb15m-tbb + + end subroutine + + + subroutine find_std(nlongitude,nlatitude,n, arr, std_dev, miss1, miss2) + implicit none + integer, intent(in):: nlongitude,nlatitude + integer, intent(in) :: n + real, intent(in) :: arr(n) + real, intent(in) :: miss1, miss2 + real, intent(out):: std_dev + real :: variance, avg + real :: arr_temp(n) + integer :: i, m + + avg = 0. + m=0 + do i=1, n + if(arr(i)/=miss1 .and. arr(i)/=miss2)then + avg = avg + arr(i) + arr_temp(m+1)=arr(i) + m=m+1 + end if + end do + + if (m>=6) then + avg=avg/m + variance=0. + do i=1,m + variance=variance+(arr_temp(i)-avg)**2 + end do + variance=variance/m + std_dev=sqrt(variance) + else + std_dev=missing_r + end if + end subroutine + + subroutine calc_rtmt(npixs_rtmt,temp_tb12,temp_tb13,max_tb12,max_tb13) + implicit none + integer, intent(in) :: npixs_rtmt + real, intent(in) :: temp_tb12(npixs_rtmt),temp_tb13(npixs_rtmt) + real, intent(out):: max_tb12, max_tb13 + real :: variance, avg + real :: arr_temp12(npixs_rtmt),arr_temp13(npixs_rtmt) + integer :: i, m + + m=0 + do i=1, npixs_rtmt + if( temp_tb12(i) /=missing_r .and. temp_tb12(i) /=missing_r .and. & + temp_tb13(i) /=missing_r .and. temp_tb13(i) /=missing_r )then + arr_temp12(m+1)=temp_tb12(i) + arr_temp13(m+1)=temp_tb13(i) + m=m+1 + end if + end do + ! 80 per 121 would be ok + if (m>=80) then + max_tb12=maxval(arr_temp12(1:m)) + max_tb13=maxval(arr_temp13(1:m)) + else + max_tb12=missing_r + max_tb13=missing_r + end if + end subroutine + + subroutine calc_correlation( n, arr1, arr2, r) + implicit none + + integer, intent(in) :: n + real, intent(in) :: arr1(n) + real, intent(in) :: arr2(n) + real, intent(out) :: r +! below variables can be output: m, t(significant) + integer, parameter :: lag = 0 + real :: t + integer :: m +! temp variables + real :: arr1_tem(n), arr2_tem(n) +! real, allocatable :: x(:) +! real, allocatable :: y(:) +! real, allocatable :: xdev(:) +! real, allocatable :: ydev(:) +! real, allocatable :: xdevydev(:) +! real, allocatable :: xdevxdev(:) +! real, allocatable :: ydevydev(:) + real :: x(3*n) + real :: y(3*n) + real :: xdev(3*n) + real :: ydev(3*n) + real :: xdevydev(3*n) + real :: xdevxdev(3*n) + real :: ydevydev(3*n) + real :: xmn + real :: ymn + real :: COVxy + real :: Sx + real :: Sy + integer :: i + + arr1_tem=arr1 + arr2_tem=arr2 +! allocate( x(1:3*N) ) +! allocate( y(1:3*N) ) + + x(:) = missing_r + y(:) = missing_r + x(N+1:2*N) = arr1_tem(:) + y(N+1:2*N) = arr2_tem(:) + + y = EOSHIFT( y, SHIFT = -lag, BOUNDARY = missing_r ) + + WHERE ( x == missing_r ) y = missing_r + WHERE ( y == missing_r ) x = missing_r + + m = COUNT( x /= missing_r ) + r = missing_r + if (m >= 16) then + + xmn = sum( x, dim = 1, mask = x /= missing_r ) / real(m) + ymn = sum( y, dim = 1, mask = y /= missing_r ) / real(m) +! allocate ( xdev(1:3*N) ) +! allocate ( ydev(1:3*N) ) + xdev(:) = x(:) - xmn + where ( x == missing_r ) xdev(:) = missing_r + ydev(:) = y(:) - ymn + where ( y == missing_r ) ydev(:) = missing_r +! allocate ( xdevydev(1:3*N) ) + xdevydev(:) = xdev(:) * ydev(:) + where ( x == missing_r ) xdevydev(:) = missing_r + + COVxy = sum( xdevydev, dim = 1, mask = xdevydev /= missing_r ) / real(m) + +! allocate ( xdevxdev(1:3*N) ) + xdevxdev(:) = xdev(:) * xdev(:) + where ( x == missing_r ) xdevxdev(:) = missing_r + Sx = sqrt( sum( xdevxdev, dim = 1, mask = xdevxdev /= missing_r ) / real(m) ) + +! allocate ( ydevydev(1:3*N) ) + ydevydev(:) = ydev(:) * ydev(:) + where ( y == missing_r ) ydevydev(:) = missing_r + Sy = SQRT( SUM( ydevydev, DIM = 1, MASK = ydevydev /= missing_r ) / REAL(m) ) + if ( Sx /= missing_r .and. Sy/= missing_r .and. (Sx * Sy/=0) ) r = COVxy / ( Sx * Sy ) + else + r = missing_r + end if +! deallocate( x, y ) + + end subroutine + + + + + + +end module