Skip to content

Commit

Permalink
GitHub Issue #68. Switches assimilating Ta to Tb.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ScottSieron-NOAA
ScottSieron-NOAA committed Mar 2, 2022
1 parent 828c0cf commit 19fb52c
Show file tree
Hide file tree
Showing 4 changed files with 216 additions and 39 deletions.
1 change: 1 addition & 0 deletions scripts/exglobal_atmos_analysis.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -461,6 +461,7 @@ for file in $(awk '{if($1!~"!"){print $1}}' satinfo | sort | uniq); do
$NLN $RTMFIX/${file}.SpcCoeff.bin ./crtm_coeffs/
$NLN $RTMFIX/${file}.TauCoeff.bin ./crtm_coeffs/
done
$NLN $RTMFIX/amsua_metop-a_v2.SpcCoeff.bin ./crtm_coeffs/amsua_metop-a_v2.SpcCoeff.bin

$NLN $RTMFIX/Nalli.IRwater.EmisCoeff.bin ./crtm_coeffs/Nalli.IRwater.EmisCoeff.bin
$NLN $RTMFIX/NPOESS.IRice.EmisCoeff.bin ./crtm_coeffs/NPOESS.IRice.EmisCoeff.bin
Expand Down
101 changes: 101 additions & 0 deletions src/gsi/antcorr_application.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@
! paul.vandelst@noaa.gov
!
! 2011-04-25 A.Collard Modified to be consistent with CRTM 2.1
! 2020-12-13 S.Sieron Added Apply_Antcorr for assimilation of
! brightness temperatures instead of antenna
! temperatures
!

MODULE AntCorr_Application
Expand All @@ -35,6 +38,7 @@ MODULE AntCorr_Application
! The AntCorr structure definition
PUBLIC :: ACCoeff_type
PUBLIC :: Remove_AntCorr
PUBLIC :: Apply_AntCorr

