From 58246e82af34eb2ee1e82406129991ffb3d9f2dc Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Thu, 25 Jul 2024 12:16:19 +0000 Subject: [PATCH] These changes restructure the imager information for CrIS (VIIRS) and IASI (AVHRR) to avoid zero values for cluster size and radiance values. I also removed an unused variable. --- src/gsi/read_cris.f90 | 87 ++++++++++++++++++++++++------------------- src/gsi/read_iasi.f90 | 79 ++++++++++++++++++++++----------------- 2 files changed, 93 insertions(+), 73 deletions(-) diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 190ca24c5..84288a7f0 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -112,9 +112,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile + character(len=*) ,intent(in ) :: infile character(len=*) ,intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_cris @@ -200,7 +200,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(83,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev, imager_conversion - real(r_kind) :: imager_cluster_tot ! bufr error codes ! real(r_kind),dimension(7,3) :: error_codes @@ -848,10 +847,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& if ( cris_cads ) then call ufbseq(lnbufr,imager_info,83,7,iret,'CRISCS') if ( iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -871,51 +869,62 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster - iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster - imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) + if ( imager_cluster_size(i) > zero .and. imager_info(76,i) > zero .and. imager_info(81,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. - imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) + iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster + imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. + imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) - iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster - imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster + imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) - iexponent = -(nint(imager_info(82,i))-5 ) ! channel 16 radiance std dev for each cluster. - iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. - imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. + imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) - - - end do imager_cluster_info + call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else + data_all(maxinfo+j,itx) = zero ! something is wrong + data_all(maxinfo+7+j,itx) = zero ! set everything to zero + data_all(maxinfo+14+j,itx) = zero + endif + end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average - imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average - imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! cris_cads diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index a7e2e742b..75dc50bb7 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -142,9 +142,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile + character(len=*) ,intent(in ) :: infile character(len=*) ,intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_iasi @@ -233,7 +233,6 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(33,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev - real(r_kind) :: imager_cluster_tot ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" @@ -813,10 +812,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& if ( iasi_cads ) then call ufbseq(lnbufr,imager_info,33,7,iret,'IASIL1CS') if (iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -834,49 +832,62 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster + if ( imager_cluster_size(i) > zero .and. imager_info(26,i) > zero .and. imager_info(31,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. - imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. + imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. - imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. + imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster + imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster - imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. + imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. - imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - - call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else ! something is wrong + data_all(maxinfo+j,itx) = zero ! set everything to zero + data_all(maxinfo+7+j,itx) = zero + data_all(maxinfo+14+j,itx) = zero + endif end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(26,:)**2 + imager_info(28,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(31,:)**2 + imager_info(33,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! iasi_cads = .true. !