subroutine obs_para(ndata,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: obs_para assign and distribute observations to subdomains
! prgmmr: weiyu yang org: np20 date: 1998-05-27
!
! abstract: Based on observation location and domain decomposition, assign
! and distribute observations to subdomains
!
! program history log:
! 1998-05-27 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-06-15 treadon, reformat documenation
! 2004-07-23 derber - modify to include conventional sst
! 2004-07-28 treadon - add only to module use, add intent in/out
! 2004-11-19 derber - modify to eliminate file and change to use logical
! rather than weights
! 2005-06-14 wu - add OMI total ozone
! 2005-09-08 derber - simplify data set handling
! 2005-10-17 treadon - fix bug in defnition of mype_diaghdr
! 2006-02-03 derber - modify for new obs control and obs count
! 2006-07-28 derber - eliminate counting of data and clean up
! 2008-04-29 safford - rm unused vars and uses
! 2008-06-30 derber - optimize calculation of criterion
! 2008-09-08 lueken - merged ed's changes into q1fy09 code
! 2009-04-21 derber - reformulate to remove communication
! 2008-05-10 meunier - handle for lagrangian data
! 2013-01-26 parrish - attempt fix for bug flagged by WCOSS debug compiler.
! Replace
! "call dislag(.....,nobs_s)"
! with
! "call dislag(.....,nobs_s(mm1))"
! nobs_s is an array in current subroutine, but is a
! scalar inside subroutine dislag.
! 2014-10-03 carley - add creation mpi subcommunicator needed for
! buddy check QC to distinguish among pe subdomains
! with and without obs (only for t obs and twodvar_regional
! at the moment)
! 2015-10-28 guo - added ioid_s(:) to support split-observer GSI with a
! control vector grid different from the guess state
! grid.
!
! input argument list:
! ndata(*,1)- number of prefiles retained for further processing
! ndata(*,2)- number of observations read
! ndata(*,3)- number of observations keep after read
! mype - mpi task number
! ipoint - pointer in array containing information about all obs type to process
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: i_kind
use constants, only: zero
use jfunc, only: factqmin,factqmax
use mpimod, only: npe,mpi_itype,mpi_comm_world,ierror
use obsmod, only: obs_setup,dtype,mype_diaghdr,ndat,nsat1
use obsmod, only: obsfile_all,dplat,nobs_sub,obs_sub_comm
use gridmod, only: twodvar_regional
use qcmod, only: buddycheck_t,buddydiag_save
use gsi_io, only: verbose, print_obs_para
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: mype
integer(i_kind),dimension(ndat,3),intent(in ) :: ndata
! Declare local variables
integer(i_kind) lunout,is,ii
integer(i_kind) mm1
integer(i_kind) ndatax_all,ikey_yes,ikey_no,newprocs,newrank
integer(i_kind),dimension(npe):: ikey,icolor
integer(i_kind) nobs_s
logical print_verbose
print_verbose=.false.
if(verbose)print_verbose=.true.
!
!****************************************************************
! Begin obs_para here
!
! Distribute observations as a function of pe number.
nsat1=0
obs_sub_comm=0
mype_diaghdr = -999
mm1=mype+1
ndatax_all=0
lunout=55
!
! Set number of obs on each pe for each data type
open(lunout,file=obs_setup,form='unformatted')
rewind lunout
do is=1,ndat
if(dtype(is) /= ' ' .and. ndata(is,1) > 0)then
ndatax_all=ndatax_all + ndata(is,1)
if (dtype(is)=='lag') then ! lagrangian data
call dislag(ndata(is,1),mm1,lunout,obsfile_all(is),dtype(is),&
nobs_s)
nsat1(is)= nobs_sub(mm1,is)
else
obproc:do ii=1,npe
if(nobs_sub(ii,is) > 0)then
mype_diaghdr(is) = ii-1
exit obproc
end if
end do obproc
if(nobs_sub(mm1,is) > zero)then ! classical observations
call disobs(ndata(is,1),nobs_sub(mm1,is),mm1,lunout, &
obsfile_all(is),dtype(is))
end if
end if
nsat1(is)= nobs_sub(mm1,is)
if(mm1 == 1 .and. (print_verbose .or. print_obs_para))then
write(6,1000)dtype(is),dplat(is),(nobs_sub(ii,is),ii=1,npe)
1000 format('OBS_PARA: ',2A10,8I10,/,(10X,10I10))
end if
! Simple logic to organize which tasks do and do not have obs.
! Needed for buddy check QC.
if (twodvar_regional .and. dtype(is) == 't' .and. buddycheck_t) then
ikey_yes=0
ikey_no=0
ikey=0
do ii=1,npe
if (nobs_sub(ii,is)>0) then
icolor(ii)=1
ikey(ii)=ikey_yes
ikey_yes=ikey_yes+1
else
icolor(ii)=2
ikey(ii)=ikey_no
ikey_no=ikey_no+1
end if
end do
! With organized colors and keys, now create the new MPI communicator
! which only talks to pe's who have obs on their subdomains. This is
! needed for MPI communication within the setup* routines (e.g. a buddy check).
call mpi_comm_split(mpi_comm_world,icolor(mm1),ikey(mm1),obs_sub_comm(is),ierror)
CALL MPI_COMM_SIZE(obs_sub_comm(is), newprocs, ierror)
CALL MPI_COMM_RANK(obs_sub_comm(is), newrank, ierror)
if (buddydiag_save) write(6,'(A,I3,I10,A,I20,A,I3,A,I3)') 'obs_para: mype/myobs=',&
mype,nobs_sub(mm1,is),'newcomm=',obs_sub_comm(is),'newprocs=', &
newprocs,'newrank=',newrank
end if
end if
end do
close(lunout)
! If there are no obs available, turn off moisture constraint.
! If the user still wants the moisture constraint active when no obs are
! present, comment out the block of code below.
if (ndatax_all == 0) then
factqmin=zero
factqmax=zero
if (mype==0) write(6,*)'OBS_PARA: ***WARNING*** no observations to be ',&
' assimilated. reset factqmin,factqmax=',factqmin,factqmax
endif
return
end subroutine obs_para
subroutine disobs(ndata,nobs,mm1,lunout,obsfile,obstypeall)
!$$$ subprogram documentation block
! . . . .
! subprogram: disobs distribute observations into each pe subdomain
! prgmmr: weiyu yang date: 1998-04-06
!
! abstract: distribute observations into each pe subdomain
!
! program history log:
! 1998-04-06 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-06-15 treadon, reformat documenation
! 2004-07-28 treadon - add only to module use, add intent in/out
! 2004-11-19 derber - change to eliminate additional file and use logical
! rather than weights
! 2005-10-17 treadon - remove conversion to grid-relative cooridnates
! 2006-04-06 middlecoff - changed lunin from 15 to 11
! 2008-04-29 safford - rm unused vars and uses
! 2008-09-08 lueken - merged ed's changes into q1fy09 code
! 2009-04-21 derber - reformulate to remove communication
!
! input argument list:
! nn_obs - number of an observation data
! ndata - number of observations
! lunin - unit from which to read all obs of a given type
! lunout - unit to which to write subdomain specific observations
! mm1 - mpi task number + 1
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: periodic_s,nlon,nlat,jlon1,ilat1,istart,jstart
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: ndata,lunout,mm1,nobs
character(len=*) ,intent(in ) :: obsfile
character(len=*) ,intent(in ) :: obstypeall
! Declare local variables
integer(i_kind) lon,lat,lat_data,lon_data,n,k,lunin
integer(i_kind) jj,nreal,nchanl,nn_obs
integer(i_kind) ndata_s
integer(i_kind),dimension(mm1):: ibe,ibw,ibn,ibs
logical,allocatable,dimension(:):: luse_s
integer(i_kind),allocatable,dimension(:):: ioid_s
! IOID represents Initial-Obs-ID, which is the serial index of given
! observation record of a given observation source, before they are
! distributed according to the guess-grid partition. This ID is saved
! for all future needs of obs. redistribution, such that the sorted
! sequences can be reproduced.
real(r_kind),allocatable,dimension(:,:):: obs_data,data1_s
character(10):: obstype
character(20):: isis
! Read and write header
do k=1,mm1
! ibw,ibe,ibs,ibn west,east,south and north boundaries of total region
ibw(k)=jstart(k)-1
ibe(k)=jstart(k)+jlon1(k)-1
ibs(k)=istart(k)-1
ibn(k)=istart(k)+ilat1(k)-1
end do
lunin=11
open(lunin,file=trim(obsfile),form='unformatted')
read(lunin)obstype,isis,nreal,nchanl,lat_data,lon_data
if(trim(obstype) /=trim(obstypeall)) &
write(6,*)'DISOBS: ***ERROR*** obstype,obstypeall=',trim(obstype),trim(obstypeall)
nn_obs = nreal + nchanl
allocate(obs_data(nn_obs,ndata))
! Read in all observations of a given type along with subdomain flags
read(lunin) obs_data
close(lunin)
ndata_s=0
allocate(data1_s(nn_obs,nobs),luse_s(nobs),ioid_s(nobs))
! Loop over all observations. Locate each observation with respect
! to subdomains.
obs_loop: do n=1,ndata
lat=obs_data(lat_data,n)
lat=min(max(1,lat),nlat)
if(lat>=ibs(mm1).and.lat<=ibn(mm1))then
lon=obs_data(lon_data,n)
lon=min(max(0,lon),nlon)
if((lon >= ibw(mm1).and. lon <=ibe(mm1)) .or. &
(lon == 0 .and. ibe(mm1) >=nlon) .or. &
(lon == nlon .and. ibw(mm1) <=1) .or. periodic_s(mm1)) then
ndata_s=ndata_s+1
ioid_s(ndata_s)=n
luse_s(ndata_s)= .true.
do jj= 1,nn_obs
data1_s(jj,ndata_s) = obs_data(jj,n)
end do
prec_loop: do k=1,mm1-1
if(lat>=ibs(k).and.lat<=ibn(k)) then
if((lon >= ibw(k).and. lon <=ibe(k)) .or. &
(lon == 0 .and. ibe(k) >=nlon) .or. &
(lon == nlon .and. ibw(k) <=1) .or. periodic_s(k)) then
luse_s(ndata_s)= .false.
exit prec_loop
end if
end if
end do prec_loop
if(nobs == ndata_s) exit obs_loop
end if
end if
end do obs_loop
deallocate(obs_data)
! Write observations for given task to output file
write(lunout) obstypeall,isis,nreal,nchanl
write(lunout) data1_s,luse_s,ioid_s
deallocate(data1_s,luse_s,ioid_s)
return
end subroutine disobs
subroutine count_obs(ndata,nn_obs,lat_data,lon_data,obs_data,nobs_s)
!$$$ subprogram documentation block
! . . . .
! subprogram: count_obs counts number of observations on each subdomain
! prgmmr: derber date: 2014-12-16
!
! abstract: count observations in each pe subdomain
!
! program history log:
! 2014-12-16 derber
!
! input argument list:
! ndata - number of observations
! nn_obs - number of an observation data for each ob
! lat_data - location of lattitude
! lon_data - location of longitude
! obs_data - observation information array
!
! output argument list:
! nobs_s - number of observations on all subdomains
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: periodic_s,nlon,nlat,jlon1,ilat1,istart,jstart
use mpimod, only: npe
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: ndata,lat_data,lon_data
integer(i_kind) ,intent(in ) :: nn_obs
integer(i_kind),dimension(npe),intent(inout) :: nobs_s
real(r_kind),dimension(nn_obs,ndata),intent(in) :: obs_data
! Declare local variables
integer(i_kind) lon,lat,n,k
integer(i_kind),dimension(npe):: ibe,ibw,ibn,ibs
! Read and write header
do k=1,npe
! ibw,ibe,ibs,ibn west,east,south and north boundaries of total region
ibw(k)=jstart(k)-1
ibe(k)=jstart(k)+jlon1(k)-1
ibs(k)=istart(k)-1
ibn(k)=istart(k)+ilat1(k)-1
end do
nobs_s=0
! Loop over all observations. Locate each observation with respect
! to subdomains.
do n=1,ndata
lat=obs_data(lat_data,n)
lat=min(max(1,lat),nlat)
do k=1,npe
if(lat>=ibs(k).and.lat<=ibn(k)) then
lon=obs_data(lon_data,n)
lon=min(max(0,lon),nlon)
if((lon >= ibw(k).and. lon <=ibe(k)) .or. &
(lon == 0 .and. ibe(k) >=nlon) .or. &
(lon == nlon .and. ibw(k) <=1) .or. periodic_s(k)) then
nobs_s(k)=nobs_s(k)+1
end if
end if
end do
end do
return
end subroutine count_obs
! ------------------------------------------------------------------------
subroutine dislag(ndata,mm1,lunout,obsfile,obstypeall,ndata_s)
!$$$ subprogram documentation block
! . . . .
! subprogram: dislag distribute lagrangian observations into each pe
! subdomain
! prgmmr: lmeunier date: 2009-03-12
!
! abstract: distribute lagrangian observations into each pe subdomain
! (based on disobs). All observations of one balloon are
! are associated with a given processor
!
! program history log:
! 2009-03-12 lmeunier
!
! input argument list:
! ndata_s - number of observations in each pe sub-domain
! ndata - number of observations
! obsfile - unit from which to read all obs of a given type
! lunout - unit to which to write subdomain specific observations
! mm1 - mpi task number + 1
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use lag_fields, only: orig_lag_num
implicit none
! Declare passed variables
integer(i_kind),intent(in ) :: ndata,lunout,mm1
integer(i_kind),intent(inout) :: ndata_s
character(14) ,intent(in ) :: obsfile
character(10) ,intent(in ) :: obstypeall
! Declare local variables
integer(i_kind) num_data,n,lunin,num
integer(i_kind) jj,nreal,nchanl,nn_obs,ndatax
logical,allocatable,dimension(:):: luse,luse_s,luse_x
real(r_kind),allocatable,dimension(:,:):: obs_data,data1_s
character(10):: obstype
character(20):: isis
ndata_s=0
lunin=11
open(lunin,file=obsfile,form='unformatted')
read(lunin)obstype,isis,nreal,nchanl,num_data
if(obstype /=obstypeall) &
write(6,*)'DISLAG: ***ERROR*** obstype,obstypeall=',obstype,obstypeall
nn_obs = nreal + nchanl
allocate(obs_data(nn_obs,ndata))
! Read in all observations of a given type along with subdomain flags
read(lunin) obs_data
close(lunin)
allocate(luse(ndata),luse_x(ndata))
luse=.false.
luse_x=.false.
! Loop over all observations. Locate each observation with respect
! to subdomains.
use: do n=1,ndata
! Does the observation belong to the subdomain for this task?
num=int(obs_data(num_data,n),i_kind)
if ((mm1-1)==orig_lag_num(num,2)) then
ndata_s=ndata_s+1
luse(n)=.true.
luse_x(n)=.true. ! Never on another subdomain
end if
end do use
if(ndata_s > 0)then
allocate(data1_s(nn_obs,ndata_s),luse_s(ndata_s))
ndatax=0
do n=1,ndata
if(luse(n))then
ndatax=ndatax+1
luse_s(ndatax)=luse_x(n)
do jj= 1,nn_obs
data1_s(jj,ndatax) = obs_data(jj,n)
end do
end if
end do
! Write observations for given task to output file
write(lunout) obstypeall,isis,nreal,nchanl
write(lunout) data1_s,luse_s
deallocate(data1_s,luse_s)
endif
deallocate(obs_data,luse,luse_x)
return
end subroutine dislag