! -----------------
! Module parameters
Expand Down Expand Up @@ -153,5 +157,102 @@ SUBROUTINE Remove_AntCorr( AC , & ! Input
T(l) = AC%A_earth(iFOV,l)*T(l) + AC%A_platform(iFOV,l)*T(l) + AC%A_space(iFOV,l)*TSPACE
END DO
END SUBROUTINE Remove_AntCorr

!--------------------------------------------------------------------------------
!
! NAME:
! Apply_AntCorr
!
! PURPOSE:
! Subroutine to apply an antenna correction to microwave instrument
! antenna temperatures, Ta, to produce brightness temperatures, Tb.
!
! CALLING SEQUENCE:
! CALL Apply_AntCorr( AC , & ! Input
! iFOV, & ! Input
! T ) ! In/Output
!
! INPUT ARGUMENTS:
! AC: Structure containing the antenna correction coefficients
! for the sensor of interest.
! UNITS: N/A
! TYPE: TYPE(ACCoeff_type)
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! iFOV: The FOV index for a scanline of the sensor of interest.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! T: On input, this argument contains the antenna
! temperatures for the sensor channels.
! UNITS: Kelvin
! TYPE: REAL(fp)
! DIMENSION: Rank-1 (n_Channels)
! ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUT ARGUMENTS:
! T: On output, this argument contains the brightness
! temperatures for the sensor channels.
! If an error occurs, the return values are all -1.
! UNITS: Kelvin
! TYPE: REAL(fp)
! DIMENSION: Rank-1 (n_Channels)
! ATTRIBUTES: INTENT(IN OUT)
!
! SIDE EFFECTS:
! The temperature array argument, T, is modified.
!
! PROCEDURE:
! For every FOV and channel, the antenna temperature, Ta, is converted
! to brightness temperature, Tb, using,
!
! Tb = (Ta - As.Ts) / (Ae + Ap)
!
! where Ae == antenna efficiency for the Earth view
! Ap == antenna efficiency for satellite platform view
! As == antenna efficiency for cold space view
! Ts == cosmic background temperature.
!
! Note that the observed earth view brightness temperature is used as a
! proxy for the platform temperature for the (Ap.Tb) component since
! there is no measurement of the platform temperature in-flight.
!
!--------------------------------------------------------------------------------

SUBROUTINE Apply_AntCorr( AC , & ! Input
iFOV, & ! Input
T ) ! In/Output
implicit none

! Arguments
TYPE(ACCoeff_type), INTENT(IN) :: AC
INTEGER , INTENT(IN) :: iFOV
REAL(fp) , INTENT(IN OUT) :: T(:)
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Apply_AntCorr'
! Local variables
INTEGER :: l

! Check input
IF ( iFOV < 1 .OR. iFOV > AC%n_FOVS ) THEN
CALL Display_Message( ROUTINE_NAME, 'Input iFOV inconsistent with AC data', FAILURE )
T = INVALID
RETURN
END IF
IF ( SIZE(T) /= AC%n_Channels ) THEN
CALL Display_Message( ROUTINE_NAME, 'Size of T() inconsistent with AC data', FAILURE )
T = INVALID
RETURN
END IF
! Compute the brightness temperature
DO l = 1, AC%n_Channels
T(l) = (T(l) - AC%A_space(iFOV,l)*TSPACE) / (AC%A_earth(iFOV,l) + &
AC%A_platform(iFOV,l))
END DO
END SUBROUTINE Apply_AntCorr


END MODULE AntCorr_Application
8 changes: 5 additions & 3 deletions src/gsi/read_atms.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
! 2018-02-05 collard - get orbit height from BUFR file
! 2018-04-19 eliu - allow data selection for precipitation-affected data
! 2018-05-21 j.jin - added time-thinning, to replace thin4d
! 2020-12-13 s.sieron - change from reading antenna to brightness temperature
!
! input argument list:
! mype - mpi task id
Expand Down Expand Up @@ -479,10 +480,11 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
solazi_save(iob)=bfr2bhdr(4)

! Read data record. Increment data counter
! TMBR is actually the antenna temperature for most microwave sounders but for
! ATMS it is stored in TMANT.
! Though TMBR is actually the antenna temperature for most microwave sounders,
! for ATMS, TMBR is brightness temperature. (Antenna temperature is in
! TMANT.)
! ATMS is assumed not to come via EARS
call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMANT')
call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR')

bt_save(1:nchanl,iob) = data1b8(1:nchanl)

Expand Down
145 changes: 109 additions & 36 deletions src/gsi/read_bufrtovs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,9 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
! channels are missing.
! 2018-04-19 eliu - allow data selection for precipitation-affected data
! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90.
! 2020-12-13 s.sieron - change to converting antenna temperatures to brightness temperautres.
! added support for multiple SpcCoeff files (antenna correction
! coefficients).
!
! input argument list:
! mype - mpi task id
Expand Down Expand Up @@ -141,9 +144,10 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
crtm_kind => fp, &
MAX_SENSOR_ZENITH_ANGLE
use crtm_spccoeff, only: sc,crtm_spccoeff_load,crtm_spccoeff_destroy
use ACCoeff_Define, only: ACCoeff_type
use calc_fov_crosstrk, only : instrument_init, fov_cleanup, fov_check
use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen
use antcorr_application, only: remove_antcorr
use antcorr_application, only: remove_antcorr, apply_antcorr
use mpeu_util, only: getindex
use deter_sfc_mod, only: deter_sfc_fov,deter_sfc
use gsi_nstcouplermod, only: nst_gsi,nstinfo
Expand Down Expand Up @@ -174,7 +178,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
! Declare local parameters

character(8),parameter:: fov_flag="crosstrk"
integer(i_kind),parameter:: n1bhdr=13
integer(i_kind),parameter:: n1bhdr=14
integer(i_kind),parameter:: n2bhdr=4
real(r_kind),parameter:: r360=360.0_r_kind
real(r_kind),parameter:: tbmin=50.0_r_kind
Expand All @@ -190,6 +194,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&

integer(i_kind) ireadsb,ireadmg,irec,next,nrec_startx
integer(i_kind) i,j,k,ifov,ntest,llll
integer(i_kind) sacv
integer(i_kind) iret,idate,nchanl,n,idomsfc(1)
integer(i_kind) ich1,ich2,ich8,ich15,ich16,ich17
integer(i_kind) kidsat,instrument,maxinfo
Expand Down Expand Up @@ -234,6 +239,11 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
logical :: critical_channels_missing,quiet
logical :: print_verbose

logical :: spc_coeff_found
integer(i_kind) :: spc_coeff_versions
character(len=80) :: spc_filename
type(ACCoeff_type),dimension(5) :: accoeff_sets

!**************************************************************************
! Initialize variables

Expand Down Expand Up @@ -476,7 +486,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
if(dval_use) maxinfo=maxinfo+2
nreal = maxinfo + nstinfo
nele = nreal + nchanl
hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS'
hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS SACV'
hdr2b ='SAZA SOZA BEARAZ SOLAZI'
allocate(data_all(nele,itxmax),data1b8(nchanl),data1b4(nchanl),nrec(itxmax))

Expand Down Expand Up @@ -506,34 +516,73 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&

call openbf(lnbufr,'IN',lnbufr)

if(llll >= 2 .and. (amsua .or. amsub .or. mhs))then
if(amsua .or. amsub .or. mhs)then
quiet=.not.verbose
allocate(data1b8x(nchanl))
sensorlist(1)=sis
if( crtm_coeffs_path /= "" ) then
if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_BUFRTOVS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"'
error_status = crtm_spccoeff_load(sensorlist,&
File_Path = crtm_coeffs_path, quiet=quiet )
spc_coeff_versions = 0
spc_coeff_found = .true.
do while (spc_coeff_found)

! determine the name of next SpcCoeff file to test for existence
if (spc_coeff_versions == 0) then
sensorlist(1)=sis
else
error_status = crtm_spccoeff_load(sensorlist,quiet=quiet)
endif
if (error_status /= success) then
write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,&
' TERMINATE PROGRAM EXECUTION'
call stop2(71)
endif
ninstruments = size(sc)
instrument_loop: do n=1,ninstruments
if(sis == sc(n)%sensor_id)then
instrument=n
exit instrument_loop
i = spc_coeff_versions+1
write(sensorlist(1),'(a,a,i1)') trim(sis),'_v',i
end if
end do instrument_loop
if(instrument == 0)then
write(6,*)' failure to find instrument in read_bufrtovs ',sis
end if
end if

if( crtm_coeffs_path /= "" ) then
if(mype_sub==mype_root .and. print_verbose) &
write(6,*)'READ_BUFRTOVS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"'
end if

! test SpcCoeff file for existence
spc_filename = trim(crtm_coeffs_path) // trim(sensorlist(1)) // '.SpcCoeff.bin'
INQUIRE(FILE=trim(spc_filename), EXIST=spc_coeff_found)

if (.NOT. spc_coeff_found) then
if (spc_coeff_versions == 0) then
write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load no versiosn of SpcCoeff found for ', trim(sis)
call stop2(71)
else
write(6,*)'READ_BUFRTOVS: ', spc_coeff_versions, ' versions of SpcCoeff found for ', trim(sis)
end if
else
spc_coeff_versions = spc_coeff_versions+1

if( crtm_coeffs_path /= "" ) then
error_status = crtm_spccoeff_load(sensorlist,&
File_Path = crtm_coeffs_path, quiet=quiet )
else
error_status = crtm_spccoeff_load(sensorlist,quiet=quiet)
endif
if (error_status /= success) then
write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,&
' despite file ',trim(spc_filename),' existing, TERMINATE PROGRAM EXECUTION'
call stop2(71)
endif

ninstruments = size(sc)
instrument_loop: do n=1,ninstruments
if(sis == sc(n)%sensor_id)then
instrument=n
exit instrument_loop
end if
end do instrument_loop
if(instrument == 0)then
write(6,*)' failure to find instrument in read_bufrtovs ',sis
end if

accoeff_sets(spc_coeff_versions) = sc(instrument)%ac

! deallocate crtm info
error_status = crtm_spccoeff_destroy()
if (error_status /= success) &
write(6,*)'OBSERVER: ***ERROR*** crtm_spccoeff_destroy error_status=',error_status
end if
end do

end if

! Loop to read bufr file
irecx=0
Expand Down Expand Up @@ -633,6 +682,17 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh)
if(.not. iuse)cycle read_loop

