From 05a00cf22e3857bf5c613539f6b8f03722d80c37 Mon Sep 17 00:00:00 2001 From: Craig Schwartz Date: Wed, 7 Dec 2022 19:42:36 -0700 Subject: [PATCH] Enhancements for AHI radiance DA (#1774) TYPE: Enhancement KEYWORDS: WRFDA, AHI, Himawari-8 SOURCE: Craig Schwartz (NCAR/MMM) DESCRIPTION OF CHANGES: This PR makes several enhancements for assimilating Himawari-8 radiance data, including (1) Introduction of an all-sky obs error model (Harnisch et al., 2016) for all-sky AHI DA; (2) Optional read and use of AHI level-2 product (e.g., cloud mask); (3) More efficient read of a sub-area of the full disk data; (4) Allow the use of offline statistics of constant bias correction values; (5) More diagnostic output (peak of weighting functions, cloud ice water path, cloud flag) in omb/oma files. LIST OF MODIFIED FILES: var/da/da_radiance/module_radiance.f90 Registry/registry.var var/da/da_define_structures/da_define_structures.f90 var/da/da_monitor/da_rad_diags.f90 var/da/da_radiance/da_allocate_rad_iv.inc var/da/da_radiance/da_deallocate_radiance.inc var/da/da_radiance/da_get_innov_vector_crtm.inc var/da/da_radiance/da_initialize_rad_iv.inc var/da/da_radiance/da_qc_ahi.inc var/da/da_radiance/da_radiance.f90 var/da/da_radiance/da_radiance1.f90 var/da/da_radiance/da_radiance_init.inc var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc var/da/da_radiance/da_write_iv_rad_ascii.inc var/da/da_radiance/module_radiance.f90 var/run/ahi_info var/run/radiance_info/himawari-8-ahi.info TESTING: (1) WRFDA regression test passed; (2) AHI all-sky DA runs also Ok. Release Note: Enhancements for AHI radiance DA, including all-sky observation error model, Level-2 AHI product read, and more diagnostic output. Xu, D. M., Z. Q. Liu, S. Y. Fan, M. Chen, and F. F. Shen, 2021: Assimilating all-sky infrared radiances from Himawari-8 using the 3DVar method for the prediction of a severe storm over North China. Adv. Atmos. Sci., 38(4), 661-676. --- Registry/registry.var | 2 + .../da_define_structures.f90 | 1 + var/da/da_monitor/da_rad_diags.f90 | 45 ++++++++- var/da/da_radiance/da_allocate_rad_iv.inc | 1 + var/da/da_radiance/da_deallocate_radiance.inc | 13 +++ .../da_radiance/da_get_innov_vector_crtm.inc | 22 ++++- var/da/da_radiance/da_initialize_rad_iv.inc | 1 + var/da/da_radiance/da_qc_ahi.inc | 49 ++++++++-- var/da/da_radiance/da_radiance.f90 | 2 +- var/da/da_radiance/da_radiance1.f90 | 2 +- var/da/da_radiance/da_radiance_init.inc | 34 +++++++ .../da_read_obs_netcdf4ahi_jaxa.inc | 91 ++++++++++++++++--- var/da/da_radiance/da_write_iv_rad_ascii.inc | 56 +++++++++++- var/da/da_radiance/module_radiance.f90 | 4 + var/run/ahi_info | 2 +- var/run/radiance_info/himawari-8-ahi.info | 22 ++--- 16 files changed, 297 insertions(+), 50 deletions(-) diff --git a/Registry/registry.var b/Registry/registry.var index d8937a9581..e3c6c9cfa3 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -467,6 +467,8 @@ rconfig character crtm_irwater_coef namelist,wrfvar14 1 "Nalli.IRwater rconfig character crtm_mwwater_coef namelist,wrfvar14 1 "FASTEM5.MWwater.EmisCoeff.bin" - "crtm_mwwater_coef" "" "" rconfig character crtm_irland_coef namelist,wrfvar14 1 "USGS.IRland.EmisCoeff.bin" - "crtm_irland_coef" "" "" rconfig character crtm_visland_coef namelist,wrfvar14 1 "USGS.VISland.EmisCoeff.bin" - "crtm_visland_coef" "" "" +rconfig logical ahi_use_symm_obs_err namelist,wrfvar14 1 .false. - "ahi_use_symm_obs_err" "" "" +rconfig logical ahi_apply_clrsky_bias namelist,wrfvar14 1 .false. - "ahi_apply_clrsky_bias" "" "" rconfig integer num_pseudo namelist,wrfvar15 1 0 - "num_pseudo" "" "" rconfig real pseudo_x namelist,wrfvar15 1 1.0 - "pseudo_x" "" "" rconfig real pseudo_y namelist,wrfvar15 1 1.0 - "pseudo_y" "" "" diff --git a/var/da/da_define_structures/da_define_structures.f90 b/var/da/da_define_structures/da_define_structures.f90 index a703ff0aff..7d3249e4c0 100644 --- a/var/da/da_define_structures/da_define_structures.f90 +++ b/var/da/da_define_structures/da_define_structures.f90 @@ -636,6 +636,7 @@ module da_define_structures real, pointer :: vegtyp(:) real, pointer :: vegfra(:) real, pointer :: clwp(:) ! model/guess clwp + real, pointer :: cip(:) ! model/guess cloud-ice path real, pointer :: clw(:) ! currently AMSR2 only real, pointer :: ps_jacobian(:,:) ! only RTTOV real, pointer :: ts_jacobian(:,:) ! only over water CRTM diff --git a/var/da/da_monitor/da_rad_diags.f90 b/var/da/da_monitor/da_rad_diags.f90 index 15a4b8d997..af42a488ff 100644 --- a/var/da/da_monitor/da_rad_diags.f90 +++ b/var/da/da_monitor/da_rad_diags.f90 @@ -60,8 +60,12 @@ program da_rad_diags real*4, dimension(:), allocatable :: ret_clw 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 + real*4, dimension(:), allocatable :: cip ! cloud-ice path + integer, dimension(:), allocatable :: cloudflag ! cloudflag from L2 AHI 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 :: tb_bak_clr ! clear-sky brightness temp + real*4, dimension(:,:), allocatable :: weightfunc_peak ! peak of weighting function real*4, dimension(:,:), allocatable :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water real*4, dimension(:,:), allocatable :: prf_ice, prf_rain, prf_snow, prf_grau, prf_hail real*4, dimension(:,:), allocatable :: prf_water_reff, prf_ice_reff, prf_rain_reff @@ -248,9 +252,13 @@ program da_rad_diags allocate ( vegfra(1:total_npixel) ) allocate ( elev(1:total_npixel) ) allocate ( clwp(1:total_npixel) ) !model/guess clwp + allocate ( cip(1:total_npixel) ) ! model/guess cloud-ice path + allocate ( cloudflag(1:total_npixel) ) ! AHI Level-2 cloud flag allocate ( cloud_frac(1:total_npixel) ) allocate ( tb_obs(1:nchan,1:total_npixel) ) allocate ( tb_bak(1:nchan,1:total_npixel) ) + allocate ( tb_bak_clr(1:nchan,1:total_npixel) ) + allocate ( weightfunc_peak(1:nchan,1:total_npixel) ) allocate ( tb_inv(1:nchan,1:total_npixel) ) allocate ( tb_oma(1:nchan,1:total_npixel) ) allocate ( tb_err(1:nchan,1:total_npixel) ) @@ -320,11 +328,14 @@ program da_rad_diags ! initialize tb_obs = missing_r tb_bak = missing_r + tb_bak_clr = missing_r + weightfunc_peak = missing_r tb_inv = missing_r tb_oma = missing_r tb_err = missing_r ncname = 'diags_'//trim(instid(iinst))//"_"//datestr1(itime)//'.nc' - ios = NF_CREATE(trim(ncname), NF_CLOBBER, ncid) ! NF_CLOBBER specifies the default behavior of + ios = NF_CREATE(trim(ncname), NF_NETCDF4, ncid) ! Change to output netcdf4 files + !ios = NF_CREATE(trim(ncname), NF_CLOBBER, ncid) ! NF_CLOBBER specifies the default behavior of ! overwritting any existing dataset with the ! same file name if ( ios /= 0 ) then @@ -349,15 +360,23 @@ program da_rad_diags 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,f15.5)',iostat=ios) & + read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2,f15.5,f15.5,i7)',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), cloud_frac(ipixel) + vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel), cloud_frac(ipixel), cip(ipixel), cloudflag(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 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_bak(:,ipixel) - read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! IVBC + read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! IVBC or BAK_clr or WEIGHTFUNC_PEAK + if ( buf(1:7) == "BAK_clr" ) then + read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_bak_clr(:,ipixel) + read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! IVBC or WEIGHTFUNC_PEAK + end if + if ( buf(1:15) == "WEIGHTFUNC_PEAK" ) then + read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) weightfunc_peak(:,ipixel) + read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! IVBC + end if read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_inv(:,ipixel) read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OMA or EMS if ( buf(1:3) == "OMA" ) then @@ -490,6 +509,10 @@ program da_rad_diags ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r) ios = NF_DEF_VAR(ncid, 'tb_bak', NF_FLOAT, 2, ishape(1:2), varid) ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r) + ios = NF_DEF_VAR(ncid, 'tb_bak_clr', NF_FLOAT, 2, ishape(1:2), varid) + ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r) + ios = NF_DEF_VAR(ncid, 'weightfunc_peak', NF_FLOAT, 2, ishape(1:2), varid) + ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r) ios = NF_DEF_VAR(ncid, 'tb_inv', NF_FLOAT, 2, ishape(1:2), varid) ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r) ios = NF_DEF_VAR(ncid, 'tb_oma', NF_FLOAT, 2, ishape(1:2), varid) @@ -601,6 +624,8 @@ program da_rad_diags ios = NF_DEF_VAR(ncid, 'vegfra', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'elev', NF_FLOAT, 1, ishape(1), varid) ios = NF_DEF_VAR(ncid, 'clwp', NF_FLOAT, 1, ishape(1), varid) + ios = NF_DEF_VAR(ncid, 'cip', NF_FLOAT, 1, ishape(1), varid) + ios = NF_DEF_VAR(ncid, 'cloudflag', NF_INT, 1, ishape(1), varid) if ( amsr2 ) then ios = NF_DEF_VAR(ncid, 'ret_clw', NF_FLOAT, 1, ishape(1), varid) end if @@ -626,6 +651,10 @@ program da_rad_diags ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_obs) ios = NF_INQ_VARID (ncid, 'tb_bak', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_bak) + ios = NF_INQ_VARID (ncid, 'tb_bak_clr', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_bak_clr) + ios = NF_INQ_VARID (ncid, 'weightfunc_peak', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), weightfunc_peak) ios = NF_INQ_VARID (ncid, 'tb_inv', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_inv) ios = NF_INQ_VARID (ncid, 'tb_oma', varid) @@ -809,6 +838,10 @@ program da_rad_diags ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), elev) ios = NF_INQ_VARID (ncid, 'clwp', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), clwp) + ios = NF_INQ_VARID (ncid, 'cip', varid) + ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), cip) + ios = NF_INQ_VARID (ncid, 'cloudflag', varid) + ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), cloudflag) if ( amsr2 ) then ios = NF_INQ_VARID (ncid, 'ret_clw', varid) ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), ret_clw) @@ -848,10 +881,14 @@ program da_rad_diags deallocate ( vegfra ) deallocate ( elev ) deallocate ( clwp ) + deallocate ( cip ) + deallocate ( cloudflag ) deallocate ( cloud_frac ) deallocate ( ret_clw ) deallocate ( tb_obs ) deallocate ( tb_bak ) + deallocate ( tb_bak_clr ) + deallocate ( weightfunc_peak ) deallocate ( tb_inv ) deallocate ( tb_oma ) deallocate ( ems ) diff --git a/var/da/da_radiance/da_allocate_rad_iv.inc b/var/da/da_radiance/da_allocate_rad_iv.inc index 062422cb38..d5b5eb61ad 100644 --- a/var/da/da_radiance/da_allocate_rad_iv.inc +++ b/var/da/da_radiance/da_allocate_rad_iv.inc @@ -88,6 +88,7 @@ subroutine da_allocate_rad_iv (i, nchan, iv) allocate (iv%instid(i)%tb_xb(nchan,iv%instid(i)%num_rad)) if ( crtm_cloud ) then allocate (iv%instid(i)%tb_xb_clr(nchan,iv%instid(i)%num_rad)) + allocate (iv%instid(i)%cip(iv%instid(i)%num_rad)) end if allocate (iv%instid(i)%tb_qc(nchan,iv%instid(i)%num_rad)) allocate (iv%instid(i)%tb_inv(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 602ada8da3..e0e9f71b55 100644 --- a/var/da/da_radiance/da_deallocate_radiance.inc +++ b/var/da/da_radiance/da_deallocate_radiance.inc @@ -30,6 +30,18 @@ deallocate (satinfo(i) % bcoef0) deallocate (satinfo(i) % error_std) + ! Deallocate extra variables for AHI + if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then + deallocate ( satinfo(i) % BTLim) + deallocate ( satinfo(i) % ca1) + deallocate ( satinfo(i) % ca2) + deallocate ( satinfo(i) % clearSkyBias) + endif + + if (use_error_factor_rad) then + deallocate (satinfo(i) % error_factor) + endif + deallocate (ob%instid(i) % ichan) deallocate (iv%instid(i) % ichan) @@ -111,6 +123,7 @@ deallocate (iv%instid(i)%tb_xb) if ( crtm_cloud ) then deallocate (iv%instid(i)%tb_xb_clr) + deallocate (iv%instid(i)%cip) end if deallocate (iv%instid(i)%tb_qc) deallocate (iv%instid(i)%tb_inv) 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 31b18b2ef0..d41260953d 100644 --- a/var/da/da_radiance/da_get_innov_vector_crtm.inc +++ b/var/da/da_radiance/da_get_innov_vector_crtm.inc @@ -41,6 +41,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) integer :: model_isflg real :: model_qcw(kms:kme) + real :: model_qci(kms:kme) real :: model_rho(kms:kme) real, allocatable :: model_snow(:) ! snow water equivalent, different from model_snowh, ! used in calculating reff_water @@ -55,9 +56,9 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) integer :: inst, nchanl, n1,n2 integer :: ipred, npred, gammapred - ! variables for computing clwp - real :: clw(kms:kme), dpf(kms:kme) - real :: clwp + ! variables for computing clwp (clw/clwp) and cloud ice (ci/cip) + real :: clw(kms:kme), dpf(kms:kme), ci(kms:kme) + real :: clwp, cip real :: total_od @@ -241,7 +242,8 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) calc_tb_clr = .false. if ( crtm_cloud .and. & - trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' ) then + ( trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' .or. & + trim( crtm_sensor_name(rtminit_sensor(inst))) == 'ahi') ) then !Tb_clear_sky is only needed for symmetric obs error model !symmetric obs error model only implemented for amsr2 for now calc_tb_clr = .true. @@ -389,6 +391,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) if (crtm_cloud) then call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,qci) + call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,model_qci(kte-k+1:kte-k+1)) call da_interp_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n,qrn) @@ -511,6 +514,16 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) clwp = clwp + clw(k) end do + ! computing ice path from guess + if (crtm_cloud) then + cip = 0.0 + do k = kts,kte ! from top to bottom + dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1)) + ci(k) = model_qci(k)*dpf(k)/gravity ! kg/m2 or mm + cip = cip + ci(k) + end do + endif + ! CRTM GeometryInfo Structure GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n) GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n) @@ -853,6 +866,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv ) iv%instid(inst)%vegtyp(n) = model_ivgtyp iv%instid(inst)%vegfra(n) = model_vegfra(n) iv%instid(inst)%clwp(n) = clwp + if (crtm_cloud) iv%instid(inst)%cip(n) = cip iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage iv%instid(inst)%land_coverage(n) = Surface(1)%land_coverage iv%instid(inst)%ice_coverage(n) = Surface(1)%ice_coverage diff --git a/var/da/da_radiance/da_initialize_rad_iv.inc b/var/da/da_radiance/da_initialize_rad_iv.inc index 8caa0ee00a..8c6de31102 100644 --- a/var/da/da_radiance/da_initialize_rad_iv.inc +++ b/var/da/da_radiance/da_initialize_rad_iv.inc @@ -84,6 +84,7 @@ subroutine da_initialize_rad_iv (i, n, iv, p) iv%instid(i)%tb_xb(:,n) = 0.0 if ( crtm_cloud ) then iv%instid(i)%tb_xb_clr(:,n) = 0.0 + iv%instid(i)%cip(n) = 0.0 end if iv%instid(i)%tb_inv(:,n) = p%tb_inv(:) iv%instid(i)%tb_qc(:,n) = 0 diff --git a/var/da/da_radiance/da_qc_ahi.inc b/var/da/da_radiance/da_qc_ahi.inc index 474ea603ab..32108a02c4 100644 --- a/var/da/da_radiance/da_qc_ahi.inc +++ b/var/da/da_radiance/da_qc_ahi.inc @@ -58,6 +58,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) num_clddet_tests) integer :: superob_center integer :: isuper, jsuper + real :: cm, co, ca ! variables for all-sky obs error real, pointer :: tb_ob(:,:), tb_xb(:,:), tb_inv(:,:), tb_xb_clr(:,:), ca_mean(:,:) integer :: tb_qc(nchan), tb_qc_clddet(nchan) @@ -215,7 +216,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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) + iv%instid(i)%tb_xb_clr(:,n) end if if (print_cld_debug) write(stdout,'(A,I8,8F12.4,2x,A)') 'PIXEL_DEBUG5: ', n, & @@ -410,7 +411,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) qual_clddet ) if (reject_clddet) then - tb_qc_clddet = qc_bad + tb_qc_clddet = qc_bad ! CSS do we want to set it bad for a given pixel within a superob? 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) @@ -422,12 +423,10 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) 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 + iv%instid(i)%cloud_frac(n) = & + real( count(sum(clddet_tests,3) > 0),8) / real( iv%instid(i)%superob_width**2,8) - if (.not. crtm_cloud ) tb_qc = tb_qc_clddet + if (.not. crtm_cloud ) tb_qc = tb_qc_clddet ! CSS logic here isn't quite right if superobbing !JJGDEBUG if (print_cld_debug) write(stdout,'(A,I8,*(2x,I1:))') 'PIXEL_DEBUG6: ', n, clddet_tests @@ -476,7 +475,7 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) ! ---------------------------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 + ! tb_xb_clr => iv%instid(i)%tb_xb_clr ! currently not used if (.not. crtm_cloud ) then do k = 1, nchan if (use_error_factor_rad) then @@ -486,6 +485,35 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k) end if end do ! nchan + + else ! Added this else block...until now obs error = 500 if crtm_cloud = T...not good! + + if ( ahi_use_symm_obs_err ) then + do k = 1, nchan + ! Okamato et al. (2014) + !cm = iv%instid(i)%tb_xb(k,n) - iv%instid(i)%tb_xb_clr(k,n) + !co = ob%instid(i)%tb(k,n) - iv%instid(i)%tb_xb_clr(k,n) + + ! Harnisch et al. (2016) + cm = max(0.0, satinfo(i)%BTLim(k) - iv%instid(i)%tb_xb(k,n)) + co = max(0.0, satinfo(i)%BTLim(k) - ob%instid(i)%tb(k,n)) + + ! Symmetric cloud amount + ca = 0.5*( abs(cm) + abs(co) ) + + ! Figure out observation error as a function of ca + if (ca.lt.satinfo(i)%ca1(k)) then + iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k) + else if (ca.ge.satinfo(i)%ca1(k) .and. ca.lt.satinfo(i)%ca2(k)) then + iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)+ & + (ca-satinfo(i)%ca1(k))*(satinfo(i)%error_cld(k)-satinfo(i)%error_std(k))/(satinfo(i)%ca2(k)-satinfo(i)%ca1(k)) + else + iv%instid(i)%tb_error(k,n)= satinfo(i)%error_cld(k) + end if + end do ! nchan + else + iv%instid(i)%tb_error(:,n)= 500.0 ! this is the default + endif end if ! 6.0 gross 8k check -clr,sdobs-clr @@ -528,7 +556,10 @@ subroutine da_qc_ahi (it, i, nchan, ob, iv) end do ! nchan !final QC decsion - if (crtm_cloud ) tb_qc(2:4) = qc_good ! no qc for crtm_cloud + ! CSS comment out below. Dangerous, especially if satinfo(i)%iuse(k) = -1 + ! CSS also okay to fail 3std check if using symmetric error model for obs errors + !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 diff --git a/var/da/da_radiance/da_radiance.f90 b/var/da/da_radiance/da_radiance.f90 index 99883ba663..167d0480b5 100644 --- a/var/da/da_radiance/da_radiance.f90 +++ b/var/da/da_radiance/da_radiance.f90 @@ -58,7 +58,7 @@ 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, ahi_superob_halfwidth + use_goesimgobs, pi, earth_radius, satellite_height,use_clddet_zz, ahi_superob_halfwidth, ahi_apply_clrsky_bias #ifdef CRTM use da_crtm, only : da_crtm_init, da_get_innov_vector_crtm diff --git a/var/da/da_radiance/da_radiance1.f90 b/var/da/da_radiance/da_radiance1.f90 index 9ee1a791e1..e4690c086b 100644 --- a/var/da/da_radiance/da_radiance1.f90 +++ b/var/da/da_radiance/da_radiance1.f90 @@ -22,7 +22,7 @@ module da_radiance1 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, & - ahi_superob_halfwidth + ahi_superob_halfwidth, ahi_use_symm_obs_err 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, superob_type diff --git a/var/da/da_radiance/da_radiance_init.inc b/var/da/da_radiance/da_radiance_init.inc index ea920a8182..3773b40122 100644 --- a/var/da/da_radiance/da_radiance_init.inc +++ b/var/da/da_radiance/da_radiance_init.inc @@ -34,6 +34,7 @@ subroutine da_radiance_init(iv,ob) integer :: iunit character(len=filename_len) :: filename character(len=20) :: cdum + character(len=12) :: cdum12 real :: error_cld ! local variables for tuning error factor @@ -195,6 +196,14 @@ subroutine da_radiance_init(iv,ob) satinfo(n) % error_cld(:) = 500.0 !initialize + ! Allocate additional fields for AHI + if ( index(iv%instid(n)%rttovid_string, 'ahi') > 0 ) then + allocate ( satinfo(n) % BTLim(nchanl(n)) ) + allocate ( satinfo(n) % ca1(nchanl(n)) ) + allocate ( satinfo(n) % ca2(nchanl(n)) ) + allocate ( satinfo(n) % clearSkyBias(nchanl(n)) ) + endif + read(iunit,*) do j = 1, nchanl(n) read(iunit,'(1x,5i5,2e18.10,a20)') & @@ -224,6 +233,31 @@ subroutine da_radiance_init(iv,ob) satinfo(n)%error_cld(j) = error_cld end if end if + + ! If AHI, read some extra things + ! Unfortunately, we need to read everything again... + if ( index(iv%instid(n)%rttovid_string, 'ahi') > 0 ) then + backspace(iunit) + read(iunit,'(1x,5i5,2e18.10,a12,f8.2,2f6.2,f9.3)') & + wmo_sensor_id, & + satinfo(n)%ichan(j), & + sensor_type, & + satinfo(n)%iuse(j) , & + idum, & + satinfo(n)%error(j), & + satinfo(n)%polar(j), & + cdum12, & + satinfo(n)%BTLim(j), & + satinfo(n)%ca1(j) , & + satinfo(n)%ca2(j) , & + satinfo(n)%clearSkyBias(j) + if ( j == 1 ) then + write(*,*)'Reading extra data for AHI' + write(*,*)'Channel BTLim ca1 ca2 clearSkyBias' + endif + write(*,fmt='(i7,6x,4f9.3)') satinfo(n)%ichan(j), satinfo(n)%BTLim(j), satinfo(n)%ca1(j), satinfo(n)%ca2(j), satinfo(n)%clearSkyBias(j) + endif + iv%instid(n)%ichan(j) = satinfo(n)%ichan(j) ob%instid(n)%ichan(j) = satinfo(n)%ichan(j) end do 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 3ad59da7c5..73a7e9a9b3 100644 --- a/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc +++ b/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc @@ -65,6 +65,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) integer(i_kind) :: itx, itt character(80) :: filename1,filename2 integer :: nchan,nlongitude,nlatitude,ilongitude,ilatitude,ichannels + integer :: lonstart,latstart,ahi_info_unit integer :: LatDimID,LonDimID integer :: latid,lonid,tbb_id,sazid,cltyid,ter_id integer :: nfile @@ -85,6 +86,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) 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(:,:) + real(r_kind), allocatable :: cloud_mask(:,:) ! For reading L2 cloud information character(200) :: filenameStatic ! Allocatable arrays integer(i_kind),allocatable :: ptotal(:) @@ -112,6 +114,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) num_ahi_used = 0 num_ahi_thinned = 0 + inst = 0 do i = 1, rtminit_nsensor if (platform_id == rtminit_platform(i) & .and. satellite_id == rtminit_satid(i) & @@ -135,6 +138,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) !------------------------------------------------------------------------- nfile = 0 !initialize fname_tb(:) = '' !initialize + fname_tbp(:) = '' ! initialize ! first check if ahi nc file is available filename1 = trim(infile_tb) @@ -169,6 +173,24 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) return end if + ! Added to get area information from the file, in case you don't + ! want to read the entire area, which is large and memory intensive + !open the data area info file + call da_get_unit(ahi_info_unit) + open(unit=ahi_info_unit,file='ahi_info',status='old',iostat=iret) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__,(/"area_info file read error"/)) + endif + !read date information + read(ahi_info_unit,*) + read(ahi_info_unit,*) + read(ahi_info_unit,*) + read(ahi_info_unit,*) + read(ahi_info_unit,*) + read(ahi_info_unit,*) lonstart,latstart,nlongitude,nlatitude + close(ahi_info_unit) + call da_free_unit(ahi_info_unit) + infile_loop: do ifile = 1, nfile num_ahi_file_local = 0 num_ahi_local_local = 0 @@ -182,14 +204,22 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) cycle infile_loop endif - ! read dimensions: latitude and longitude - iret = nf90_inq_dimid(ncid, "latitude", LatDimID) - iret = nf90_inquire_dimension(ncid, LatDimID, len=nlatitude) - - iret = nf90_inq_dimid(ncid, "longitude", LonDimID) - iret = nf90_inquire_dimension(ncid, LonDimID, len=nlongitude) + ! read dimensions: latitude and longitude + ! If any value in the file is < 0, that means to ignore them + ! and get the number of lats/lons from file, and use everything + if ( nlongitude < 0 .or. nlatitude < 0 .or. lonstart < 0 .or. latstart < 0 ) then + iret = nf90_inq_dimid(ncid, "latitude", LatDimID) + iret = nf90_inquire_dimension(ncid, LatDimID, len=nlatitude) + + iret = nf90_inq_dimid(ncid, "longitude", LonDimID) + iret = nf90_inquire_dimension(ncid, LonDimID, len=nlongitude) + + lonstart = 1 + latstart = 1 + endif write(unit=stdout,fmt=*)'nlongitude,nlatitude: ',nlongitude,nlatitude + write(unit=stdout,fmt=*)'lonstart,latstart: ',lonstart,latstart ! read array: time iret = nf90_get_att(ncid, nf90_global, "id", filename) @@ -206,7 +236,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read lat iret = nf90_inq_varid(ncid, 'latitude', latid) allocate( vlatitude(nlatitude)) - iret = nf90_get_var(ncid, latid, vlatitude) + iret = nf90_get_var(ncid, latid, vlatitude,start=(/latstart/),count=(/nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: Latitude of Observation Point"/)) @@ -214,7 +244,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read lon iret = nf90_inq_varid(ncid, 'longitude', lonid) allocate( vlongitude(nlongitude)) - iret = nf90_get_var(ncid, lonid, vlongitude) + iret = nf90_get_var(ncid, lonid, vlongitude,start=(/lonstart/),count=(/nlongitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: Longitude of Observation Point"/)) @@ -228,7 +258,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) allocate( tbb(nlongitude,nlatitude,nchan)) do k=1,nchan iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id) - iret = nf90_get_var(ncid, tbb_id, tbb(:,:,k)) + iret = nf90_get_var(ncid, tbb_id, tbb(:,:,k),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: Brightness Temperature"/)) @@ -248,7 +278,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read iret = nf90_inq_varid(ncid, 'SAZ', sazid) allocate( sat_zenith(nlongitude,nlatitude)) - iret = nf90_get_var(ncid, sazid, sat_zenith) + iret = nf90_get_var(ncid, sazid, sat_zenith,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: satellite zenith angle"/)) @@ -267,7 +297,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read iret = nf90_inq_varid(ncid, 'SAA', sazid) allocate( sat_azi(nlongitude,nlatitude)) - iret = nf90_get_var(ncid, sazid, sat_azi) + iret = nf90_get_var(ncid, sazid, sat_azi,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: satellite azimuth angle"/)) @@ -285,7 +315,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read iret = nf90_inq_varid(ncid, 'SOZ', sazid) allocate( sol_zenith(nlongitude,nlatitude)) - iret = nf90_get_var(ncid, sazid, sol_zenith) + iret = nf90_get_var(ncid, sazid, sol_zenith,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: satellite zenith angle"/)) @@ -304,7 +334,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) ! read iret = nf90_inq_varid(ncid, 'SOA', sazid) allocate( sol_azi(nlongitude,nlatitude)) - iret = nf90_get_var(ncid, sazid, sol_azi) + iret = nf90_get_var(ncid, sazid, sol_azi,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: satellite azimuth angle"/)) @@ -319,6 +349,28 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) 'solar azimuth angle(pixel=1,scan=1): ',sol_azi(1,1) ! close infile_tb file iret = nf90_close(ncid) + + ! Read the level-2 cloud information file from JAXA + ! This must be on the same grid as everything else + allocate(cloud_mask(nlongitude,nlatitude)) + cloud_mask = -9999.0 + inquire(file='L2AHICLP', exist=fexist) + if ( fexist ) then + iret = nf90_open('L2AHICLP', nf90_NOWRITE, ncid) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/)) + endif + + iret = nf90_inq_varid(ncid, "CLTYPE", cltyid) + iret = nf90_get_var(ncid,cltyid,cloud_mask,start=(/lonstart,latstart/), count=(/nlongitude,nlatitude/)) + if(iret /= 0)then + call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/)) + endif + iret = nf90_close(ncid) + end if + + ! Note, these are only deallocated if use_clddet_zz = true + ! so their allocation should probably be moved inside the conditional... allocate( tbbp(nlongitude,nlatitude,nchan)) allocate( ter_sat(nlongitude,nlatitude)) allocate( SDob_10(nlongitude,nlatitude)) @@ -349,7 +401,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) do k=1,nchan iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id) - iret = nf90_get_var(ncid, tbb_id, tbbp(:,:,k)) + iret = nf90_get_var(ncid, tbb_id, tbbp(:,:,k),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0)then call da_warning(__FILE__,__LINE__, & (/"NETCDF4 read error for: Brightness Temperature"/)) @@ -379,7 +431,7 @@ subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp) 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(:,:)) + iret = nf90_get_var(ncid, ter_id, ter_sat(:,:),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/)) if(iret /= 0) print*,"NETCDF4 read error for: SatTer" where (ter_sat == -999.) ter_sat = missing_r iret = nf90_close(ncid) @@ -543,6 +595,12 @@ end if tb = tbb(ilongitude,ilatitude,k) end if + ! Apply clear-sky bias correction if requested + if ( ahi_apply_clrsky_bias ) then + tb = tb - satinfo(inst)%clearSkyBias(k) + !write(*,*)'channel, clearSkyBias = ',satinfo(inst)%ichan(k),satinfo(inst)%clearSkyBias(k) + endif + if (tb > 0.0) p % tb_inv(k) = tb end do @@ -594,6 +652,7 @@ end if p%solazi = sol_azi(ilongitude,ilatitude) p%sensor_index = inst p%ifgat = ifgat + p%cloudflag = int(cloud_mask(ilongitude,ilatitude)) ! AHI L-2 cloud flag allocate (p%next) ! add next data p => p%next nullify (p%next) @@ -611,6 +670,8 @@ end if deallocate(sol_zenith) deallocate(sat_azi) deallocate(sol_azi) + deallocate(cloud_mask) + if( use_clddet_zz ) deallocate(tbbp) if( use_clddet_zz ) deallocate(ter_sat) if( use_clddet_zz ) deallocate(SDob_10) 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 6e8a6acfeb..c5a6fa84dd 100644 --- a/var/da/da_radiance/da_write_iv_rad_ascii.inc +++ b/var/da/da_radiance/da_write_iv_rad_ascii.inc @@ -18,7 +18,11 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) character(len=filename_len) :: filename character(len=7) :: surftype integer :: ndomain - logical :: amsr2 + logical :: amsr2, ahi + real :: cip ! to output cloud-ice path + integer :: cloudflag ! to output cloudflag + integer, dimension(1) :: maxl + integer, allocatable :: level_max(:) real, allocatable :: dtransmt(:,:), transmt_jac(:,:), transmt(:,:), lod(:,:), lod_jac(:,:) @@ -48,7 +52,14 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) allocate ( lod_jac(iv%instid(i)%nchan,iv%instid(i)%nlevels) ) end if + ! Get variables for maximum level of weighting function + ! Only filled if calc_weightfunc .and. use_crtm_kmatrix both are true + if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then + allocate(level_max(iv%instid(i)%nchan)) + endif + amsr2 = index(iv%instid(i)%rttovid_string,'amsr2') > 0 + ahi = index(iv%instid(i)%rttovid_string,'ahi') > 0 write(unit=filename, fmt='(i2.2,a,i4.4)') it,'_inv_'//trim(iv%instid(i)%rttovid_string)//'.', myproc @@ -69,7 +80,7 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) 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 cloud_frac' + & soiltyp vegtyp vegfra elev clwp cloud_frac cip cloudflag' ndomain = 0 do n =1,iv%instid(i)%num_rad if (iv%instid(i)%info%proc_domain(1,n)) then @@ -118,7 +129,31 @@ 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,f15.5)') surftype, n, & + + ! Output cloud-ice path, cloudflag, pressure of weighting function peak + if (rtm_option==rtm_option_crtm .and. crtm_cloud ) then + cip = iv%instid(i)%cip(n) + else + cip = -9999.0 + endif + if ( ahi ) then + cloudflag = iv%instid(i)%cloudflag(n) + else + cloudflag = -999 + endif + + ! %nlevels is unstaggered, so subtract 1 to get number of mass points + ! Order is top down in CRTM + if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then + level_max(:) = -1 + do l=1,iv%instid(i)%nchan + maxl(:) = maxloc(iv%instid(i)%der_trans(l,:,n)) ! returns index of maximum value; if a tie, returns first occurrence of maximum + level_max(l) = maxl(1) ! model level of weighting function peak for this pixel and channel, going from top-->down + enddo + !level_max(:) = ( iv%instid(i)%nlevels - 1 ) - level_max(:) + 1 ! make order bottom-up + endif + + write(unit=innov_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3,f15.5,f15.5,i7)') surftype, n, & iv%instid(i)%t2m(n), & iv%instid(i)%mr2m(n), & iv%instid(i)%u10(n), & @@ -134,12 +169,22 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) iv%instid(i)%vegfra(n), & iv%instid(i)%elevation(n), & iv%instid(i)%clwp(n), & - iv%instid(i)%cloud_frac(n) + iv%instid(i)%cloud_frac(n), & + cip, & + cloudflag write(unit=innov_rad_unit,fmt='(a)') 'OBS : ' write(unit=innov_rad_unit,fmt='(10f11.2)') ob%instid(i)%tb(:,n) write(unit=innov_rad_unit,fmt='(a)') 'BAK : ' write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_xb(:,n) + if (rtm_option==rtm_option_crtm .and. crtm_cloud .and. (amsr2 .or. ahi) ) then + write(unit=innov_rad_unit,fmt='(a)') 'BAK_clr : ' + write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_xb_clr(:,n) + endif + if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then + write(unit=innov_rad_unit,fmt='(a)') 'WEIGHTFUNC_PEAK : ' + write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%pm( (/level_max/), n ) + endif write(unit=innov_rad_unit,fmt='(a)') 'IVBC : ' write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_inv(:,n) write(unit=innov_rad_unit,fmt='(a)') 'EMS : ' @@ -313,6 +358,9 @@ subroutine da_write_iv_rad_ascii (it, ob, iv ) deallocate ( lod ) deallocate ( lod_jac ) end if + if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then + deallocate(level_max) + endif close(unit=innov_rad_unit) call da_free_unit(innov_rad_unit) end do ! end do instruments diff --git a/var/da/da_radiance/module_radiance.f90 b/var/da/da_radiance/module_radiance.f90 index fab0afa39d..2fbfdd0a9c 100644 --- a/var/da/da_radiance/module_radiance.f90 +++ b/var/da/da_radiance/module_radiance.f90 @@ -169,6 +169,10 @@ module module_radiance real , pointer :: bcoef(:,:) ! airmass predictor bias coefficients real , pointer :: bcoef0(:) ! airmass constant coefficient real , pointer :: error_std(:) ! error standard deviation + real , pointer :: BTLim(:) ! for all-sky radiances, "BTLim" for each channel (Harnish et al. 2016) + real , pointer :: ca1(:) ! for all-sky radiances, symmetric cloud amount below which we set obs error to clear-sky obs error. + real , pointer :: ca2(:) ! for all-sky radiances, symmetric cloud amount above which we set obs error to fully cloudy obs error. + real , pointer :: clearSkyBias(:) ! for all-sky radiances, bias correction determined offline based on only clear pixels. end type satinfo_type type (satinfo_type), pointer :: satinfo(:) diff --git a/var/run/ahi_info b/var/run/ahi_info index 2d6da8f68b..f993ef833d 100644 --- a/var/run/ahi_info +++ b/var/run/ahi_info @@ -2,7 +2,7 @@ data source:1.cma hdf5;2.geocat netcdf4;3.jaxa netcdf4;4.ncep bufr 2 nscan 3000 -area information for geocat netcdf4 data: lonstart latstart nlongitude nlatitude +area information for geocat netcdf4 and jaxa netcdf4 data: lonstart latstart nlongitude nlatitude 1,1,1500,1000 date infomation for cma hdf5 data 2016,07,18,19,00,00 diff --git a/var/run/radiance_info/himawari-8-ahi.info b/var/run/radiance_info/himawari-8-ahi.info index 229563e49a..9b58c78ed0 100644 --- a/var/run/radiance_info/himawari-8-ahi.info +++ b/var/run/radiance_info/himawari-8-ahi.info @@ -1,11 +1,11 @@ -sensor channel IR/MW use idum varch polarisation(0:vertical;1:horizontal) - 478 1 1 -1 0 1.0520000000E+00 1.0000000000E+00 28.30175 - 478 2 1 1 0 1.7000000000E+00 1.0000000000E+00 57.58830 - 478 3 1 -1 0 1.7000000000E+00 1.0000000000E+00 12.69287 - 478 4 1 -1 0 1.3500000000E+00 1.0000000000E+00 27.33099 - 478 5 1 -1 0 0.8140000000E+00 1.0000000000E+00 23.24269 - 478 6 1 -1 0 0.9310000000E+00 1.0000000000E+00 53.35099 - 478 7 1 -1 0 0.8710000000E+00 1.0000000000E+00 36.07700 - 478 8 1 -1 0 0.9260000000E+00 1.0000000000E+00 33.61592 - 478 9 1 -1 0 0.9330000000E+00 1.0000000000E+00 33.61592 - 478 10 1 -1 0 0.7870000000E+00 1.0000000000E+00 33.61592 +sensor channel IR/MW use idum varch polarisation(0:vertical;1:horizontal) BTLim ca1 ca2 clearSkyBias + 478 7 1 -1 0 1.0520000000E+00 1.0000000000E+00 28.30175 273.30 0.5 14.5 -1.023 + 478 8 1 1 0 1.7000000000E+00 1.0000000000E+00 57.58830 238.00 1.0 25.0 0.000 + 478 9 1 -1 0 1.7000000000E+00 1.0000000000E+00 12.69287 230.35 0.5 14.5 -1.221 + 478 10 1 -1 0 1.3500000000E+00 1.0000000000E+00 27.33099 220.72 0.5 14.5 -1.621 + 478 11 1 -1 0 0.8140000000E+00 1.0000000000E+00 23.24269 220.72 0.5 14.5 -1.621 + 478 12 1 -1 0 0.9310000000E+00 1.0000000000E+00 53.35099 220.72 0.5 14.5 -1.621 + 478 13 1 -1 0 0.8710000000E+00 1.0000000000E+00 36.07700 220.72 0.5 14.5 -1.621 + 478 14 1 -1 0 0.9260000000E+00 1.0000000000E+00 33.61592 220.72 0.5 14.5 -1.621 + 478 15 1 -1 0 0.9330000000E+00 1.0000000000E+00 33.61592 220.72 0.5 14.5 -1.621 + 478 16 1 -1 0 0.7870000000E+00 1.0000000000E+00 33.61592 220.72 0.5 14.5 -1.621