module sst_retrieval !$$$ module documentation block ! . . . . ! module: sst_retrieval sst retrieval algorithm and realted subroutines ! prgmmr: li, xu org: np23 date: 2004-12-21 ! ! abstract: SST retrieval setup and application ! ! program history log: ! 2004-12-21 xu li ! ! Subroutines Included: ! sub setup_sst_retrieval - define parameters for SST retrieval ! sub avhrr_sst_retrieval - perform SST retrieval with AVHRR radiance, save results ! sub spline_cub - cubic spline interpolation for SST dependent bias correction ! sub finish_sst_retrieval - close bufr output file ! ! Functions Included: ! ! Variable Definitions: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_single,i_kind use radinfo, only: numt implicit none ! set default to private private ! set subroutines to public public :: setup_sst_retrieval public :: avhrr_sst_retrieval public :: spline_cub public :: finish_sst_retrieval ! Define parameters integer(i_kind),parameter:: lnbufr = 12_i_kind ! Declare variables and arrays real(r_kind), dimension(0:1) :: e_ts,e_ta,e_qa real(r_kind), dimension(numt):: x_bias_sst contains subroutine setup_sst_retrieval(obstype,csatid,mype) !$$$ subprogram documentation block ! . . . . ! subprogram: setup_sst_retrieval compute rhs of oi for sst retrieval ! prgmmr: xu li org: np23 date: 2004-12-21 ! ! abstract: set up sst retrieval ! ! program history log: ! 2004-12-21 xu li ! 2005-02/16 xu li - modify the way to define bufr table for SST ! physical retrieval ! 2006-02-03 derber - modify for new obs control and obs count ! 2006-04-06 middlecoff - change lundx from 20 to lendian_in so ! can be read as little endian ! 2006-07-28 derber - modify hanling of lextra ! 2008-04-11 safford - rm unused vars and uses ! ! input argument list: ! obstype - type of tb observation ! csatid - satellite id ! mype - mpi task id ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use constants, only: ione,one,half,ttp use gsi_io, only: lendian_in implicit none ! Define parameters character(10) , intent(in ) :: obstype,csatid integer(i_kind), intent(in ) :: mype integer(i_kind) i character(9) :: string character(len=80) :: bufrtabf character(30) bufr_sst ! bufr output for sst2dvar input ! Open output bufr file and table for SST retrieval write(string,99)mype 99 format('.',I4.4) bufr_sst= 'bfsst_'//trim(obstype)//'.'//trim(csatid)//string open(lnbufr,file=trim(bufr_sst),form='unformatted') ! open bufr data file bufrtabf = 'bftab_sstphr' open(lendian_in,file=trim(bufrtabf)) ! open bufr table call openbf (lnbufr,'NODX',lendian_in) ! Assign x cooridinate for SST dependent AVHRR radiance bias correction do i = 1, numt x_bias_sst(i) = ttp - one + real(i-ione) enddo ! Assign error parameters for background (Ts, Ta, Qa) ! ! Day time e_ts(0) = half; e_ta(0) = 1.20_r_kind; e_qa(0) = 0.95_r_kind ! Night time ! e_ts(1) = 0.45_r_kind; e_ta(1) = 0.90_r_kind; e_qa(1) = 0.65_r_kind end subroutine setup_sst_retrieval subroutine avhrr_sst_retrieval(obstype,csatid,nchanl,& tnoise,varinv,ts5,sstnv,sstcu,temp,wmix,ts,tbc,obslat,obslon,zasat,& dtime,dtp_avh,tlapchn,predterms,emissivity,pangs,tbcnob,tb_obs, & rad_diagsave,sfcpct,id_qc,nadir,ireal,ipchan,duse) !subprogram: avhrr_sst_retrieval compute sst retrieval from AVHRR radiances ! prgmmr: Xu Li, John Derber org: w/nmc2 date: 04-12-29 ! ! abstract: perform sst retrieval based on input radiative transfer info ! save bufr output and satellite diagnostics files ! ! program history log: ! 2004-12-29 li ! 2005-04-18 treadon - add passed logical lextra ! 2005-09-28 derber - modify land mask stuff ! 2006-04-27 derber - modify to do single profile ! 2008-04-11 safford - rm unused vars and uses ! ! input argument list: ! obstype - observation type ! nsig - number of model layers ! csatid - satellite id ! nchanl - number of channels for instruments ! tnoise - error of observed radiance ! varinv - inverse error squared ! ts5 - sst used in Radiative transfer and first guess for SST retrieval ! sstnv - navy sst retrieval ! sstcu - ncep rtg sst analysis for the current day ! temp - d(brightness temperature)/d(temperature) ! wmix - d(brightness temperature)/d(mixing ratio) ! ts - d(brightness temperature)/d(skin temperature) ! tbc - bias corrected (observed - simulated brightness temperatures) ! obslat - latitude of observations ! obslon - longitude of observations ! iadate - analysis time window (00, 06, 12, 18) ! dtime - observed time distance from iadate ! dtp_avh - night/day mode ! rad_diagsave - logical to switch on diagnostic output (.false.=no output) ! nadir - number of FOV ! ireal - rank of array diagbuf ! ipchan - The first number (7) of elements in array diagbufchan ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use radinfo, only: npred use gridmod, only: nsig use obsmod, only: iadate,rmiss_single use constants, only: izero,ione,zero,one,tiny_r_kind,rad2deg,ttp,one_tenth implicit none ! Declare passed variables integer(i_kind) , intent(in ) :: nchanl,ireal,ipchan integer(i_kind) , intent(in ) :: nadir real(r_kind),dimension(nchanl) , intent(in ) :: tnoise real(r_kind),dimension(nchanl) , intent(in ) :: varinv,tbc,tb_obs,tbcnob,tlapchn real(r_kind),dimension(nchanl) , intent(in ) :: emissivity real(r_kind) , intent(in ) :: obslat,obslon,pangs,zasat,dtime,ts5,& sstnv,sstcu,dtp_avh real(r_kind),dimension(nchanl) , intent(in ) :: ts real(r_kind),dimension(nsig,nchanl) , intent(in ) :: wmix,temp real(r_kind),dimension(npred+ione,nchanl), intent(in ) :: predterms character(10) , intent(in ) :: obstype,csatid logical , intent(in ) :: rad_diagsave logical , intent(in ) :: duse ! Declare local parameters integer(i_kind),parameter:: nmsg=11_i_kind character(len=8), parameter :: subset='NC012017' real(r_kind),parameter:: r180=180.0_r_kind real(r_kind),parameter:: r360=360.0_r_kind ! Declare local variables real(r_kind) :: ws,wa,wq,errinv integer(i_kind) :: icount,idate,md,ch3,ch4,ch5,iret,nn,i,j,l integer(i_kind) :: iextra,jextra real(r_kind), dimension(nchanl) :: tb_ta,tb_qa real(r_kind), dimension(nchanl) :: w_avh real(r_kind) :: delt,delt1,delt2,delt3,c1x,c2x,c3x real(r_kind) :: a11,a12,a13,a23,a22,a33 real(r_kind) :: bt3,bt4,bt5,varrad,rsat real(r_kind) :: nuse logical :: lextra ! Declare local arrays real(r_kind) :: sstph,dsst,dqa,dta integer(i_kind),dimension(8) :: idat8,jdat8 real(r_kind), dimension(5) :: rinc real(r_kind), dimension(4) :: sfcpct real(r_kind), dimension(nmsg) :: bufrf integer(i_kind),dimension(nchanl):: id_qc real(r_single),dimension(ireal) :: diagbuf real(r_single),dimension(ipchan+npred+ione,nchanl) ::diagbufchan real(r_single),allocatable,dimension(:,:) :: diagbufex ch3=ione;ch4=2_i_kind;ch5=3_i_kind lextra=.false. iextra=izero jextra=izero if(obstype == 'avhrr' .or. obstype == 'avhrr_navy')then lextra=.true. iextra=3_i_kind jextra=nchanl end if allocate (diagbufex(iextra,jextra)) !********************************** ! Get tb_ta & tb_qa !********************************** do i = 1, nchanl tb_ta(i) = temp(1,i) tb_qa(i) = wmix(1,i) do l = 2, nsig tb_ta(i) = tb_ta(i) + temp(l,i) tb_qa(i) = tb_qa(i) + wmix(l,i) enddo enddo if (duse) then if (duse) then nuse = one else nuse = zero endif dsst = rmiss_single; dta = rmiss_single; dqa = rmiss_single ! Get day/night mode if ( dtp_avh /= 152.0_r_kind ) then md = izero ! Day time else md = ione ! Night time endif ws = one/e_ts(md)**2 wa = one/e_ta(md)**2 wq = one/(e_qa(md)*(max((ts5-ttp)*0.03_r_kind,zero)+one_tenth))**2 a11 = ws ! 1./tserr**2 a22 = wa ! 1./taerr**2 a33 = wq ! 1./qaerr**2 a12 = zero; a13 = zero; a23 = zero c1x = zero; c2x = zero; c3x = zero delt1 = zero; delt = one; icount = izero do l=1,nchanl ! Get coefficients for linear equations if (varinv(l) > tiny_r_kind) then icount = icount+ione w_avh(l) = (one/tnoise(l))**2 a11 = a11 + w_avh(l)*ts(l)**2 a12 = a12 + w_avh(l)*ts(l)*tb_ta(l) a13 = a13 + w_avh(l)*ts(l)*tb_qa(l) a22 = a22 + w_avh(l)*tb_ta(l)**2 a23 = a23 + w_avh(l)*tb_ta(l)*tb_qa(l) a33 = a33 + w_avh(l)*tb_qa(l)**2 varrad=w_avh(l)*tbc(l) c1x = c1x + varrad*ts(l) c2x = c2x + varrad*tb_ta(l) c3x = c3x + varrad*tb_qa(l) end if ! if(varinv(l) > tiny_r_kind) then end do ! do l=1,nchanl ! Solve linear equations with three unknowns (dsst, dta, dqa) if( (dtp_avh == 152._r_kind .and. icount > 2_i_kind ) .or. & (dtp_avh /= 152._r_kind .and. icount > ione ) ) then delt = a11*(a22*a33-a23*a23) + & a12*(a13*a23-a12*a33) + & a13*(a12*a23-a13*a22) delt1 = c1x*(a22*a33-a23*a23) + & c2x*(a13*a23-a12*a33) + & c3x*(a12*a23-a13*a22) delt2 = c1x*(a13*a23-a12*a33) + & c2x*(a11*a33-a13*a13) + & c3x*(a12*a13-a11*a23) delt3 = c1x*(a12*a23-a13*a22) + & c2x*(a13*a12-a11*a23) + & c3x*(a11*a22-a12*a12) dsst = delt1/delt dta = delt2/delt dqa = delt3/delt sstph = ts5 + dsst ! SST retrieval : SSTPH bt3 = tb_obs(ch3) - tbc(ch3) bt4 = tb_obs(ch4) - tbc(ch4) bt5 = tb_obs(ch5) - tbc(ch5) ! if ( varinv(ch3) == zero ) then ! bt3 = rmiss_single ! endif ! write(*,'(15f8.2,4F6.1,F3.0)') obslon,obslat,sstph,sstnv,ts5, & ! tb_obs(ch3),bt3,tb_obs(ch4),bt4,tb_obs(ch5),bt5, & ! sstph-ts5,dta,dqa,pangs,varinv(ch3),varinv(ch4),varinv(ch5),zasat,nuse ! ts(ch3),ts(ch4),ts(ch5), & ! tb_ta(ch3),tb_ta(ch4),tb_ta(ch5), & ! tb_qa(ch3),tb_qa(ch4),tb_qa(ch5), & ! Save bufr message for physical SST retrieval with Navy AVHRR radiance idat8(1) = iadate(1) ! 4-digit year idat8(2) = iadate(2) ! months of a year idat8(3) = iadate(3) ! days of a month idat8(4) = izero ! time zone idat8(5) = iadate(4) ! hours of a day idat8(6) = izero ! minutes of a hour idat8(7) = izero ! seconds idat8(8) = izero ! milliseconds rinc(1) = dtime/24._r_kind rinc(2) = zero rinc(3) = zero rinc(4) = zero rinc(5) = zero call w3movdat(rinc,idat8,jdat8) ! TRANSFORM THE NAVY AVHRR(GAC/NAVY) RECORD TO BUFR FORMAT ! ------------------------------------------------------------------------------ ! NC012015 | YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAID SIID SLNM FOVN | ! NC012015 | HMSL SAZA SOZA SSTTYP IREL IRMS OPTH SST1 SST2 SST3 | ! NC012015 | ATB1 ATB2 ATB3 ATB4 ATB5 DTB1 DTB2 DTB3 DTB4 DTB5 | ! ------------------------------------------------------------------------------ bufrf( 1) = jdat8( 1) ! 4-digit year bufrf( 2) = jdat8( 2) ! month of a year bufrf( 3) = jdat8( 3) ! day of a month bufrf( 4) = jdat8( 5) ! hour of a day bufrf( 5) = jdat8( 6) ! minute of a hour bufrf( 6) = jdat8( 7) ! second of a minute bufrf( 7) = obslat ! Latitude bufrf( 8) = obslon ! Longitude if(bufrf(8) > r180) bufrf(8) = bufrf(8) - r360 if(csatid == 'n15') rsat=206.0_r_kind if(csatid == 'n16') rsat=207.0_r_kind if(csatid == 'n17') rsat=208.0_r_kind if(csatid == 'n18') rsat=209.0_r_kind bufrf(9) = rsat ! Satellite ID (Data source) bufrf(10) = dtp_avh ! Type of observation/retrieval (151, 152, 159) bufrf(11) = sstph ! Physical SST retrieval ! WRITE THIS ARRAY INTO BUFR ! -------------------------- idate=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000 CALL OPENMB(lnbufr,subset,idate) CALL UFBSEQ(lnbufr,bufrf,nmsg,ione,iret,subset) CALL WRITSB(lnbufr) ! Save diagnostics at the satellite observation locations if (rad_diagsave) then ! Write diagnostics to output file. Only generate ! file on first outer iteration. diagbuf(1) = obslat ! observation latitude (degrees) diagbuf(2) = obslon ! observation longitude (degrees) diagbuf(3) = rmiss_single ! model (guess) elevation at observation location diagbuf(4) = dtime ! observation time (hours relative to analysis time) diagbuf(5) = nadir ! sensor scan position diagbuf(6) = zasat*rad2deg ! satellite zenith angle (degrees) diagbuf(7) = rmiss_single ! satellite azimuth angle (degrees) diagbuf(8) = pangs ! solar zenith angle (degrees) diagbuf(9) = rmiss_single ! solar azimuth angle (degrees) diagbuf(10) = rmiss_single ! sun glint angle (degrees) (sgagl) diagbuf(11) = sfcpct(1) ! fractional coverage by water diagbuf(12) = sfcpct(2) ! fractional coverage by land diagbuf(13) = sfcpct(3) ! fractional coverage by ice diagbuf(14) = sfcpct(4) ! fractional coverage by snow diagbuf(15) = ts5 ! SST first guess used for SST retrieval diagbuf(16) = sstcu ! NCEP SST analysis at t diagbuf(17) = sstph ! Physical SST retrieval diagbuf(18) = sstnv ! Navy SST retrieval diagbuf(19) = dta ! d(ta) corresponding to sstph diagbuf(20) = dqa ! d(qa) corresponding to sstph diagbuf(21) = dtp_avh ! data type diagbuf(22) = rmiss_single ! vegetation fraction diagbuf(23) = rmiss_single ! snow depth diagbuf(24) = rmiss_single ! surface wind speed (m/s) ! Note: The following quantities are not computed for all sensors diagbuf(25) = rmiss_single ! cloud fraction (%) diagbuf(26) = rmiss_single ! cloud top pressure (hPa) do i=1,nchanl diagbufchan(1,i)=tb_obs(i) ! observed brightness temperature (K) diagbufchan(2,i)=tbc(i) ! final observed minus simulated difference (K) diagbufchan(3,i)=tbcnob(i) ! initial observed minus simulated difference (K) errinv = sqrt(varinv(i)) diagbufchan(4,i)=errinv ! inverse observation error diagbufchan(5,i)=id_qc(i) ! quality control mark or event diagbufchan(6,i)=emissivity(i) ! surface emissivity diagbufchan(7,i)=tlapchn(i) ! stability index if ( nn == ione .and. varinv(nn) == zero ) then diagbufchan(1,nn) = rmiss_single diagbufchan(2,nn) = rmiss_single diagbufchan(3,nn) = rmiss_single diagbufchan(4,nn) = rmiss_single diagbufchan(5,nn) = rmiss_single diagbufchan(6,nn) = rmiss_single diagbufchan(7,nn) = rmiss_single endif do j=1,npred+ione diagbufchan(7_i_kind+j,i)=predterms(j,i) ! brightness temperature bias correction terms (K) end do if (lextra) then do nn=1,nchanl diagbufex(1,nn)=ts(nn) ! d(Tb)/d(Ts) diagbufex(2,nn)=tb_ta(nn) ! d(Tb)/d(Ta) diagbufex(3,nn)=tb_qa(nn) ! d(Tb)/d(Qa) end do endif end do if (.not.lextra) then write(4) diagbuf,diagbufchan else write(4) diagbuf,diagbufchan,diagbufex endif end if ! if (rad_diagsave) then end if ! if( (dtp_avh == 152. .and. icount > 2 ) .or. & end if ! if(duse) then end subroutine avhrr_sst_retrieval subroutine spline_cub(y,xs,ys) !$$$ subprogram documentation block ! . . . . ! subprogram: spline_cub one dimensional cubic spline interpolation ! prgmmr: Xu Li org: np23 date: 2003-09-01 ! ! abstract: one dimensional cubic spline interpolation ! ! program history log: ! 2003-09-01 xu li ! 2008-04-11 safford - rm unused vars and uses ! ! input argument list: ! y - cubic spline value at knots ! xs - interpolator ! ! output argument list: ! xs - interpolator ! ys - interpolatee ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use constants, only: ione,zero,half,one,two,three implicit none ! Passed variables real(r_kind) ,intent(inout) :: xs real(r_kind) ,intent( out) :: ys real(r_kind), dimension(numt),intent(in ) :: y ! Local variables integer(i_kind) i,j real(r_kind), dimension(numt) :: g,a,d,m,k ! handle the case when xs beyond [-1,31] if ( xs < x_bias_sst(1) ) xs = x_bias_sst(1) if ( xs > x_bias_sst(numt) ) xs = x_bias_sst(numt) do i = 1, numt - ione g(i) = one/(x_bias_sst(i+ione) - x_bias_sst(i)) enddo a(1) = half do i = 2, numt - ione a(i) = g(i)/(two*(g(i-ione)+g(i)) - a(i-ione)*g(i-ione)) enddo d(1) = 1.5_r_kind*(y(2) - y(1))/(x_bias_sst(2) - x_bias_sst(1)) do i = 2, numt - ione d(i) = (three*(g(i-ione)*(y(i)-y(i-ione))/(x_bias_sst(i)-x_bias_sst(i-ione))+g(i)*(y(i+ione)-y(i))/(x_bias_sst(i+ione)-x_bias_sst(i))) & -g(i-ione)*d(i-ione))/(two*(g(i-ione)+g(i))-g(i-ione)*a(i-ione)) enddo m(numt) = (three*(y(numt)-y(numt-ione))/(x_bias_sst(numt)-x_bias_sst(numt-ione)) - g(numt-ione)*d(numt-ione)) & /(g(numt-ione)*(two-a(numt-ione))) do i = numt-ione, 1, -1 m(i) = d(i) - a(i)*m(i+ione) enddo k(1) = zero k(numt) = zero do i = 2, numt - ione k(i) = 6.0_r_kind*g(i)*(y(i+ione)-y(i))/(x_bias_sst(i+ione)-x_bias_sst(i)) - two*g(i)*(two*m(i)+m(i+ione)) enddo do j = 1, numt-ione if ( xs >= x_bias_sst(j) .and. xs < x_bias_sst(j+ione) ) then i = j endif enddo ys = y(i)+m(i)*(xs-x_bias_sst(i))+half*k(i)*(xs-x_bias_sst(i))**2 & + ((k(i+ione)-k(i))/((x_bias_sst(i+ione)-x_bias_sst(i))*6.0_r_kind))*(xs-x_bias_sst(i))**3 end subroutine spline_cub subroutine finish_sst_retrieval !$$$ subprogram documentation block ! . . . . ! subprogram: finish_sst_retrieval close open bufr file ! prgmmr: treadon org: np23 date: 2005-01-07 ! ! abstract: Finish sst retrieval code by closing bufr file ! ! program history log: ! 2005-01-07 treadon ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block implicit none call closbf(lnbufr) return end subroutine finish_sst_retrieval end module sst_retrieval