diff --git a/scripts/exglobal_atmos_analysis.sh b/scripts/exglobal_atmos_analysis.sh index 87a0ab4f53..5bb01cc9b5 100755 --- a/scripts/exglobal_atmos_analysis.sh +++ b/scripts/exglobal_atmos_analysis.sh @@ -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 diff --git a/src/gsi/antcorr_application.f90 b/src/gsi/antcorr_application.f90 index db2808d08b..12eb7a4488 100644 --- a/src/gsi/antcorr_application.f90 +++ b/src/gsi/antcorr_application.f90 @@ -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 @@ -35,6 +38,7 @@ MODULE AntCorr_Application ! The AntCorr structure definition PUBLIC :: ACCoeff_type PUBLIC :: Remove_AntCorr + PUBLIC :: Apply_AntCorr ! ----------------- ! Module parameters @@ -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 diff --git a/src/gsi/read_atms.f90 b/src/gsi/read_atms.f90 index 6fc0006c93..cb87e26898 100644 --- a/src/gsi/read_atms.f90 +++ b/src/gsi/read_atms.f90 @@ -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 @@ -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) diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index cf3e7ea17c..7100b179a1 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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