! Extract satellite antenna corrections version number
if (llll > 1) then
sacv = nint(bfr1bhdr(14))
if (sacv > spc_coeff_versions) then
write(6,*) 'READ_BUFRTOVS WARNING sacv greater than spc_coeff_versions'
end if
else ! normal feed doesn't have antenna correction, so set sacv to 0
sacv = 0
end if


call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b)

! Set common predictor parameters
Expand Down Expand Up @@ -675,18 +735,36 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&

! Read data record. Increment data counter
! TMBR is actually the antenna temperature for most microwave
! sounders.
! sounders. EARS / DB has brightness temperature (TMBRST).
! Convert TMBR to brightness temperature (add_antcorr).
! Accept TMBRST. But if the SACV version isn't 1, convert to
! antenna temperature, and convert back to brightness temperature.
if (llll == 1) then
call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR')
if ( amsua .or. amsub .or. mhs )then
! convert antenna temperature to brightness temperature
data1b8x=data1b8
data1b4=data1b8
call apply_antcorr(accoeff_sets(1),ifov,data1b4)
data1b8=data1b4
do j=1,nchanl
if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind
end do
end if
else ! EARS / DB
call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST')
if ( amsua .or. amsub .or. mhs )then
if ( amsua .or. amsub .or. mhs .AND. spc_coeff_versions /= 1 .AND. sacv /= 1)then
! convert brightness temperature to antenna temperature using
! the satellite antenna correction version (sacv) used by the
! data originator, then convert back to brightness temperature
! using the version of parameters used by the CRTM (v1)
data1b8x=data1b8
data1b4=data1b8
call remove_antcorr(sc(instrument)%ac,ifov,data1b4)
call remove_antcorr(accoeff_sets(sacv),ifov,data1b4)
call apply_antcorr(accoeff_sets(1),ifov,data1b4)
data1b8=data1b4
do j=1,nchanl
if(data1b8x(j) > r1000)data1b8(j) = 1000000._r_kind
if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind
end do
end if
end if
Expand Down Expand Up @@ -947,13 +1025,8 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
call closbf(lnbufr)
close(lnbufr)

if(llll > 1 .and. (amsua .or. amsub .or. mhs))then
if(amsua .or. amsub .or. mhs)then
deallocate(data1b8x)

! deallocate crtm info
error_status = crtm_spccoeff_destroy()
if (error_status /= success) &
write(6,*)'OBSERVER: ***ERROR*** crtm_spccoeff_destroy error_status=',error_status
end if

end do ears_db_loop
Expand Down

0 comments on commit 19fb52c

Please sign in to comment.