module buddycheck_mod
!$$$ module documentation block
! . . . .
! module: adjust_cloudobs_mod
! prgmmr: carley org: np22 date: 2014-09-12
!
! abstract: Module contains routines performing buddychecks
!
! subroutines included:
! sub buddy_check_t - routine to call from setupt to perform buddy check on
! T innovations
! sub execute_buddy_check - routine which performs actual buddy check operation
! (e.g. is called by buddy_check_t)
!
! functions included:
! gc_dist: - returns the great circle distance (m) between
! two pairs of lat/lon points
!
! attributes:
! language: f90/95/2003
! machine: WCOSS
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_double
implicit none
! set default to private
private
! set subroutines to public
public :: buddy_check_t
public :: execute_buddy_check
public :: gc_dist
contains
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE: buddy_check_t --- Perform buddy check for temperature observations
!
! !INTERFACE:
!
subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse)
! !USES:
use gridmod, only: nsig,regional
use guess_grids, only: nfldsig, hrdifsig,ges_lnprsl,&
geop_hgtl,ges_tsen
use gridmod,only:pt_ll
use constants, only: zero,one,r10
use obsmod, only: bmiss,sfcmodel,time_offset
use m_dtime, only: dtime_setup, dtime_check, dtime_show
use gsi_bundlemod, only : gsi_bundlegetpointer
use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle
use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc, &
aircraft_t_bc_ext
use convinfo, only: ictype
implicit none
! !INPUT PARAMETERS:
real(r_kind) , intent(in ) :: data(nele,nobs) ! data array containing all observations
integer(i_kind) , intent(in ) :: mype ! mpi task id
integer(i_kind) , intent(in ) :: nele ! number of data elements per observation
integer(i_kind) , intent(in ) :: nobs ! number of observations
integer(i_kind) , intent(in ) :: is ! ndat index
logical , intent(in ) :: luse(nobs),muse(nobs)
! !OUTPUT PARAMETERS:
integer(i_kind) , intent( out) :: buddyuse(nobs) !Holds info pertaining to buddycheck.
! 1 = pass
! 0 = buddy check no performed
! -1 = fail
! !INPUT/OUTPUT PARAMETERS:
!
!
! !DESCRIPTION: This routine creates temperature innovations and then performs a simple two-pass
! buddy check on the innovations. An array indicating which observations have
! passed, failed, or were not checked is returned.
!
!
!
! !REVISION HISTORY:
!
! 2014-09-05 Carley - Originator. Initial version based off setupt.f90 to obtain innovations
!
! !REMARKS:
! language: f90/95
! machine: WCOSS
!
! !AUTHOR:
! Carley org: np22 date: 2014-09-05
!
!EOP
!-------------------------------------------------------------------------
! Declare local parameters
real(r_kind),parameter:: r0_001 = 0.001_r_kind
real(r_kind),parameter:: r0_7=0.7_r_kind
real(r_kind),parameter:: r8 = 8.0_r_kind
character(len=*),parameter :: myname='buddy_check_t'
! Declare external calls for code analysis
external:: SFC_WTQ_FWD
external:: get_tlm_tsfc
! A bug with the intel compiler requires that the
! externals for the tintrp* routines be commented out
! in order for GSI to compile in debug mode. (has to do
! with the fact that these are set as externals in setupt,
! which calls this routine).
! external:: tintrp2a1,tintrp2a11
! external:: tintrp31
external:: grdcrd1
external:: stop2
! Declare local variables
real(r_kind) rsig
real(r_kind) psges
real(r_kind) tges
real(r_kind) tob
real(r_kind) dlon,dlat,dtime,dpres,prest
real(r_kind),dimension(nsig):: prsltmp
real(r_kind) tgges,roges
real(r_kind),dimension(nsig):: tvtmp,qtmp,utmp,vtmp,hsges
real(r_kind),dimension(nobs,7):: vals ! innovation, lat, lon, elev., station id, report time usage from read_prepbufr
real(r_kind) u10ges,v10ges,t2ges,q2ges,psges2,f10ges,range,difmax
real(r_kind),dimension(nsig):: prsltmp2
integer(i_kind) mm1
integer(i_kind) itype,msges,iqt,i
integer(i_kind) ier,ilon,ilat,ipres,itob,id,itime,ikx,iqc,iptrb,icat,ipof,ivvlc,idx
integer(i_kind) ier2,iuse,ilate,ilone,ikxx,istnelv,iobshgt,izz,iprvd,isprvd
integer(i_kind) regime
integer(i_kind) idomsfc,iskint,iff10,isfcr
logical sfctype
logical iqtflg
logical:: in_curbin, in_anybin
logical proceed
real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps
real(r_kind),allocatable,dimension(:,:,: ) :: ges_z
real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u
real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v
real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv
real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q
! Check to see if required guess fields are available
call check_vars_(proceed)
if(.not.proceed) return ! not all vars available, simply return
! If require guess vars available, extract from bundle ...
call init_vars_
! index information for data array (see reading routine)
ier=1 ! index of obs error
ilon=2 ! index of grid relative obs location (x)
ilat=3 ! index of grid relative obs location (y)
ipres=4 ! index of pressure
itob=5 ! index of t observation
id=6 ! index of station id
itime=7 ! index of observation time in data array
ikxx=8 ! index of ob type
iqt=9 ! index of flag indicating if moisture ob available
iqc=10 ! index of quality mark
ier2=11 ! index of original-original obs error ratio
iuse=12 ! index of use parameter
idomsfc=13 ! index of dominant surface type
iskint=14 ! index of surface skin temperature
iff10=15 ! index of 10 meter wind factor
isfcr=16 ! index of surface roughness
ilone=17 ! index of longitude (degrees)
ilate=18 ! index of latitude (degrees)
istnelv=19 ! index of station elevation (m)
iobshgt=20 ! index of observation height (m)
izz=21 ! index of surface height
iprvd=22 ! index of observation provider
isprvd=23 ! index of observation subprovider
icat=24 ! index of data level category
if (aircraft_t_bc_pof .or. aircraft_t_bc .or. aircraft_t_bc_ext) then
ipof=25 ! index of data pof
ivvlc=26 ! index of data vertical velocity
idx=27 ! index of tail number
iptrb=28 ! index of t perturbation
else
iptrb=25 ! index of t perturbation
end if
rsig=float(nsig)
mm1=mype+1
!initialize buddyuse to 1, start by assuming all obs are good!
buddyuse=1
call dtime_setup()
do i=1,nobs
dtime=data(itime,i)
call dtime_check(dtime, in_curbin, in_anybin)
if( (.not.in_anybin) .or. (.not.in_curbin) .or. (.not. luse(i)) .or. (.not. muse(i)) ) then
!If outside the time window - don't run the buddy check since we will not be using this ob anyway
! Also do not run if muse or luse is false
buddyuse(i)=0
vals(i,1)=bmiss
vals(i,2)=data(ilate,i)
vals(i,3)=data(ilone,i)
vals(i,4)=data(iobshgt,i)
vals(i,5)=data(iuse,i)
vals(i,6)=data(id,i)
vals(i,7)=dtime-time_offset
cycle
end if
! Convert obs lats and lons to grid coordinates
dlat=data(ilat,i)
dlon=data(ilon,i)
dpres=data(ipres,i)
ikx=nint(data(ikxx,i))
itype=ictype(ikx)
prest=r10*exp(dpres) ! in mb
sfctype=(itype>179.and.itype<190).or.(itype>=192.and.itype<=199)
iqtflg=nint(data(iqt,i)) == 0
! Load observation value
tob=data(itob,i)
! Interpolate log(ps) & log(pres) at mid-layers to obs locations/times
call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,&
mype,nfldsig)
call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
! Put obs pressure in correct units to get grid coord. number
call grdcrd1(dpres,prsltmp(1),nsig,-1)
! Implementation of forward model ----------
if(sfctype.and.sfcmodel) then
tgges=data(iskint,i)
roges=data(isfcr,i)
msges = 0
if(itype == 180 .or. itype == 182 .or. itype == 183) then !sea
msges=0
elseif(itype == 181 .or. itype == 187 .or. itype == 188) then !land
msges=1
endif
call tintrp2a1(ges_tv,tvtmp,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
call tintrp2a1(ges_q,qtmp,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
call tintrp2a1(ges_u,utmp,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
call tintrp2a1(ges_v,vtmp,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
call tintrp2a1(geop_hgtl,hsges,dlat,dlon,dtime,hrdifsig,&
nsig,mype,nfldsig)
psges2 = psges ! keep in cb
prsltmp2 = exp(prsltmp) ! convert from ln p to cb
call SFC_WTQ_FWD (psges2, tgges,&
prsltmp2(1), tvtmp(1), qtmp(1), utmp(1), vtmp(1), &
prsltmp2(2), tvtmp(2), qtmp(2), hsges(1), roges, msges, &
f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg)
tges = t2ges
else
if(iqtflg)then
! Interpolate guess tv to observation location and time
call tintrp31(ges_tv,tges,dlat,dlon,dpres,dtime, &
hrdifsig,mype,nfldsig)
else
! Interpolate guess tsen to observation location and time
call tintrp31(ges_tsen,tges,dlat,dlon,dpres,dtime, &
hrdifsig,mype,nfldsig)
end if
endif
if(sfctype.and.sfcmodel) dpres = one ! place sfc T obs at the model sfc
if (dpres > rsig )then
if( regional .and. prest > pt_ll )then
dpres=rsig
else
!Do no use obs outside the domain
buddyuse(i)=0
vals(i,1)=bmiss
vals(i,2)=data(ilate,i)
vals(i,3)=data(ilone,i)
vals(i,4)=data(iobshgt,i)
vals(i,5)=data(iuse,i)
vals(i,6)=data(id,i)
vals(i,7)=dtime-time_offset
cycle
endif
endif
! Compute innovation and store remaining data
vals(i,1) = tob-tges
vals(i,2)=data(ilate,i)
vals(i,3)=data(ilone,i)
vals(i,4)=data(iobshgt,i)
vals(i,5)=data(iuse,i)
vals(i,6)=data(id,i)
vals(i,7)=dtime-time_offset
! End of loop over observations
end do
! Release memory of local guess arrays
call final_vars_
! - Now call buddy check routine
range = 108000.0_r_kind ! Radius within which we check for an ob's buddies (units are m)
!
! 8/31/2015: adjust the difmax value based on case: 2/29/2015 03Z.
! difmax= 7.0_r_kind ! the final codeMax difference allowed relative to buddies
!for being consistent, difmax=8, but late will be changed.
difmax= 8.0_r_kind ! Max difference allowed relative to buddies
call execute_buddy_check(mype,is,nobs,vals,range,difmax,buddyuse)
! End of routine
return
contains
subroutine check_vars_ (proceed)
logical,intent(inout) :: proceed
integer(i_kind) ivar, istatus
! Check to see if required guess fields are available
call gsi_metguess_get ('var::ps', ivar, istatus )
proceed=ivar>0
call gsi_metguess_get ('var::u' , ivar, istatus )
proceed=proceed.and.ivar>0
call gsi_metguess_get ('var::v' , ivar, istatus )
proceed=proceed.and.ivar>0
call gsi_metguess_get ('var::tv', ivar, istatus )
proceed=proceed.and.ivar>0
call gsi_metguess_get ('var::q', ivar, istatus )
proceed=proceed.and.ivar>0
end subroutine check_vars_
subroutine init_vars_
real(r_kind),dimension(:,: ),pointer:: rank2=>NULL()
real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL()
character(len=5) :: varname
integer(i_kind) ifld, istatus
! If require guess vars available, extract from bundle ...
if(size(gsi_metguess_bundle)==nfldsig) then
! get ps ...
varname='ps'
call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
if (istatus==0) then
if(allocated(ges_z))then
write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
call stop2(999)
endif
allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig))
ges_ps(:,:,1)=rank2
do ifld=2,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
ges_ps(:,:,ifld)=rank2
enddo
else
write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
call stop2(999)
endif
! get u ...
varname='u'
call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
if (istatus==0) then
if(allocated(ges_u))then
write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
call stop2(999)
endif
allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
ges_u(:,:,:,1)=rank3
do ifld=2,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
ges_u(:,:,:,ifld)=rank3
enddo
else
write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
call stop2(999)
endif
! get v ...
varname='v'
call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
if (istatus==0) then
if(allocated(ges_v))then
write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
call stop2(999)
endif
allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
ges_v(:,:,:,1)=rank3
do ifld=2,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
ges_v(:,:,:,ifld)=rank3
enddo
else
write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
call stop2(999)
endif
! get tv ...
varname='tv'
call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
if (istatus==0) then
if(allocated(ges_tv))then
write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
call stop2(999)
endif
allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
ges_tv(:,:,:,1)=rank3
do ifld=2,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
ges_tv(:,:,:,ifld)=rank3
enddo
else
write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
call stop2(999)
endif
! get q ...
varname='q'
call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
if (istatus==0) then
if(allocated(ges_q))then
write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
call stop2(999)
endif
allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
ges_q(:,:,:,1)=rank3
do ifld=2,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
ges_q(:,:,:,ifld)=rank3
enddo
else
write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
call stop2(999)
endif
else
write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',&
nfldsig,size(gsi_metguess_bundle)
call stop2(999)
endif
end subroutine init_vars_
subroutine final_vars_
if(allocated(ges_q )) deallocate(ges_q )
if(allocated(ges_tv)) deallocate(ges_tv)
if(allocated(ges_v )) deallocate(ges_v )
if(allocated(ges_u )) deallocate(ges_u )
if(allocated(ges_ps)) deallocate(ges_ps)
end subroutine final_vars_
end subroutine buddy_check_t
subroutine execute_buddy_check(mype,is,numobs,pevals,range,difmax,pebuddyuse)
! !USES:
use jfunc, only: jiter
use mpimod, only: ierror,mpi_rtype,mpi_itype,mpi_sum,mpi_comm_world
use constants, only: zero, one,one_tenth,r100,tiny_r_kind
use obsmod, only: obs_sub_comm,bmiss,dtype
use qcmod, only: buddydiag_save
implicit none
! !INPUT PARAMETERS:
real(r_kind) , intent(in ) :: pevals(numobs,7) ! data array containing all observations
real(r_kind) , intent(in ) :: range,& ! Radius within which we check for an ob's buddies (units are m)
difmax ! Max difference allowed relative to avg of buddy innovations
integer(i_kind) , intent(in ) :: mype ! mpi task id
integer(i_kind) , intent(in ) :: numobs ! number of observations on this task
integer(i_kind) , intent(in ) :: is ! ndat index
! !OUTPUT PARAMETERS:
integer(i_kind) , intent(inout) :: pebuddyuse(numobs) !Holds info pertaining to buddycheck.
! 1 = pass
! 0 = buddy check no performed
! -1 = fail
! !INPUT/OUTPUT PARAMETERS:
!
!
! !DESCRIPTION: This routine performs a simple two-pass buddy check on provided innovations.
! An array indicating which observations have passed, failed, or were not
! checked is returned.
!
! This buddy check algorithm is based upon that found in the WRFDA
! routine da_buddy_qc, written by Yong-Run Guo 10/10/2008.
!
!
! !REVISION HISTORY:
!
! 2014-09-09 Carley - Originator. Method is based on WRFVAR code from
! var/da/da_tools/da_buddy_qc.inc (Yong-Run Guo)
!
! !REMARKS:
! language: f90/95
! machine: WCOSS
!
! !AUTHOR:
! Carley org: np22 date: 2014-09-09
!
!EOP
!-------------------------------------------------------------------------
! External routines
external:: mpi_allgather
external:: mpi_allreduce
external:: mpi_allgatherv
! Declare local variables
integer(i_kind),parameter :: nvals=7
real(r_kind) :: average,myinnov,mylat,mylon,myelev,sum,&
diff_check_1,diff_check_2,diff_j,err,fact
real(r_kind) :: distance,vdist,mydhr
integer(i_kind) :: i,j,buddy_num,numobs_global,kob,buddy_num_final
integer(i_kind) :: submm1,mynpe,ji,newrank,rem,amount,mydata
integer(i_kind),allocatable,dimension(:) :: buddyuse_global,buddyuse,idisp,ircnt
integer(i_kind),allocatable,dimension(:) :: locdisp,locsendcnt
real(r_kind) :: pevals1d(numobs*nvals)
real(r_kind),allocatable,dimension(:,:) :: vals_global,vals
real(r_kind),allocatable,dimension(:) :: vals1dglob,myvals1d,diff
character(len=*),parameter :: myname='execute_buddy_check'
character(8) :: mype_file,my_jiter
real(r_double) rstation_id
character(8) myid
equivalence(rstation_id,myid)
! Obtain the number of tasks here by using the communicator setup in obs_para.f90
! Note that the only communicator that gets activated here is the one associated with nobs>0
! since only pe's with obs on their subdomains will execute any of the setup* routines within
! setuprhsall.f90
call mpi_comm_size(obs_sub_comm(is), mynpe, ierror)
call mpi_comm_rank(obs_sub_comm(is), newrank, ierror)
submm1=newrank+1
allocate(idisp(mynpe),ircnt(mynpe),locdisp(mynpe),locsendcnt(mynpe))
idisp=0;ircnt=0
if (buddydiag_save) then
! For diagnostic output create a an output ASCII file on each pe to record which
! data pass/fail the buddy check
write (mype_file, "(I4.4)") mype
write (my_jiter, "(I4.4)") jiter
open(1920+mype,file='buddycheck_'//trim(dtype(is))//'_'//trim(mype_file)//'.'//trim(my_jiter),form='formatted',&
status='unknown')
write(1920+mype,'(A)')&
'ID Lat Lon O-F |O-F-AVG(O-Fs)| #Buddies Pass?(1:pass,-1:fail,0s not recorded) time'
end if
call mpi_allgather(numobs,1,mpi_itype,ircnt,1,mpi_itype,obs_sub_comm(is),ierror)
! Now we need to calculate the displacement array here since we can't use the one from the
! global obsveration data array since that one is not ordered in any particular manner
! and nor does it necessarily matter (i.e. ob at n=1 could be in Oregon,
! ob at n=2 may be in Florida, and then ob at n=3 might be back in Oregon)
idisp(1)=0 !first spot will always be zero
if (mynpe >=2) then
do i=2,mynpe
idisp(i)=idisp(i-1)+ircnt(i-1)
end do
end if
! Get the total, global number of obs/innovations we need to consider
call mpi_allreduce(numobs,numobs_global,1,mpi_itype,mpi_sum,obs_sub_comm(is),ierror)
! Allocate the memory for everything and initialize certain fields to zero
allocate(vals_global(numobs_global,nvals),buddyuse_global(numobs_global),vals1dglob(numobs_global*nvals),diff(numobs_global))
vals_global=zero
! Gather all pevals(:,:) on subdomains into vals_global on every task, noting that pevals takes the following
! organization:
! pevals(i,1)=ob-ges
! pevals(i,2)=earth lat
! pevals(i,3)=earth lon
! pevals(i,4)=ob height
! pevals(i,5)=rusage
! pevals(i,6)=rstation_id
! pevals(i,7)=obtime relative to analysis time
!put pevals in a 1d array for easy mpi comms
ji=0
do i=1,numobs
do j=1,nvals
ji=ji+1
pevals1d(ji)=pevals(i,j)
end do
end do
call mpi_allgatherv(pevals1d,numobs*nvals,mpi_rtype,vals1dglob,ircnt*nvals,idisp*nvals,mpi_rtype,obs_sub_comm(is),ierror)
! put vals1dglob to 2d vals_glob so things are more human-readable
ji=0
do i=1,numobs_global
do j=1,nvals
ji=ji+1
vals_global(i,j)=vals1dglob(ji)
end do
end do
! Gather all pebuddyuse params into a global copy (buddyuse_global) that all tasks have
call mpi_allgatherv(pebuddyuse,numobs,mpi_itype,buddyuse_global,ircnt,idisp,mpi_itype,obs_sub_comm(is),ierror)
! Alright - now this may seem silly but we need to scatter obs to new subdomains because there
! is a severe load imbalance issue otherwise (some tasks have 10s of obs and others have
! 1000s or more).
! First figure out the best decomposition given the tasks present *in this routine*
rem=mod(numobs_global,mynpe)
amount=int(numobs_global/mynpe)
do i=1,mynpe
locsendcnt(i)=amount
end do
if (rem /= 0 .and. mynpe>1) locsendcnt(mynpe)=locsendcnt(mynpe)+rem
if (newrank==0 .and. buddydiag_save ) then
write(6,1000)myname,' obs per pe ob subdomain',(locsendcnt(i),i=1,mynpe)
1000 format(2A,1x,1x,8I10,/,(10X,10I10))
end if
locdisp=0
locdisp(1)=0 !first spot will always be zero
if (mynpe >=2) then
do i=2,mynpe
locdisp(i)=locdisp(i-1)+locsendcnt(i-1)
end do
end if
! Now scatterv vals1dglob to workers for better load balance
mydata=locsendcnt(submm1)
allocate(myvals1d(mydata*nvals),vals(mydata,nvals),buddyuse(mydata))
myvals1d=0
call mpi_scatterv(vals1dglob,locsendcnt*nvals,locdisp*nvals,mpi_rtype,myvals1d,mydata*nvals,mpi_rtype,0,obs_sub_comm(is),ierror)
ji=0
do i=1,mydata
do j=1,nvals
ji=ji+1
vals(i,j)=myvals1d(ji)
end do
end do
! Scatter the buddy use array
call mpi_scatterv(buddyuse_global,locsendcnt,locdisp,mpi_itype,buddyuse,mydata,mpi_itype,0,obs_sub_comm(is),ierror)
! Now - finally - start what we came her for originally, the buddy check!
! Begin with main loop through all obs on pe subdomain
main: do i=1,mydata
average = zero
! No buddy check is applied to bad data (rusage => 100)
! Do not run buddy check if buddyuse is already 0 - this ob may
! have luse=.false., be outside the timewindow, etc
if (vals(i,5)>= r100 .or. buddyuse(i)==0 .or. (abs(vals(i,1)-bmiss)<=tiny_r_kind) ) then
buddyuse(i)=0
cycle main
end if
! Now find the distance between all pairs of innovations (numobs_global) and get the number
! within our range, find number of observations within a given distance,
! and store that innovation
buddy_num = 0
myinnov=vals(i,1)
mylat=vals(i,2)
mylon=vals(i,3)
myelev=vals(i,4)
rstation_id=vals(i,6)
mydhr=vals(i,7)
diff=zero
inner: do j = 1, numobs_global
! Do not use buddies who are not eligible
if (buddyuse_global(j)==0 .or. vals_global(j,5) >= r100 .or. &
(abs(vals_global(j,1)-bmiss)<=tiny_r_kind) ) cycle inner
! Calculate horiz. and vertical distance between test innov. and potential buddy
distance=gc_dist(mylat,mylon,vals_global(j,2),vals_global(j,3)) !units are meters (calc. via great circle method)
vdist=abs(myelev-vals_global(j,4)) !units are meters (vertical distance between obs)
! Distance must be greater than 0.1 meters so we don't include the test ob in the buddy list
if (distance <= range .and. distance >= one_tenth .and. vdist <= 200._r_kind) then
buddy_num = buddy_num + 1 ! We have found a buddy, so increment the counter
diff(buddy_num)=vals_global(j,1) ! Save the innovation associated with this buddy
end if
end do inner
! Initialize buddy_num_final to 0
buddy_num_final=0
! Sum all the buddy innovations
sum = zero
do j = 1, buddy_num
sum = sum + diff(j)
end do
! Check to see if there are any buddies, compute the mean innovation.
if (buddy_num > 0) average = sum / buddy_num
! Check if there are any bad obs (bad buddy buddies) among the obs located within the
! the radius surrounding the test ob.
diff_check_1 = difmax*1.25_r_kind
diff_check_2 = difmax
if (buddy_num >= 2) then
kob = buddy_num
do j = 1, buddy_num
diff_j = abs(diff(j)-average)
if ( abs(diff(j)) > diff_check_1 .and. diff_j > diff_check_2 ) then
! Bad obs: Innovation itself: diff(numj) > diff_check_1
! The distance between the innovation and average > diff_check_2
kob = kob - 1
sum = sum - diff(j)
end if
end do
! Set the final number of buddies here - after we've possibly removed some bad ones
buddy_num_final = kob
! We may have removed too many buddies when we were checking for bad buddies.
! Check to see if we still have enough
if (kob > 2) then
average = sum / kob
err = myinnov - average
else
! Not enough buddies, cycle back to the beginning of the main loop
err = zero
buddyuse(i)=0
cycle main
end if
else
! Not enough buddies, cycle back to the beginning of the main loop
buddyuse(i)=0
cycle main
end if
! If the buddy number is ONLY 2, increase the tolerance value:
fact=1.0
if (buddy_num_final == 2) fact=1.2
! Now do the final buddy check
if (abs(err) > fact*difmax) then
buddyuse(i) = -1
! Most obs pass so it can be useful to print every buddy check failure to the screen
if(buddydiag_save) write(6,'(A,A,A,f10.4,A,f10.4,A,f10.4,A,I10,A,f8.4)') 'Station: ',myid,&
' Fails! innov: ',myinnov,' err: ', err,' avg buddy innovs: ',average,' buddy count:', buddy_num_final,' dhr: ',mydhr
else
buddyuse(i)=1
! Most obs pass - so printing every pass to the screen would be overwhelming
if(mod(i,100)==0 .and. buddydiag_save) write(6,'(A,A,A,f10.4,A,f10.4,A,f10.4,A,I10,A,f8.4)') 'Station: ',myid,&
' Passes! innov: ',myinnov,' err: ', err,' avg buddy innovs: ',average, ' buddy count:', buddy_num_final, ' dhr: ',mydhr
end if
if (buddydiag_save) write(1920+mype,'(A,f6.2,3x,f6.2,3x,f6.2,3x,f6.2,3x,I10,4x,I2,3x,f8.4)')&
myid,mylat,mylon,myinnov,err,buddy_num_final,buddyuse(i),mydhr
end do main
if (buddydiag_save) close (1920+mype)
! Begin MPI comms to put buddyuse back on pebuddyuse (AKA the geographic subdomains created in obs_para.f90)
buddyuse_global=0
! Receive count and displacement will be identical to prior scatter send count
call mpi_gatherv(buddyuse,mydata,mpi_itype,buddyuse_global,&
locsendcnt,locdisp,mpi_itype,0,obs_sub_comm(is),ierror)
! mpi_scatterv back to pe subdomains (pebuddyuse)
! Use the ircnt and idisp arrays obtained earlier to put arrays on all procs. They will be the same here.
pebuddyuse=0
call mpi_scatterv(buddyuse_global,ircnt,idisp,mpi_itype,&
pebuddyuse,numobs,mpi_itype,0,obs_sub_comm(is),ierror)
deallocate(vals_global,buddyuse_global,idisp,ircnt,vals1dglob,&
locdisp,locsendcnt,myvals1d,vals,buddyuse,diff)
return
end subroutine execute_buddy_check
real(r_kind) function gc_dist(inlat1,inlon1,inlat2,inlon2)
!$$$ subprogram documentation block
! . . . .
! subprogram: gc_dist
! prgmmr: org: date:
!
! abstract: Return the great circle distance (m) between a two pairs of lat/lon points
!
! program history log:
! 2014-09-09 carley - added subprogram doc block
!
! input argument list:
! lat1,lon1,lat2,lon2 in degrees
!
! output argument list:
! gc_dist
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: rearth,deg2rad,one,two
implicit none
real(r_kind),intent (in) :: inlat1,inlon1,inlat2,inlon2
real(r_kind) :: lat1,lon1,lat2,lon2
real(r_kind) :: dLat,dLon,a,c
lat2=inlat2*deg2rad
lat1=inlat1*deg2rad
lon1=inlon1*deg2rad
lon2=inlon2*deg2rad
dLat = lat2 - lat1
dLon = lon2 - lon1
a = sin(dLat / two) * sin(dLat / two) + cos(lat1) * cos(lat2) * &
sin(dLon / two) * sin(dLon / two)
c = two * atan2(sqrt(a), sqrt(one - a))
gc_dist = rearth * c
end function gc_dist
end module buddycheck_mod