Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

GitHub Issue NOAA-EMC/GSI#68. Switches assimilating Ta to Tb. #271

Closed
wants to merge 1 commit into from
+216 −39
Closed

GitHub Issue NOAA-EMC/GSI#68. Switches assimilating Ta to Tb. #271

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions scripts/exglobal_atmos_analysis.sh
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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