Skip to content

Commit

Permalink
These changes restructure the imager information for CrIS (VIIRS) and…
Browse files Browse the repository at this point in the history 
… IASI (AVHRR) to avoid zero values for cluster size and radiance values. I also removed an unused variable.
  • Loading branch information
@wx20jjung
wx20jjung committed Jul 25, 2024
1 parent 0db28da commit 58246e8
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 73 deletions.
87 changes: 48 additions & 39 deletions src/gsi/read_cris.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(:))
Expand All @@ -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

Expand Down
79 changes: 45 additions & 34 deletions src/gsi/read_iasi.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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(:))
Expand All @@ -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.
!
Expand Down

0 comments on commit 58246e8

Please sign in to comment.