diff --git a/Registry/registry.var b/Registry/registry.var index b7147ffbe7..516239f451 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -452,6 +452,7 @@ rconfig integer varbc_nbgerr namelist,wrfvar14 1 5000 - "va 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 integer ahi_superob_halfwidth namelist,wrfvar14 1 0 - "ahi_superob_halfwidth" "" "" rconfig logical airs_warmest_fov namelist,wrfvar14 1 .false. - "airs_warmest_fov" "" "" rconfig logical use_satcv namelist,wrfvar14 2 .false. - "use_satcv" "" "" rconfig logical use_blacklist_rad namelist,wrfvar14 1 .true. - "use_blacklist_rad" "" "" diff --git a/var/da/da_define_structures/da_define_structures.f90 b/var/da/da_define_structures/da_define_structures.f90 index 5c117803f8..bb0e2cefaa 100644 --- a/var/da/da_define_structures/da_define_structures.f90 +++ b/var/da/da_define_structures/da_define_structures.f90 @@ -514,6 +514,10 @@ module da_define_structures real :: CIRH2O !real, allocatable :: CIRH2O(:,:,:) end type clddet_geoir_type + type superob_type + real, allocatable :: tb_obs(:,:) + type(clddet_geoir_type), allocatable :: cld_qc(:) + end type superob_type type cv_index_type integer :: ts integer :: nclouds @@ -549,6 +553,7 @@ module da_define_structures integer, pointer :: cloud_flag(:,:) integer, pointer :: cloudflag(:) integer, pointer :: rain_flag(:) + real, allocatable :: cloud_frac(:) real, pointer :: satzen(:) real, pointer :: satazi(:) real, pointer :: solzen(:) @@ -629,11 +634,12 @@ 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(:) + integer :: superob_width = 1 type (varbc_info_type) :: varbc_info type (varbc_type),pointer :: varbc(:) type (cv_index_type), pointer :: cv_index(:) type (infa_type) :: info + type (superob_type), allocatable :: superob(:,:) end type instid_type type iv_type diff --git a/var/da/da_monitor/da_rad_diags.f90 b/var/da/da_monitor/da_rad_diags.f90 index 2f8ebdff7d..57d3c4fda5 100644 --- a/var/da/da_monitor/da_rad_diags.f90 +++ b/var/da/da_monitor/da_rad_diags.f90 @@ -58,8 +58,8 @@ program da_rad_diags integer, dimension(:), allocatable :: landsea_mask, soiltyp, vegtyp real*4, dimension(:), allocatable :: lat, lon, elv, elev real*4, dimension(:), allocatable :: ret_clw - real*4, dimension(:), allocatable :: satzen, satazi, t2m, mr2m, u10, v10, ps, ts - real*4, dimension(:), allocatable :: smois, tslb, snowh, vegfra, clwp + real*4, dimension(:), allocatable :: satzen, satazi, solzen, solazi, t2m, mr2m, u10, v10, ps, ts + real*4, dimension(:), allocatable :: smois, tslb, snowh, vegfra, clwp, cloud_frac integer, dimension(:,:), allocatable :: tb_qc real*4, dimension(:,:), allocatable :: tb_obs, tb_bak, tb_inv, tb_oma, tb_err, ems, ems_jac real*4, dimension(:,:), allocatable :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water @@ -230,6 +230,8 @@ program da_rad_diags allocate ( lon(1:total_npixel) ) allocate ( satzen(1:total_npixel) ) allocate ( satazi(1:total_npixel) ) + allocate ( solzen(1:total_npixel) ) + allocate ( solazi(1:total_npixel) ) allocate ( ret_clw(1:total_npixel) ) !obs retrieved clw allocate ( t2m(1:total_npixel) ) allocate ( mr2m(1:total_npixel) ) @@ -246,6 +248,7 @@ program da_rad_diags allocate ( vegfra(1:total_npixel) ) allocate ( elev(1:total_npixel) ) allocate ( clwp(1:total_npixel) ) !model/guess clwp + allocate ( cloud_frac(1:total_npixel) ) allocate ( tb_obs(1:nchan,1:total_npixel) ) allocate ( tb_bak(1:nchan,1:total_npixel) ) allocate ( tb_inv(1:nchan,1:total_npixel) ) @@ -341,14 +344,15 @@ program da_rad_diags npixel_loop: do ipixel = ips, ipe - read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i6,i3,f6.0,4f8.2,f8.3)',iostat=ios) & + read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i6,i3,f6.0,6f8.2,f8.3)',iostat=ios) & n, datestr2(ipixel), scanpos(ipixel), landsea_mask(ipixel), elv(ipixel), & - lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel), ret_clw(ipixel) + lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel), solzen(ipixel), & + solazi(ipixel), ret_clw(ipixel) - read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2)',iostat=ios) & + read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2,f15.5)',iostat=ios) & t2m(ipixel), mr2m(ipixel), u10(ipixel), v10(ipixel), ps(ipixel), ts(ipixel), & smois(ipixel), tslb(ipixel), snowh(ipixel), isflg(ipixel), soiltyp(ipixel), & - vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel) + vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel), cloud_frac(ipixel) read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OBS read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_obs(:,ipixel) read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! BAK @@ -580,6 +584,8 @@ program da_rad_diags ios = NF_DEF_VAR(ncid, 'lon', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'satzen', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'satazi', NF_FLOAT, 1, ishape(1), varid) + ios = NF_DEF_VAR(ncid, 'solzen', NF_FLOAT, 1, ishape(1), varid) + ios = NF_DEF_VAR(ncid, 'solazi', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 't2m', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'mr2m', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'u10', NF_FLOAT, 1, ishape(1), varid) @@ -598,6 +604,7 @@ program da_rad_diags if ( amsr2 ) then ios = NF_DEF_VAR(ncid, 'ret_clw', NF_FLOAT, 1, ishape(1), varid) end if + ios = NF_DEF_VAR(ncid, 'cloud_frac', NF_FLOAT, 1, ishape(1), varid) ios = NF_ENDDEF(ncid) ! @@ -768,6 +775,10 @@ program da_rad_diags ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), satzen) ios = NF_INQ_VARID (ncid, 'satazi', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), satazi) + ios = NF_INQ_VARID (ncid, 'solzen', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), solzen) + ios = NF_INQ_VARID (ncid, 'solazi', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), solazi) ios = NF_INQ_VARID (ncid, 't2m', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), t2m) ios = NF_INQ_VARID (ncid, 'mr2m', varid) @@ -802,6 +813,8 @@ program da_rad_diags ios = NF_INQ_VARID (ncid, 'ret_clw', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), ret_clw) end if + ios = NF_INQ_VARID (ncid, 'cloud_frac', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), cloud_frac) ! ios = NF_CLOSE(ncid) @@ -818,6 +831,8 @@ program da_rad_diags deallocate ( lon ) deallocate ( satzen ) deallocate ( satazi ) + deallocate ( solzen ) + deallocate ( solazi ) deallocate ( t2m ) deallocate ( mr2m ) deallocate ( u10 ) @@ -833,6 +848,7 @@ program da_rad_diags deallocate ( vegfra ) deallocate ( elev ) deallocate ( clwp ) + deallocate ( cloud_frac ) deallocate ( ret_clw ) deallocate ( tb_obs ) deallocate ( tb_bak ) diff --git a/var/da/da_radiance/da_allocate_rad_iv.inc b/var/da/da_radiance/da_allocate_rad_iv.inc index a34bbb0c65..1f8b825637 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 + integer :: n, ix, iy call da_trace_entry("da_allocate_rad_iv") allocate (iv%instid(i)%info%date_char(iv%instid(i)%num_rad)) @@ -106,8 +106,17 @@ subroutine da_allocate_rad_iv (i, nchan, iv) 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)) + allocate (iv%instid(i)%cloud_frac(iv%instid(i)%num_rad)) if ( use_clddet_zz ) then - allocate (iv%instid(i)%cld_qc(iv%instid(i)%num_rad)) + iv%instid(i)%superob_width = 2*ahi_superob_halfwidth+1 + allocate (iv%instid(i)%superob(iv%instid(i)%superob_width, & + iv%instid(i)%superob_width)) + do iy = 1, iv%instid(i)%superob_width + do ix = 1, iv%instid(i)%superob_width + allocate (iv%instid(i)%superob(ix,iy)%cld_qc(iv%instid(i)%num_rad)) + allocate (iv%instid(i)%superob(ix,iy)%tb_obs(nchan,iv%instid(i)%num_rad)) + end do + end do end if if ( use_rttov_kmatrix .or. use_crtm_kmatrix ) then allocate(iv%instid(i)%ts_jacobian(nchan,iv%instid(i)%num_rad)) diff --git a/var/da/da_radiance/da_deallocate_radiance.inc b/var/da/da_radiance/da_deallocate_radiance.inc index 76d8249887..d785b62c1f 100644 --- a/var/da/da_radiance/da_deallocate_radiance.inc +++ b/var/da/da_radiance/da_deallocate_radiance.inc @@ -11,7 +11,7 @@ type (iv_type), intent(inout) :: iv ! Obs. increment structure. type (j_type), intent(inout) :: j ! Cost function. - integer :: i,n,ichan + integer :: i,n,ichan,ix,iy if (trace_use) call da_trace_entry("da_deallocate_radiance") @@ -127,10 +127,17 @@ deallocate (iv%instid(i)%satazi) deallocate (iv%instid(i)%solzen) deallocate (iv%instid(i)%solazi) - deallocate (iv%instid(i)%tropt) + deallocate (iv%instid(i)%tropt) deallocate (iv%instid(i)%gamma_jacobian) + deallocate(iv%instid(i)%cloud_frac) if ( use_clddet_zz ) then - deallocate (iv%instid(i)%cld_qc) + do iy = 1, iv%instid(i)%superob_width + do ix = 1, iv%instid(i)%superob_width + deallocate (iv%instid(i)%superob(ix,iy)%cld_qc) + deallocate (iv%instid(i)%superob(ix,iy)%tb_obs) + end do + end do + deallocate (iv%instid(i)%superob) end if if (ANY(use_satcv)) then if (use_satcv(2)) then diff --git a/var/da/da_radiance/da_initialize_rad_iv.inc b/var/da/da_radiance/da_initialize_rad_iv.inc index 68aa250be6..bb03dd7f02 100644 --- a/var/da/da_radiance/da_initialize_rad_iv.inc +++ b/var/da/da_radiance/da_initialize_rad_iv.inc @@ -11,6 +11,7 @@ subroutine da_initialize_rad_iv (i, n, iv, p) integer, intent(in) :: i, n type(datalink_type), intent(in) :: p type(iv_type), intent(inout) :: iv + integer :: ix, iy call da_trace_entry("da_initialize_rad_iv") @@ -101,16 +102,24 @@ subroutine da_initialize_rad_iv (i, n, iv, p) iv%instid(i)%solzen(n) = p%solzen iv%instid(i)%solazi(n) = p%solazi iv%instid(i)%tropt(n) = 0.0 + iv%instid(i)%cloud_frac(n) = missing_r ! 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 + if ( allocated ( p % superob ) ) then + do iy = 1, iv%instid(i)%superob_width + do ix = 1, iv%instid(i)%superob_width + iv%instid(i)%superob(ix,iy)%tb_obs(:,n) = p % superob(ix,iy) % tb_obs(:,1) + iv%instid(i)%superob(ix,iy)%cld_qc(n)%RTCT = p % superob(ix,iy) % cld_qc(1) % RTCT + iv%instid(i)%superob(ix,iy)%cld_qc(n)%RFMFT = p % superob(ix,iy) % cld_qc(1) % RFMFT + iv%instid(i)%superob(ix,iy)%cld_qc(n)%TEMPIR = p % superob(ix,iy) % cld_qc(1) % TEMPIR + iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_10 = p % superob(ix,iy) % cld_qc(1) % tb_stddev_10 + iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_13 = p % superob(ix,iy) % cld_qc(1) % tb_stddev_13 + iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_14 = p % superob(ix,iy) % cld_qc(1) % tb_stddev_14 + iv%instid(i)%superob(ix,iy)%cld_qc(n)%terr_hgt = p % superob(ix,iy) % cld_qc(1) % terr_hgt + iv%instid(i)%superob(ix,iy)%cld_qc(n)%CIRH2O = p % superob(ix,iy) % cld_qc(1) % CIRH2O + end do + end do + end if end if if ( rtm_option == rtm_option_rttov ) then iv%instid(i)%surftype(n) = 0 diff --git a/var/da/da_radiance/da_qc_ahi.inc b/var/da/da_radiance/da_qc_ahi.inc index 0790c1c6e3..474ea603ab 100644 --- a/var/da/da_radiance/da_qc_ahi.inc +++ b/var/da/da_radiance/da_qc_ahi.inc @@ -53,7 +53,11 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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) + integer*2 :: clddet_tests(iv%instid(i)%superob_width, & + iv%instid(i)%superob_width, & + num_clddet_tests) + integer :: superob_center + integer :: isuper, jsuper real, pointer :: tb_ob(:,:), tb_xb(:,:), tb_inv(:,:), tb_xb_clr(:,:), ca_mean(:,:) integer :: tb_qc(nchan), tb_qc_clddet(nchan) @@ -132,10 +136,14 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) nrej_clddet(:,:)= 0 num_proc_domain = 0 - tb_ob => ob%instid(i)%tb + superob_center = ahi_superob_halfwidth + 1 + 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 + AHIPixelQCLoop: do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2 + + tb_ob => ob%instid(i)%tb + if (iv%instid(i)%info%proc_domain(1,n)) & num_proc_domain = num_proc_domain + 1 @@ -214,12 +222,17 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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)%tropt(n), iv%instid(i)%superob(superob_center,superob_center)%cld_qc(n)%terr_hgt, & iv%instid(i)%info%date_char(n) !JJGDEBUG clddet_tests = 0 + do jsuper = 1, iv%instid(i)%superob_width + do isuper = 1, iv%instid(i)%superob_width + + tb_ob => iv%instid(i)%superob(isuper,jsuper)%tb_obs + 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) / & @@ -258,7 +271,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) tb_qc_clddet = tb_qc - do itest = 1, num_clddet_tests + AHICloudTestLoop: do itest = 1, num_clddet_tests qual_clddet = .true. offset_clddet = 0 crit_clddet = missing_r @@ -267,7 +280,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) !-------------------------------------------------------------------------- ! 4.1 Relative Thermal Contrast Test (RTCT) !-------------------------------------------------------------------------- - crit_clddet = iv%instid(i)%cld_qc(n)%RTCT + crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%RTCT qual_clddet(3:4) = .false. case (2) !-------------------------------------------------------------------------- @@ -284,7 +297,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) ! 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. & + iv%instid(i)%superob(isuper,jsuper)%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 @@ -321,18 +334,18 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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 + crit_clddet = iv%instid(i)%superob(isuper,jsuper)%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 + iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%terr_hgt .le. 2000. & + .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%terr_hgt .ge. 0. & + .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_10 .gt. 0.5 & + .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_14 .gt. 0.5 + crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%CIRH2O case (7) !-------------------------------------------------------------------------- ! 4.7 Modified 4um Emissivity Test (M-4EMISS) @@ -381,7 +394,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) !-------------------------------------------------------------------------- ! 4.10 Temporal Infrared Test (TEMPIR) !-------------------------------------------------------------------------- - crit_clddet = iv%instid(i)%cld_qc(n)%TEMPIR + crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%TEMPIR case default cycle end select @@ -403,9 +416,17 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) ! 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 - clddet_tests(itest) = 1 + clddet_tests(isuper, jsuper, itest) = 1 end if - end do + end do AHICloudTestLoop + end do ! isuper + end do ! jsuper + + if ( iv%instid(i)%superob_width > 1 ) then + iv%instid(i)%cloud_frac(n) = & + real( count(sum(clddet_tests,3) > 0),8) / real( iv%instid(i)%superob_width**2,8) + end if + if (.not. crtm_cloud ) tb_qc = tb_qc_clddet !JJGDEBUG @@ -450,6 +471,9 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) if (iv%instid(i)%cloud_flag(k,n) == qc_bad) tb_qc(k) = qc_bad end do end if ahi_clddet + + tb_ob => ob%instid(i)%tb + ! ---------------------------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 @@ -477,13 +501,17 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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 + if (use_clddet_zz) then + ! SDob cloud inhomogeneous check + do isuper = 1, iv%instid(i)%superob_width + do jsuper = 1, iv%instid(i)%superob_width + if (iv%instid(i)%superob(isuper,jsuper)%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 do + end do end if end if @@ -512,7 +540,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) ngood(k) = ngood(k) + 1 end if end do ! nchan - end do ! end loop pixel + end do AHIPixelQCLoop ! Do inter-processor communication to gather statistics. call da_proc_sum_int (num_proc_domain) diff --git a/var/da/da_radiance/da_radiance.f90 b/var/da/da_radiance/da_radiance.f90 index 587a75348c..e41030262e 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_clddet_zz + use_goesimgobs, pi, earth_radius, satellite_height,use_clddet_zz, ahi_superob_halfwidth #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,clddet_geoir_type + airsr_type,info_type, model_loc_type, varbc_info_type, varbc_type,clddet_geoir_type, superob_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 diff --git a/var/da/da_radiance/da_radiance1.f90 b/var/da/da_radiance/da_radiance1.f90 index 1327930964..441752d34d 100644 --- a/var/da/da_radiance/da_radiance1.f90 +++ b/var/da/da_radiance/da_radiance1.f90 @@ -21,10 +21,11 @@ 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,use_clddet_zz + cv_size_domain_js, calc_weightfunc, deg_to_rad, rad_to_deg,use_clddet_zz, & + ahi_superob_halfwidth use da_define_structures, only : info_type,model_loc_type,maxmin_type, & iv_type, y_type, jo_type,bad_data_type,bad_data_type,number_type, & - be_type, clddet_geoir_type + be_type, clddet_geoir_type, superob_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 @@ -48,6 +49,7 @@ module da_radiance1 type (info_type) :: info type (model_loc_type) :: loc type (clddet_geoir_type) :: cld_qc + type (superob_type), allocatable :: superob(:,:) 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 a4869e3d73..3ad59da7c5 100644 --- a/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc +++ b/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc @@ -77,18 +77,19 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) integer, parameter :: ch13 = 7 integer, parameter :: ch14 = 8 integer, parameter :: ch15 = 9 + integer :: superob_width ! Must be ≥ 0 + integer :: first_boxcenter, last_boxcenter_x, last_boxcenter_y, box_bottom, box_upper, box_left, box_right + integer :: isup, jsup, ixsup, iysup, ibox, jbox, nkeep 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(:) - real(r_kind),allocatable :: data_all(:) - + real(r_kind) :: temp1 = 0.0 character(len=80) tbb_name(10) data tbb_name/'tbb_07','tbb_08','tbb_09','tbb_10','tbb_11', & 'tbb_12','tbb_13','tbb_14','tbb_15','tbb_16'/ @@ -128,7 +129,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) nchan = iv%instid(inst)%nchan write(unit=stdout,fmt=*)'AHI nchan: ',nchan - allocate(data_all(1:nchan)) + superob_width = 2*ahi_superob_halfwidth+1 ! 1.0 Assign file names and prepare to read ahi files !------------------------------------------------------------------------- @@ -319,7 +320,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! close infile_tb file iret = nf90_close(ncid) allocate( tbbp(nlongitude,nlatitude,nchan)) - allocate( ter_sat(nlongitude,nlatitude)) + allocate( ter_sat(nlongitude,nlatitude)) allocate( SDob_10(nlongitude,nlatitude)) allocate( SDob_13(nlongitude,nlatitude)) allocate( SDob_14(nlongitude,nlatitude)) @@ -452,7 +453,11 @@ end if end if ! start scan_loop - scan_loop: do ilatitude=1, nlatitude + first_boxcenter = ahi_superob_halfwidth + 1 + last_boxcenter_y = superob_width * (nlatitude / superob_width) - ahi_superob_halfwidth + last_boxcenter_x = superob_width * (nlongitude / superob_width) - ahi_superob_halfwidth + + scan_loop: do ilatitude=first_boxcenter, nlatitude, superob_width 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 @@ -461,23 +466,24 @@ end if obs_time < time_slots(ifgat) ) exit end do ! start fov_loop - fov_loop: do ilongitude=1, nlongitude - if ( sat_zenith(ilongitude,ilatitude) > 65.0 ) cycle fov_loop - + + jbox = ilatitude/superob_width + 1 + if ( ahi_superob_halfwidth .gt. 0 ) then + box_bottom = superob_width * (jbox-1) +1 + box_upper = superob_width * jbox + else + box_bottom = superob_width * (jbox-1) + box_upper = superob_width * (jbox-1) + end if + + fov_loop: do ilongitude=first_boxcenter, nlongitude, superob_width + + 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 @@ -492,9 +498,9 @@ end if ':', idate5(5), ':', idate5(6) info%elv = 0.0 -! 3.0 Make Thinning -! Map obs to thinning grid -!------------------------------------------------------------------- + ! 3.0 Make Thinning + ! Map obs to thinning grid + !------------------------------------------------------------------- if (thinning) then dlat_earth = info%lat !degree dlon_earth = info%lon @@ -510,17 +516,74 @@ end if end if end if num_ahi_used = num_ahi_used + 1 - data_all = missing_r - - do k=1,nchan - tb = tbb(ilongitude,ilatitude,k) - 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 )) + + ibox = ilongitude/superob_width + 1 + if ( ahi_superob_halfwidth .gt. 0 ) then + box_left = superob_width * (ibox-1) +1 ! will exceed nlatitude/nlongitude if ahi_superob_halfwidth = 0 + box_right = superob_width * ibox + else + box_left = superob_width * (ibox-1) + box_right = superob_width * (ibox-1) + end if + + ! Super-ob BT for this channel + + allocate ( p % tb_inv (1:nchan) ) + p % tb_inv = missing_r + + do k = 1, nchan + nkeep = count(tbb(box_left:box_right,box_bottom:box_upper,k) > 0.0 ) + temp1 = sum ( tbb(box_left:box_right,box_bottom:box_upper,k), & + tbb(box_left:box_right,box_bottom:box_upper,k) > 0.0 ) + + if (ahi_superob_halfwidth .gt.0 .and. nkeep .gt. 0) then + tb = temp1 / real(nkeep,8) + else + ! Extract single pixel BT and radiance value for this channel + tb = tbb(ilongitude,ilatitude,k) + end if + + if (tb > 0.0) p % tb_inv(k) = tb + + end do + + ! Preprocessing for cloud detection including + ! extracting Tb values from cloud QC buffer + + if (use_clddet_zz) then + + if (.not. allocated(p % superob)) then + allocate( p % superob(superob_width,superob_width) ) + end if + + ! Loops over superob pixels + do jsup = 1, superob_width + do isup = 1, superob_width + iysup = ilatitude + jsup-1-ahi_superob_halfwidth + ixsup = ilongitude + isup-1-ahi_superob_halfwidth + if (ixsup < 1 .or. ixsup > nlongitude .or. iysup < 1 .or. iysup > nlatitude ) then + cycle + end if + if ( .not. allocated(p % superob(isup,jsup) % tb_obs ) ) then + allocate ( p % superob(isup,jsup) % tb_obs (1:nchan,1) ) + allocate ( p % superob(isup,jsup) % cld_qc(1) ) + end if + p % superob(isup,jsup) % tb_obs(1:nchan,1) = tbb( ixsup, iysup, 1:nchan ) + p % superob(isup,jsup) % cld_qc(1) % TEMPIR = tempir(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % RTCT = rtct(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % RFMFT = rtmt(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % CIRH2O = cwvt1(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % tb_stddev_10 = SDob_10(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % tb_stddev_13 = SDob_13(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % tb_stddev_14 = SDob_14(ixsup, iysup) + p % superob(isup,jsup) % cld_qc(1) % terr_hgt = ter_sat(ixsup, iysup) + + end do !isup + end do !jsup + + end if + ! 4.0 assign information to sequential radiance structure + !-------------------------------------------------------------------------- p%info = info p%loc = loc p%landsea_mask = 1 @@ -529,10 +592,8 @@ end if 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%cld_qc = cld_qc allocate (p%next) ! add next data p => p%next nullify (p%next) @@ -543,7 +604,6 @@ end if write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_local : ',num_ahi_local_local - deallocate(data_all) ! Deallocate data arrays deallocate(vlatitude) deallocate(vlongitude) deallocate(tbb) @@ -629,6 +689,15 @@ end if prev => p end if deallocate ( current % tb_inv ) + if ( allocated ( current % superob ) ) then + do jsup = 1, superob_width + do isup = 1, superob_width + deallocate ( current % superob(isup,jsup) % tb_obs ) + deallocate ( current % superob(isup,jsup) % cld_qc ) + end do + end do + deallocate ( current % superob ) + end if deallocate ( current ) num_ahi_thinned = num_ahi_thinned + 1 num_ahi_used = num_ahi_used - 1 @@ -697,6 +766,15 @@ end if p => p%next ! free current data deallocate ( current % tb_inv ) + if ( allocated ( current % superob ) ) then + do jsup = 1, superob_width + do isup = 1, superob_width + deallocate ( current % superob(isup,jsup) % tb_obs ) + deallocate ( current % superob(isup,jsup) % cld_qc ) + end do + end do + deallocate ( current % superob ) + end if deallocate ( current ) end do deallocate ( p ) diff --git a/var/da/da_radiance/da_write_iv_rad_ascii.inc b/var/da/da_radiance/da_write_iv_rad_ascii.inc index 4c14d51466..6e8a6acfeb 100644 --- a/var/da/da_radiance/da_write_iv_rad_ascii.inc +++ b/var/da/da_radiance/da_write_iv_rad_ascii.inc @@ -64,18 +64,18 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) ' index-of-channels : ' write(unit=innov_rad_unit,fmt='(10i5)') iv%instid(i)%ichan if ( amsr2 ) then - write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi clw' + write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi solzen solazi clw' else - write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi' + write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi solzen solazi' end if write(unit=innov_rad_unit,fmt='(a)') ' grid%xb-surf-info : i t2m mr2m(ppmv) u10 v10 ps ts smois tslb snowh isflg & - & soiltyp vegtyp vegfra elev clwp' + & soiltyp vegtyp vegfra elev clwp cloud_frac' ndomain = 0 do n =1,iv%instid(i)%num_rad if (iv%instid(i)%info%proc_domain(1,n)) then ndomain=ndomain+1 if ( amsr2 ) then ! write out clw - write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,4f8.2,f8.3)') 'INFO : ', ndomain, & + write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2,f8.3)') 'INFO : ', ndomain, & iv%instid(i)%info%date_char(n), & iv%instid(i)%scanpos(n), & iv%instid(i)%landsea_mask(n), & @@ -84,9 +84,11 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) 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)%clw(n) else ! no clw info - write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,4f8.2)') 'INFO : ', ndomain, & + write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2)') 'INFO : ', ndomain, & iv%instid(i)%info%date_char(n), & iv%instid(i)%scanpos(n), & iv%instid(i)%landsea_mask(n), & @@ -94,7 +96,9 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) 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)%satazi(n), & + iv%instid(i)%solzen(n), & + iv%instid(i)%solazi(n) end if select case (iv%instid(i)%isflg(n)) case (0) ; @@ -114,7 +118,7 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) case (7) ; surftype = 'MSNO : ' end select - write(unit=innov_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3)') surftype, n, & + write(unit=innov_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3,f15.5)') surftype, n, & iv%instid(i)%t2m(n), & iv%instid(i)%mr2m(n), & iv%instid(i)%u10(n), & @@ -129,7 +133,8 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) nint(iv%instid(i)%vegtyp(n)), & iv%instid(i)%vegfra(n), & iv%instid(i)%elevation(n), & - iv%instid(i)%clwp(n) + iv%instid(i)%clwp(n), & + iv%instid(i)%cloud_frac(n) write(unit=innov_rad_unit,fmt='(a)') 'OBS : ' write(unit=innov_rad_unit,fmt='(10f11.2)') ob%instid(i)%tb(:,n) diff --git a/var/da/da_radiance/da_write_oa_rad_ascii.inc b/var/da/da_radiance/da_write_oa_rad_ascii.inc index 8f6d5f1bfc..2f058839df 100644 --- a/var/da/da_radiance/da_write_oa_rad_ascii.inc +++ b/var/da/da_radiance/da_write_oa_rad_ascii.inc @@ -55,9 +55,9 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) ' index-of-channels : ' write(unit=oma_rad_unit,fmt='(10i5)') iv%instid(i)%ichan if ( amsr2 ) then - write(unit=oma_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi clw' + write(unit=oma_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi solzen solazi clw' else - write(unit=oma_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi ' + write(unit=oma_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask elv lat lon satzen satazi solzen solazi' end if write(unit=oma_rad_unit,fmt='(a)') ' xb-surf-info : i t2m mr2m(ppmv) u10 v10 ps ts smois tslb snowh isflg & & soiltyp vegtyp vegfra elev clwp' @@ -66,7 +66,7 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) if (iv%instid(i)%info%proc_domain(1,n)) then ndomain=ndomain+1 if ( amsr2 ) then !write out clw - write(unit=oma_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,4f8.2,f8.3)') 'INFO : ', ndomain, & + write(unit=oma_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2,f8.3)') 'INFO : ', ndomain, & iv%instid(i)%info%date_char(n), & iv%instid(i)%scanpos(n), & iv%instid(i)%landsea_mask(n), & @@ -75,9 +75,11 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) 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)%clw(n) else !no clw info - write(unit=oma_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,4f8.2)') 'INFO : ', ndomain, & + write(unit=oma_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2)') 'INFO : ', ndomain, & iv%instid(i)%info%date_char(n), & iv%instid(i)%scanpos(n), & iv%instid(i)%landsea_mask(n), & @@ -85,7 +87,9 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) 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)%satazi(n), & + iv%instid(i)%solzen(n), & + iv%instid(i)%solazi(n) end if select case (iv%instid(i)%isflg(n)) case (0) ; @@ -105,7 +109,7 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) case (7) ; surftype = 'MSNO : ' end select - write(unit=oma_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3)') surftype, n, & + write(unit=oma_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3,f15.5)') surftype, n, & iv%instid(i)%t2m(n), & iv%instid(i)%mr2m(n), & iv%instid(i)%u10(n), & @@ -120,7 +124,8 @@ subroutine da_write_oa_rad_ascii (it, ob, iv, re ) nint(iv%instid(i)%vegtyp(n)), & iv%instid(i)%vegfra(n), & iv%instid(i)%elevation(n), & - iv%instid(i)%clwp(n) + iv%instid(i)%clwp(n), & + iv%instid(i)%cloud_frac(n) write(unit=oma_rad_unit,fmt='(a)') 'OBS : ' write(unit=oma_rad_unit,fmt='(10f11.2)') ob%instid(i)%tb(:,n)