subroutine read_amsre(mype,val_amsre,ithin,isfcalc,rmesh,gstime,&
infile,lunout,obstype,nread,ndata,nodata,twind,sis,&
mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs,nrec_start,dval_use)
! subprogram: read_amsre read bufr format amsre data
! prgmmr : okamoto org: np20 date: 2004-10-12
!
! abstract: This routine reads BUFR format AQUA radiance (brightness
! temperature) files. Optionally, the data are thinned to
! a specified resolution using simple quality control checks.
!
! When running the gsi in regional mode, the code only
! retains those observations that fall within the regional
! domain
!
! Special Notes:
! 10/14/04 okamoto looks like AMSRE overlap problem is not as bad as SSM/I
!
! program history log:
! 2004-10-12 okamoto
! 2005-10-07 Xu and Pawlak - modify the code related to ityp determination to
! use routine deter_ityp, created special notes section, removed
! GrADS specific code, fixed indentation
! 2005-10-10 treadon - replace deter_ityp with deter_sfc, modify rlndsea to be
! consistent with other read_* routines
! 2005-10-17 treadon - add grid and earth relative obs location to output file
! 2005-10-18 treadon - remove array obs_load and call to sumload
! 2005-10-20 kazumori - modify to read real AMSR-E bufr data and add zensun
! 2005-11-17 kazumori - add deter_sfc_amsre_low for AMSR-E low frequency channel
! 2005-11-29 parrish - modify getsfc to work for different regional options
! 2006-02-01 parrish - remove getsfc (different version called now in read_obs)
! 2006-02-02 kazumori - modify the threshold of surface determination and change
! the origin of satellite azimuth angle for diag file
! 2006-02-03 derber - modify for new obs control and obs count
! 2006-04-26 kazumori - bug fix of order of polarization, timedif, tbmax
! 2006-05-17 kazumori - modify for new bufrtable change
! 2006-05-19 eliu - add logic to reset relative weight when all channels not used
! 2006-08-05 kazumori - add good fov selection for amsre low channel data
! 2006-09-20 treadon - remove good fov selection for amsre low channel
! data in order to add mpi_io for data read
! 2006-10-22 kazumori - bug fix for the type of zensun subroutine argument
! 2007-03-01 tremolet - measure time from beginning of assimilation window
! 2008-05-28 safford - rm unused vars
! 2009-04-18 woollen - improve mpi_io interface with bufrlib routines
! 2009-04-21 derber - add ithin to call to makegrids
! 2009-09-01 li - add to handle nst fields and read/use data over water only
! 2009-12-20 gayno - add option to calculate surface fields based on
! the size/shape of field of view.
! 2010-01-28 derber - move calculation of sun glint angle to setuprad
! 2010-02-25 collard - move where nread is calculated to before thinning
! step (more consistent with other obs).
! 2010-03-22 collard - ensure solar azimuth is in the range 0-360 degrees.
! 2011-04-08 li - (1) use nst_gsi, nstinfo, fac_dtl, fac_tsl and add NSST vars
! (2) get zob, tz_tr (call skindepth and cal_tztr)
! (3) interpolate NSST Variables to Obs. location (call deter_nst)
! (4) add more elements (nstinfo) in data array
! 2011-08-01 lueken - move deter_sfc_amsre_low to new module deter_sfc_mod,
! fix indentation
! 2011-09-13 gayno - improve error handling for FOV-based sfc calculation
! (isfcalc=1)
! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS)
! 2013-02-13 eliu - bug fix for solar zenith calculation
! 2012-03-05 akella - nst now controlled via coupler
! 2015-02-23 Rancic/Thomas - add thin4d to time window logical
! 2015-10-01 guo - consolidate use of ob location (in deg)
!
! input argument list:
! mype - mpi task id
! val_amsre- weighting factor applied to super obs
! ithin - flag to thin data
! isfcalc - flag to specify method to calculate sfc fields within FOV
! when set to one, account for size/shape of FOV. otherwise
! use bilinear interpolation.
! rmesh - thinning mesh size (km)
! gstime - analysis time in minutes from reference date
! infile - unit from which to read BUFR data
! lunout - unit to which to write data for further processing
! obstype - observation type to process
! twind - input group time window (hours)
! sis - satellite/instrument/sensor indicator
! mype_root - "root" task for sub-communicator
! mype_sub - mpi task id within sub-communicator
! npe_sub - number of data read tasks
! mpi_comm_sub - sub-communicator for data read
! nrec_start - first subset with useful information
!
! output argument list:
! nread - number of BUFR AQUA observations read
! ndata - number of BUFR AQUA profiles retained for further processing
! nodata - number of BUFR AQUA observations retained for further processing
! nobs - array of observations on each subdomain for each processor
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_double,i_kind
use satthin, only: super_val,itxmax,makegrids,map2tgrid,destroygrids, &
checkob,finalcheck,score_crit
use radinfo, only: iuse_rad,nusis,jpch_rad
use gridmod, only: diagnostic_reg,regional,nlat,nlon,rlats,rlons,&
tll2xy
use constants, only: deg2rad,zero,one,three,r60inv
use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,thin4d
use calc_fov_conical, only: instrument_init
use deter_sfc_mod, only: deter_sfc_fov,deter_sfc,deter_sfc_amsre_low
use gsi_nstcouplermod, only: nst_gsi,nstinfo
use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter
use mpimod, only: npe
! use radiance_mod, only: rad_obs_type
implicit none
! Input variables
character(len=*) ,intent(in ) :: infile
character(len=*) ,intent(in ) :: obstype
integer(i_kind) ,intent(in ) :: mype,nrec_start
integer(i_kind) ,intent(inout) :: isfcalc
integer(i_kind) ,intent(in ) :: ithin
integer(i_kind) ,intent(in ) :: lunout
real(r_kind) ,intent(inout) :: val_amsre
real(r_kind) ,intent(in ) :: gstime,twind
real(r_kind) ,intent(in ) :: rmesh
character(len=20),intent(in ) :: sis
integer(i_kind) ,intent(in ) :: mype_root
integer(i_kind) ,intent(in ) :: mype_sub
integer(i_kind) ,intent(in ) :: npe_sub
integer(i_kind) ,intent(in ) :: mpi_comm_sub
logical ,intent(in ) :: dval_use
! Output variables
integer(i_kind) ,intent(inout) :: nread
integer(i_kind),dimension(npe) ,intent(inout) :: nobs
integer(i_kind) ,intent(inout) :: ndata,nodata
! Number of channels for sensors in BUFR
integer(i_kind),parameter :: N_AMSRCH = 12
! integer(i_kind),parameter :: N_MAXCH = 20
integer(i_kind) :: said, AQUA_SAID = 784 !WMO satellite identifier
integer(i_kind) :: siid, AMSRE_SIID = 345 !WMO instrument identifier
! BUFR file sequencial number
character(len=8) :: subset
character(len=4) :: senname
integer(i_kind) :: lnbufr = 10
integer(i_kind) :: nchanl
integer(i_kind) :: iret,isflg,idomsfc(1)
! Work variables for time
integer(i_kind) :: idate
integer(i_kind) :: idate5(5)
integer(i_kind) :: nmind
real(r_kind) :: sstime, tdiff, t4dv
! Other work variables
logical :: outside,iuse,assim
integer(i_kind) :: nreal, kidsat
integer(i_kind) :: itx, k, nele, itt
integer(i_kind) :: ifov, ilat, ilon
integer(i_kind) :: i, l, n
integer(i_kind),dimension(n_amsrch) :: kchamsre
real(r_kind) :: sfcr
real(r_kind) :: dlon, dlat
real(r_kind) :: dlon_earth,dlat_earth
real(r_kind) :: timedif, pred, crit1, dist1
real(r_kind),allocatable,dimension(:,:):: data_all
integer(i_kind),allocatable,dimension(:)::nrec
integer(i_kind):: irec,next
real(r_kind),dimension(0:3):: sfcpct
real(r_kind),dimension(0:4):: rlndsea
real(r_kind),dimension(0:3):: ts
real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10
real(r_kind) :: zob,tref,dtw,dtc,tz_tr
character(len=7),parameter:: fov_flag="conical"
character(len=3) :: fov_satid
integer(i_kind) :: ichan, instr, idum
logical :: valid
real(r_kind) :: clath_sun_glint_calc , clonh_sun_glint_calc
real(r_kind) :: date5_4_sun_glint_calc
real(r_kind) :: expansion, dlat_earth_deg, dlon_earth_deg
! Set standard parameters
logical :: amsre_low
logical :: amsre_mid
logical :: amsre_hig
integer(i_kind) ntest,maxinfo
integer(i_kind) :: nscan,iskip,kskip,kch,kchanl
! real(r_kind),parameter :: TEN = 10._r_kind
! real(r_kind),parameter :: R45 = 45._r_kind
real(r_kind),parameter :: R90 = 90._r_kind
! real(r_kind),parameter :: R180 = 180._r_kind
real(r_kind),parameter :: R360 = 360._r_kind
real(r_kind),parameter :: tbmin = 70._r_kind
real(r_kind),parameter :: tbmax = 330._r_kind !tbmax is larger than same as ssmiqc
real(r_kind) disterrmax
real(r_kind),dimension(N_AMSRCH) :: tbob_org
real(r_kind) :: clath, clonh, fovn, saza, soza
real(r_kind) :: flgch !used for thinning priority range:1-36
! AMSR-E-bufr
! BUFR format for AQUASPOT
! integer(i_kind),parameter :: N_AQUASPOT_LIST = 25
! BUFR format for AMSRSPOT
integer(i_kind),parameter :: N_AMSRSPOT_LIST = 12
! BUFR format for AMSRCHAN
integer(i_kind),parameter :: N_AMSRCHAN_LIST = 4
! Variables for BUFR IO
real(r_double),dimension(3):: aquaspot_d
real(r_double),dimension(12):: amsrspot_d
real(r_double),dimension(4,12):: amsrchan_d
real(r_double),dimension(2,5):: amsrdice_latlon
real(r_double),dimension(20) :: amsrdice_tmbr
! ---- sun glint ----
integer(i_kind) doy,mlen(12),mday(12),mon,m
real(r_kind) bearaz,sun_zenith,sun_azimuth
data mlen/31,28,31,30,31,30, &
31,31,30,31,30,31/
! Orbit
! logical :: remove_ovlporbit = .true. !looks like AMSRE overlap problem is not as bad as SSM/I 10/14/04 kozo
integer(i_kind) :: orbit, old_orbit, iorbit, ireadsb, ireadmg
real(r_kind) :: saz
! data selection
maxinfo = 31
! Initialize variables
ilon = 3
ilat = 4
if (nst_gsi > 0 ) then
call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype
endif
m = 0
do mon=1,12
mday(mon) = m
m = m + mlen(mon)
end do
disterrmax=zero
ntest = 0
if(dval_use) maxinfo = maxinfo + 2
nreal = maxinfo+nstinfo
ndata = 0
nodata = 0
amsre_low= obstype == 'amsre_low'
amsre_mid= obstype == 'amsre_mid'
amsre_hig= obstype == 'amsre_hig'
orbit = -1
old_orbit=-1
iorbit = 0
sstime = zero
if(amsre_low)then
kchanl=4
kchamsre(1:4)=(/1,2,3,4/)
else if(amsre_mid) then
kchanl=6
kchamsre(1:6)=(/5,6,7,8,9,10/)
else if(amsre_hig)then
kchanl=2
kchamsre(1:2)=(/11,12/)
end if
if(amsre_low .or. amsre_mid .or. amsre_hig)then
senname = 'AMSR'
nchanl = N_AMSRCH
nscan = 196 !for low frequency ch
! nscan = 392 !for 89.0GHz ch
kidsat = 549
rlndsea(0) = zero
rlndsea(1) = 15._r_kind
rlndsea(2) = 10._r_kind
rlndsea(3) = 15._r_kind
rlndsea(4) = 100._r_kind
endif
! If all channels of a given sensor are set to monitor or not
! assimilate mode (iuse_rad<1), reset relative weight to zero.
! We do not want such observations affecting the relative
! weighting between observations within a given thinning group.
assim=.false.
search: do i=1,jpch_rad
if ((nusis(i)==sis) .and. (iuse_rad(i)>0)) then
assim=.true.
exit search
endif
end do search
if (.not.assim) val_amsre=zero
! note, fov-based surface code does not have equations for amsr-e, but
! the fov size/shape is similar to ssmi/s. so call
! fov code using f16 specs.
if(isfcalc==1)then
instr=26
fov_satid='f16'
idum = -999 ! dummy variable for fov number. not used for conical instr.
if(amsre_hig)then
ichan=18
expansion=1.5_r_kind
elseif(amsre_mid)then
ichan=3
expansion=2.9_r_kind
elseif(amsre_low)then
ichan=12
expansion=2.9_r_kind
endif
call instrument_init(instr, fov_satid, expansion, valid)
if (.not. valid) then
if (assim) then
write(6,*)'READ_AMSRE: ***ERROR*** IN SETUP OF FOV-SFC CODE. STOP'
call stop2(71)
else
isfcalc = 0
write(6,*)'READ_AMSRE: ***ERROR*** IN SETUP OF FOV-SFC CODE'
endif
endif
endif
! Make thinning grids
call makegrids(rmesh,ithin)
! Open BUFR file
open(lnbufr,file=infile,form='unformatted')
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
! Allocate local array to contain observation information
nele=nreal+nchanl
allocate(data_all(nele,itxmax),nrec(itxmax))
! Big loop to read data file
next=0
nrec=999999
irec=0
read_msg: do while(ireadmg(lnbufr,subset,idate)>=0)
irec=irec+1
if(irec < nrec_start) cycle read_msg
next=next+1
if(next == npe_sub)next=0
if(next /= mype_sub)cycle
read_loop: do while (ireadsb(lnbufr)==0)
! Retrieve bufr 1/4 :get aquaspot (said,orbn,soza)
call ufbint(lnbufr,aquaspot_d,3,1,iret,'SAID ORBN SOZA')
said = nint(aquaspot_d(1))
if(said /= AQUA_SAID) cycle read_loop
! Retrieve bufr 2/4 :get amsrspot (siid,ymdhs,lat,lon)
call ufbrep(lnbufr,amsrspot_d,N_AMSRSPOT_LIST,1,iret, &
'SIID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ FOVN')
siid = nint(amsrspot_d(1))
if(siid /= AMSRE_SIID) cycle read_loop
! Check obs time
idate5(1) = amsrspot_d(02)! year
idate5(2) = amsrspot_d(03)! month
idate5(3) = amsrspot_d(04)! day
idate5(4) = amsrspot_d(05)! hour
idate5(5) = amsrspot_d(06)! min
if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
idate5(2) < 1 .or. idate5(2) > 12 .or. &
idate5(3) < 1 .or. idate5(3) > 31 .or. &
idate5(4) < 0 .or. idate5(4) > 24 .or. &
idate5(5) < 0 .or. idate5(5) > 60 )then
write(6,*)'READ_AMSRE: ### ERROR IN READING BUFR DATA:', &
' STRANGE OBS TIME (YMDHM):', idate5(1:5)
cycle read_loop
endif
call w3fs21(idate5,nmind)
t4dv = (real((nmind-iwinbgn),r_kind) + amsrspot_d(7)*r60inv)*r60inv ! add in seconds
sstime = real(nmind,r_kind) + amsrspot_d(7)*r60inv ! add in seconds
tdiff = (sstime - gstime)*r60inv
if (l4dvar.or.l4densvar) then
if (t4dv<zero .OR. t4dv>winlen) cycle read_loop
else
if (abs(tdiff)>twind) cycle read_loop
endif
if (thin4d) then
timedif = zero
else
timedif = 6.0_r_kind*abs(tdiff) ! range: 0 to 18
endif
! --- Check observing position -----
if(amsre_low .or. amsre_mid) then
clath= amsrspot_d(08)
clonh= amsrspot_d(09)
else if(amsre_hig)then
call ufbrep(lnbufr,amsrdice_latlon,2, 5,iret,'CLATH CLONH')
clath = amsrdice_latlon(1,4)
clonh = amsrdice_latlon(2,4)
endif
if( abs(clath) > R90 .or. abs(clonh) > R360 .or. &
( abs(clath) == R90 .and. clonh /= ZERO) ) then
! write(6,*)'READ_AMSRE: ### ERROR IN READING BUFR DATA:',&
! ' STRANGE OBS POINT (LAT,LON):', clath, clonh
cycle read_loop
endif
! Pick up every three scene 3,6,9,,,,195 (num=65)
! because of too small scene size and too many scene numbers
! (low-freq ch FOV are overlapped side by side)
fovn = amsrspot_d(12)
! Set position in a given region
if(clonh >= R360)then
clonh = clonh - R360
else if(clonh < ZERO)then
clonh = clonh + R360
endif
! If regional, map obs lat,lon to rotated grid.
dlat_earth = clath * deg2rad
dlon_earth = clonh * deg2rad
dlat_earth_deg = clath
dlon_earth_deg = clonh
if(regional)then
! Convert to rotated coordinate. dlon centered on 180 (pi),
! so always positive for limited area
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
! Check to see if in domain. outside=.true. if dlon_earth,
! dlat_earth outside domain, =.false. if inside
if(outside) cycle read_loop
! Global case
else
dlat=dlat_earth
dlon=dlon_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
! Sum number of read obs before thinning step. Note that this number will contain
! some observations that may be rejected later due to bad BTs.
nread=nread+kchanl
crit1 = 0.01_r_kind+timedif
call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)
if (.not.iuse) cycle read_loop
! QC: "Score" observation. We use this information to identify "best" obs
! Locate the observation on the analysis grid. Get sst and land/sea/ice
! mask
! isflg - surface flag
! 0 sea
! 1 land
! 2 sea ice
! 3 snow
! 4 mixed
! When isfcalc is one, calculate surface fields based on size/shape of fov.
! Otherwise, use bilinear interpolation.
if (isfcalc==1)then
call deter_sfc_fov(fov_flag,idum,instr,ichan,real(amsrspot_d(11),r_kind),dlat_earth_deg,&
dlon_earth_deg,expansion,t4dv,isflg,idomsfc(1), &
sfcpct,vfr,sty,vty,stp,sm,ff10,sfcr,zz,sn,ts,tsavg)
else
call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc(1),sfcpct, &
ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr)
if (amsre_low) then
call deter_sfc_amsre_low(dlat_earth,dlon_earth,isflg,sfcpct)
endif
if(isflg/=0) cycle read_loop ! use data over water only
endif
crit1 = crit1 +rlndsea(isflg)
call checkob(dist1,crit1,itx,iuse)
if(.not. iuse)cycle read_loop
! Retrieve bufr 3/4 : get amsrchan (chnm,tbb)
call ufbrep(lnbufr,amsrchan_d,N_AMSRCHAN_LIST,12,iret,'CHNM LOGRCW ACQF TMBR')
! Retrieve bufr 4/4 : get amsrfovn (fovn)
call ufbrep(lnbufr,amsrdice_tmbr,1,20,iret,'TMBR')
tbob_org(1)=amsrchan_d(4,2)
tbob_org(2)=amsrchan_d(4,1)
tbob_org(3)=amsrchan_d(4,4)
tbob_org(4)=amsrchan_d(4,3)
tbob_org(5)=amsrchan_d(4,6)
tbob_org(6)=amsrchan_d(4,5)
tbob_org(7)=amsrchan_d(4,8)
tbob_org(8)=amsrchan_d(4,7)
tbob_org(9)=amsrchan_d(4,10)
tbob_org(10)=amsrchan_d(4,9)
tbob_org(11)=amsrdice_tmbr(18)
tbob_org(12)=amsrdice_tmbr(17)
! Set obs information
iskip = 0
do l=1,nchanl
if(tbob_org(l)<tbmin .or. tbob_org(l)>tbmax)then
iskip = iskip + 1
end if
end do
kskip = 0
do l=1,kchanl
kch=kchamsre(l)
if(tbob_org(kch)<tbmin .or. tbob_org(kch)>tbmax)then
kskip = kskip + 1
endif
end do
if(kskip == kchanl .or. iskip == nchanl) cycle read_loop
flgch=iskip*three !used for thin, range 0 to 36
crit1 = crit1 + flgch
! Set data quality predictor ***NEED TO COME UP WITH THIS***
pred = zero
! Compute "score" for observation. All scores>=0.0. Lowest score is "best"
! Map obs to grids
crit1 = crit1+pred
call finalcheck(dist1,crit1,itx,iuse)
if(.not. iuse)cycle read_loop
soza = aquaspot_d(3) !solar zenith angle
! Check observational info
if( soza < -180._r_kind .or. soza > 180._r_kind )then
write(6,*)'READ_AMSRE: ### ERROR IN READING BUFR DATA:', &
' STRANGE OBS INFO(FOV,SAZA,SOZA):', fovn, saza, soza
cycle read_loop
endif
! -------- Retreive Sun glint angle -----------
doy = mday( int(idate5(2)) ) + int(idate5(3))
if ((mod( int(idate5(1)),4)==0).and.( int(idate5(2)) > 2)) then
doy = doy + 1
end if
ifov = nint(fovn)
bearaz=amsrspot_d(11)-180.0_r_kind
clath_sun_glint_calc = clath
clonh_sun_glint_calc = clonh
if(clonh>180_r_kind) clonh_sun_glint_calc = clonh -360.0_r_kind
! date5_4_sun_glint_calc = idate5(4)
date5_4_sun_glint_calc = &
real(idate5(4),r_kind)+real(idate5(5),r_kind)*r60inv+real(amsrspot_d(7),r_kind)*r60inv*r60inv
call zensun(doy,date5_4_sun_glint_calc,clath_sun_glint_calc,clonh_sun_glint_calc,sun_zenith,sun_azimuth)
if(amsre_low .or. amsre_mid) then
saz = 55.0_r_kind*deg2rad ! satellite zenith angle (rad)
else if (amsre_hig) then
saz = 54.5_r_kind*deg2rad ! satellite zenith angle (rad)
end if
! saz = amsrspot(10,1)*deg2rad ! satellite zenith angle (rad)
! ==> not use this value but fixed values(55.0 deg) 10/12/04
! because BUFR saza value looks strange (raging -3 to 25),
!
! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr
!
if ( nst_gsi > 0 ) then
tref = ts(0)
dtw = zero
dtc = zero
tz_tr = one
if ( sfcpct(0) > zero ) then
call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr)
endif
endif
data_all(1,itx) = 49 ! satellite ID
data_all(2,itx) = t4dv ! time diff (obs - anal) (hours)
data_all(3,itx) = dlon ! grid relative longitude
data_all(4,itx) = dlat ! grid relative latitude
data_all(5,itx) = saz ! satellite zenith angle (rad)
data_all(6,itx) = amsrspot_d(11) ! satellite azimuth angle
data_all(7,itx) = zero ! look angle (rad)
! data_all(8,itx) = ifov ! fov number 1-196
data_all(8,itx) = ifov/3 + 1 ! fov number/3 1-65 !kozo
data_all(9,itx) = sun_zenith ! solar zenith angle (deg)
data_all(10,itx)= sun_azimuth ! solar azimuth angle (deg)
data_all(11,itx) = sfcpct(0) ! sea percentage of
data_all(12,itx) = sfcpct(1) ! land percentage
data_all(13,itx) = sfcpct(2) ! sea ice percentage
data_all(14,itx) = sfcpct(3) ! snow percentage
data_all(15,itx)= ts(0) ! ocean skin temperature
data_all(16,itx)= ts(1) ! land skin temperature
data_all(17,itx)= ts(2) ! ice skin temperature
data_all(18,itx)= ts(3) ! snow skin temperature
data_all(19,itx)= tsavg ! average skin temperature
data_all(20,itx)= vty ! vegetation type
data_all(21,itx)= vfr ! vegetation fraction
data_all(22,itx)= sty ! soil type
data_all(23,itx)= stp ! soil temperature
data_all(24,itx)= sm ! soil moisture
data_all(25,itx)= sn ! snow depth
data_all(26,itx)= zz ! surface height
data_all(27,itx)= idomsfc(1) + 0.001_r_kind ! dominate surface type
data_all(28,itx)= sfcr ! surface roughness
data_all(29,itx)= ff10 ! ten meter wind factor
data_all(30,itx)= dlon_earth_deg ! earth relative longitude (degrees)
data_all(31,itx)= dlat_earth_deg ! earth relative latitude (degrees)
if(dval_use)then
data_all(32,itx)= val_amsre
data_all(33,itx)= itt
end if
if ( nst_gsi > 0 ) then
data_all(maxinfo+1,itx) = tref ! foundation temperature
data_all(maxinfo+2,itx) = dtw ! dt_warm at zob
data_all(maxinfo+3,itx) = dtc ! dt_cool at zob
data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr)
endif
do l=1,nchanl
data_all(l+nreal,itx) = tbob_org(l)
end do
nrec(itx)=irec
enddo read_loop
enddo read_msg
call closbf(lnbufr)
! If multiple tasks read input bufr file, allow each tasks to write out
! information it retained and then let single task merge files together
call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
nele,itxmax,nread,ndata,data_all,score_crit,nrec)
! Allow single task to check for bad obs, update superobs sum,
! and write out data to scratch file for further processing.
if (mype_sub==mype_root.and.ndata>0) then
! Identify "bad" observation (unreasonable brightness temperatures).
! Update superobs sum according to observation location
do n=1,ndata
do i=1,nchanl
if(data_all(i+nreal,n) > tbmin .and. &
data_all(i+nreal,n) < tbmax)nodata=nodata+1
end do
end do
if(dval_use .and. assim)then
do n=1,ndata
itt=nint(data_all(33,n))
super_val(itt)=super_val(itt)+val_amsre
end do
end if
! Write final set of "best" observations to output file
call count_obs(ndata,nele,ilat,ilon,data_all,nobs)
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata)
endif
deallocate(data_all,nrec) ! Deallocate data arrays
call destroygrids ! Deallocate satthin arrays
if(diagnostic_reg.and.ntest>0 .and. mype_sub==mype_root) &
write(6,*)'READ_AMSRE: ',&
'mype,ntest,disterrmax=',mype,ntest,disterrmax
return
end subroutine read_amsre
subroutine zensun(day,time,lat,lon,sun_zenith,sun_azimuth)
!$$$ subprogram documentation block
! . . . .
! subprogram: zensun make sun zenith and sun azimuth angle
!
! prgmmr: Paul Ricchiazzi org: Earth Space Research Group,UCSB date: 1992-10-23
!
! abstract:
! Compute solar position information as a function of
! geographic coordinates, date and time.
!
! program history log:
! 2005-10-21 kazumori - reformatted for GSI
! 2013-02-13 eliu - bug fix for solar zenith calculation
! 2014-06-04 sienkiewicz - take out doy**3 and replace with simple
! linear interpolation of tables
!
! input argument list:
! day - Julian day (positive scalar or vector)
! (spring equinox = 80)
! (summer solstice= 171)
! (fall equinox = 266)
! (winter solstice= 356)
! time - Universal Time in hours (scalar or vector)
! lat - geographic latitude of point on earth's surface (degrees)
! lon - geographic longitude of point on earth's surface (degrees)
!
! output argument list:
! sun_zenith - solar zenith angle
! sun_azimuth - solar azimuth angle
!
! comments:
!
!
! PROCEDURE:
!
! 1. Calculate the subsolar point latitude and longitude, based on
! DAY and TIME. Since each year is 365.25 days long the exact
! value of the declination angle changes from year to year. For
! precise values consult THE AMERICAN EPHEMERIS AND NAUTICAL
! ALMANAC published yearly by the U.S. govt. printing office. The
! subsolar coordinates used in this code were provided by a
! program written by Jeff Dozier.
!
! 2. Given the subsolar latitude and longitude, spherical geometry is
! used to find the solar zenith, azimuth and flux multiplier.
!
! eqt = equation of time (minutes) ! solar longitude correction = -15*eqt
! dec = declination angle (degrees) = solar latitude
!
! LOWTRAN v7 data (25 points)
! The LOWTRAN solar position data is characterized by only 25 points.
! This should predict the subsolar angles within one degree. For
! increased accuracy add more data points.
!
!nday=[ 1., 9., 21., 32., 44., 60., 91., 121., 141., 152.,$
! 160., 172., 182., 190., 202., 213., 244., 274., 305., 309.,$
! 325., 335., 343., 355., 366.]
!
!eqt=[ -3.23, -6.83,-11.17,-13.57,-14.33,-12.63, -4.2, 2.83, 3.57, 2.45,$
! 1.10, -1.42, -3.52, -4.93, -6.25, -6.28,-0.25, 10.02, 16.35, 16.38,$
! 14.3, 11.27, 8.02, 2.32, -3.23]
!
!dec=[-23.07,-22.22,-20.08,-17.32,-13.62, -7.88, 4.23, 14.83, 20.03, 21.95,$
! 22.87, 23.45, 23.17, 22.47, 20.63, 18.23, 8.58, -2.88,-14.18,-15.45,$
! -19.75,-21.68,-22.75,-23.43,-23.07]
!
! Analemma information from Jeff Dozier
! This data is characterized by 74 points
!
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
!
use kinds, only: r_single,r_kind,i_kind
use constants, only: deg2rad,rad2deg,one,r60inv,zero
implicit none
integer(i_kind), intent(in ) :: day
real(r_kind) , intent(in ) :: time,lat,lon
real(r_kind) , intent( out) :: sun_zenith,sun_azimuth
integer(i_kind) di
real(r_kind) ut,noon
real(r_kind) y(5),y2(5),x(2,5),Tx(5,2),xTx(2,2),aTx(5,2),det
real(r_kind) tt,eqtime,decang,latsun,lonsun
real(r_kind) nday(74),eqt(74),dec(74)
real(r_kind) beta(2), beta2(2), a(2,2)
real(r_kind) t0,t1,p0,p1,zz,xx,yy
data nday/one, 6.0_r_kind, 11.0_r_kind, 16.0_r_kind, 21.0_r_kind, 26.0_r_kind, &
31.0_r_kind, 36.0_r_kind, 41.0_r_kind, 46.0_r_kind,&
51.0_r_kind, 56.0_r_kind, 61.0_r_kind, 66.0_r_kind, 71.0_r_kind, &
76.0_r_kind, 81.0_r_kind, 86.0_r_kind, 91.0_r_kind, 96.0_r_kind,&
101.0_r_kind, 106.0_r_kind, 111.0_r_kind, 116.0_r_kind, 121.0_r_kind, &
126.0_r_kind, 131.0_r_kind, 136.0_r_kind, 141.0_r_kind, 146.0_r_kind,&
151.0_r_kind, 156.0_r_kind, 161.0_r_kind, 166.0_r_kind, 171.0_r_kind, &
176.0_r_kind, 181.0_r_kind, 186.0_r_kind, 191.0_r_kind, 196.0_r_kind,&
201.0_r_kind, 206.0_r_kind, 211.0_r_kind, 216.0_r_kind, 221.0_r_kind, &
226.0_r_kind, 231.0_r_kind, 236.0_r_kind, 241.0_r_kind, 246.0_r_kind,&
251.0_r_kind, 256.0_r_kind, 261.0_r_kind, 266.0_r_kind, 271.0_r_kind, &
276.0_r_kind, 281.0_r_kind, 286.0_r_kind, 291.0_r_kind, 296.0_r_kind,&
301.0_r_kind, 306.0_r_kind, 311.0_r_kind, 316.0_r_kind, 321.0_r_kind, &
326.0_r_kind, 331.0_r_kind, 336.0_r_kind, 341.0_r_kind, 346.0_r_kind,&
351.0_r_kind, 356.0_r_kind, 361.0_r_kind, 366.0_r_kind/
data eqt/ -3.23_r_kind, -5.49_r_kind, -7.60_r_kind, -9.48_r_kind,-11.09_r_kind,&
-12.39_r_kind,-13.34_r_kind,-13.95_r_kind,-14.23_r_kind,-14.19_r_kind,&
-13.85_r_kind,-13.22_r_kind,-12.35_r_kind,-11.26_r_kind,-10.01_r_kind, &
-8.64_r_kind, -7.18_r_kind, -5.67_r_kind, -4.16_r_kind, -2.69_r_kind,&
-1.29_r_kind, -0.02_r_kind, 1.10_r_kind, 2.05_r_kind, 2.80_r_kind, &
3.33_r_kind, 3.63_r_kind, 3.68_r_kind, 3.49_r_kind, 3.09_r_kind,&
2.48_r_kind, 1.71_r_kind, 0.79_r_kind, -0.24_r_kind, -1.33_r_kind, &
-2.41_r_kind, -3.45_r_kind, -4.39_r_kind, -5.20_r_kind, -5.84_r_kind,&
-6.28_r_kind, -6.49_r_kind, -6.44_r_kind, -6.15_r_kind, -5.60_r_kind, &
-4.82_r_kind, -3.81_r_kind, -2.60_r_kind, -1.19_r_kind, 0.36_r_kind,&
2.03_r_kind, 3.76_r_kind, 5.54_r_kind, 7.31_r_kind, 9.04_r_kind, &
10.69_r_kind, 12.20_r_kind, 13.53_r_kind, 14.65_r_kind, 15.52_r_kind,&
16.12_r_kind, 16.41_r_kind, 16.36_r_kind, 15.95_r_kind, 15.19_r_kind, &
14.09_r_kind, 12.67_r_kind, 10.93_r_kind, 8.93_r_kind, 6.70_r_kind,&
4.32_r_kind, 1.86_r_kind, -0.62_r_kind, -3.23_r_kind/
data dec/ -23.06_r_kind,-22.57_r_kind,-21.91_r_kind,-21.06_r_kind,-20.05_r_kind,&
-18.88_r_kind,-17.57_r_kind,-16.13_r_kind,-14.57_r_kind,-12.91_r_kind,&
-11.16_r_kind, -9.34_r_kind, -7.46_r_kind, -5.54_r_kind, -3.59_r_kind,&
-1.62_r_kind, 0.36_r_kind, 2.33_r_kind, 4.28_r_kind, 6.19_r_kind,&
8.06_r_kind, 9.88_r_kind, 11.62_r_kind, 13.29_r_kind, 14.87_r_kind,&
16.34_r_kind, 17.70_r_kind, 18.94_r_kind, 20.04_r_kind, 21.00_r_kind,&
21.81_r_kind, 22.47_r_kind, 22.95_r_kind, 23.28_r_kind, 23.43_r_kind,&
23.40_r_kind, 23.21_r_kind, 22.85_r_kind, 22.32_r_kind, 21.63_r_kind,&
20.79_r_kind, 19.80_r_kind, 18.67_r_kind, 17.42_r_kind, 16.05_r_kind,&
14.57_r_kind, 13.00_r_kind, 11.33_r_kind, 9.60_r_kind, 7.80_r_kind,&
5.95_r_kind, 4.06_r_kind, 2.13_r_kind, 0.19_r_kind, -1.75_r_kind,&
-3.69_r_kind, -5.62_r_kind, -7.51_r_kind, -9.36_r_kind,-11.16_r_kind,&
-12.88_r_kind,-14.53_r_kind,-16.07_r_kind,-17.50_r_kind,-18.81_r_kind,&
-19.98_r_kind,-20.99_r_kind,-21.85_r_kind,-22.52_r_kind,-23.02_r_kind,&
-23.33_r_kind,-23.44_r_kind,-23.35_r_kind,-23.06_r_kind/
!
! compute the subsolar coordinates
!
tt= mod(real((int(day)+time/24._r_kind-one)),365.25_r_single) + one ! fractional day number
! with 12am 1jan = 1.
do di = 1, 73
if ((tt >= nday(di)) .and. (tt <= nday(di+1))) exit
end do
!============== Perform a least squares regression on doy**3 ==============
x(1,:) = one
if ((di >= 3) .and. (di <= 72)) then
y(:) = eqt(di-2:di+2)
y2(:) = dec(di-2:di+2)
x(2,:) = nday(di-2:di+2)**3
end if
if (di == 2) then
y(1) = eqt(73)
y(2:5) = eqt(di-1:di+2)
y2(1) = dec(73)
y2(2:5) = dec(di-1:di+2)
x(2,1) = nday(73)**3
x(2,2:5) = (365._r_kind+nday(di-1:di+2))**3
end if
if (di == 1) then
y(1:2) = eqt(72:73)
y(3:5) = eqt(di:di+2)
y2(1:2) = dec(72:73)
y2(3:5) = dec(di:di+2)
x(2,1:2) = nday(72:73)**3
x(2,3:5) = (365._r_kind+nday(di:di+2))**3
end if
if (di == 73) then
y(1:4) = eqt(di-2:di+1)
y(5) = eqt(2)
y2(1:4) = dec(di-2:di+1)
y2(5) = dec(2)
x(2,1:4) = nday(di-2:di+1)**3
x(2,5) = (365._r_kind+nday(2))**3
end if
if (di == 74) then
y(1:3) = eqt(di-2:di)
y(4:5) = eqt(2:3)
y2(1:3) = dec(di-2:di)
y2(4:5) = dec(2:3)
x(2,1:3) = nday(di-2:di)**3
x(2,4:5) = (365._r_kind+nday(2:3))**3
end if
! Tx = transpose(x)
Tx(1:5,1)=x(1,1:5)
Tx(1:5,2)=x(2,1:5)
! xTx = MATMUL(x,Tx)
xTx(1,1)=x(1,1)*Tx(1,1)+x(1,2)*Tx(2,1)+x(1,3)*Tx(3,1)+x(1,4)*Tx(4,1)+x(1,5)*Tx(5,1)
xTx(1,2)=x(1,1)*Tx(1,2)+x(1,2)*Tx(2,2)+x(1,3)*Tx(3,2)+x(1,4)*Tx(4,2)+x(1,5)*Tx(5,2)
xTx(2,1)=x(2,1)*Tx(1,1)+x(2,2)*Tx(2,1)+x(2,3)*Tx(3,1)+x(2,4)*Tx(4,1)+x(2,5)*Tx(5,1)
xTx(2,2)=x(2,1)*Tx(1,2)+x(2,2)*Tx(2,2)+x(2,3)*Tx(3,2)+x(2,4)*Tx(4,2)+x(2,5)*Tx(5,2)
det = xTx(1,1)*xTx(2,2) - xTx(1,2)*xTx(2,1)
a(1,1) = xTx(2,2)/det
a(1,2) = -xTx(1,2)/det
a(2,1) = -xTx(2,1)/det
a(2,2) = xTx(1,1)/det
! aTx = MATMUL(Tx,a)
aTx(1,1)=Tx(1,1)*a(1,1)+Tx(1,2)*a(2,1)
aTx(2,1)=Tx(2,1)*a(1,1)+Tx(2,2)*a(2,1)
aTx(3,1)=Tx(3,1)*a(1,1)+Tx(3,2)*a(2,1)
aTx(4,1)=Tx(4,1)*a(1,1)+Tx(4,2)*a(2,1)
aTx(5,1)=Tx(5,1)*a(1,1)+Tx(5,2)*a(2,1)
aTx(1,2)=Tx(1,1)*a(1,2)+Tx(1,2)*a(2,2)
aTx(2,2)=Tx(2,1)*a(1,2)+Tx(2,2)*a(2,2)
aTx(3,2)=Tx(3,1)*a(1,2)+Tx(3,2)*a(2,2)
aTx(4,2)=Tx(4,1)*a(1,2)+Tx(4,2)*a(2,2)
aTx(5,2)=Tx(5,1)*a(1,2)+Tx(5,2)*a(2,2)
! beta = MATMUL(y,aTx)
beta(1) = y(1)*aTx(1,1)+y(2)*aTx(2,1)+y(3)*aTx(3,1)+y(4)*aTx(4,1)+y(5)*aTx(5,1)
beta(2) = y(1)*aTx(1,2)+y(2)*aTx(2,2)+y(3)*aTx(3,2)+y(4)*aTx(4,2)+y(5)*aTx(5,2)
! beta2 = MATMUL(y2,aTx)
beta2(1) = y2(1)*aTx(1,1)+y2(2)*aTx(2,1)+y2(3)*aTx(3,1)+y2(4)*aTx(4,1)+y2(5)*aTx(5,1)
beta2(2) = y2(1)*aTx(1,2)+y2(2)*aTx(2,2)+y2(3)*aTx(3,2)+y2(4)*aTx(4,2)+y2(5)*aTx(5,2)
!============== finished least squares regression on doy**3 ==============
! if ((di < 3) .or. (di > 72)) tt = tt + 365._r_kind
eqtime=(beta(1) + beta(2)*tt**3)*r60inv
decang=beta2(1) + beta2(2)*tt**3
latsun=decang
ut=time
noon=12._r_kind-lon/15._r_kind ! universal time of noon
lonsun=-15._r_kind*(ut-12._r_kind+eqtime)
t0=(90._r_kind-lat)*deg2rad
t1=(90._r_kind-latsun)*deg2rad
p0=lon*deg2rad
p1=lonsun*deg2rad
zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0)
! zz2=sin(t0)*sin(t1)+cos(t0)*cos(t1)*cos(p1-p0)
xx=sin(t1)*sin(p1-p0)
yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
sun_zenith=90_r_kind-acos(zz)*rad2deg
sun_azimuth=atan2(xx,yy)*rad2deg
if (sun_azimuth < zero) sun_azimuth = sun_azimuth + 360.0_r_kind
return
end subroutine zensun