module read_obsmod
!$$$ module documentation block
! . . . .
! module: read_obsmod extra inquire routine for reading obs
! prgmmr: parrish org: np22 date: 2005-06-06
!
! abstract:
!
! program history log:
! 2009-01-05 todling - add gsi_inquire
! 2015-05-01 Liu Ling - Add call to read_rapidscat
! 2015-08-20 zhu - add flexibility for enabling all-sky and using aerosol info in radiance
! assimilation. Use radiance_obstype_search from radiance_mod.
! 2017-05-12 Y. Wang and X. Wang - add dbz to be read in, POC: xuguang.wang@ou.edu
!
! subroutines included:
! sub gsi_inquire - inquire statement supporting fortran earlier than 2003
! sub read_obs - read, select, reformat obs data
!
! Variable Definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
! set default to private
private
! set subroutines to public
public :: gsi_inquire
public :: read_obs
contains
subroutine gsi_inquire (lbytes,lexist,filename,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsi_inquire inquire file presence and size
! prgmmr: todling org: np22 date: 2009-01-05
!
! abstract: Inquire file presence and size; to be used when fortran
! 2003 not available or non-compliant.
!
! program history log:
! 2009-01-05 todling
! 2013-05-14 guo - changed the compiler #ifdef dependency on _size_
! inquire to a more portable way.
!
! input argument list:
! mype - mpi task id
! filename - input filename
!
! output argument list:
! lexist - file presence flag
! lbytes - file size (bytes)
!
! attributes:
! language: f90
! machine: Linux-cluster
!
!$$$ end documentation block
use kinds, only: i_kind,i_llong
implicit none
integer(i_llong),intent( out) :: lbytes
logical ,intent( out) :: lexist
character(len=*),intent(in ) :: filename
integer(i_kind) ,intent(in ) :: mype
character(len=256) command, fname
#if defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1110)
#define __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#endif
#ifdef __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#define _DO_NOT_SUPPORT_SIZE_INQUIRE_
#endif
lbytes=-1 ! in case that _size_ specifier is not supported.
#ifndef _DO_NOT_SUPPORT_SIZE_INQUIRE_
inquire(file=trim(filename),exist=lexist,size=lbytes)
! Note that the value of _size_ is defined by Fortran in "file storage units",
! which is not neccesary in byte units. It is not clear if this code had
! assumed the size value to be in byte units, or in whatever units.
#else
inquire(file=trim(filename),exist=lexist)
#endif
if(lexist) then
! Even with a compiler supporting 'size=' specifier, the size value may
! return -1, if a compiler considers that the size can not be determined.
! In that case, the size may be obtained through a user supported
! mechanism, such as reading the size from a system("wc -c") call result.
if(lbytes<0) then
write(fname,'(2a,i4.4)') 'fsize_',trim(filename),mype
write(command,'(4a)') 'wc -c ', trim(filename),' > ', trim(fname)
call system(command)
open(unit=999,file=trim(fname),form='formatted')
read(999,*) lbytes
close(999)
lexist = lbytes>0_i_llong ! skip this file if lbytes <=0
endif
endif
return
end subroutine gsi_inquire
subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_obs_check inquire file presence and size
! prgmmr: todling org: np22 date: 2010-03-05
!
! abstract: Reset file status depending on whether observation time
! matches analysis time and how offtime_date is set. This
! also checks for consistency in satellite data files and
! known types.
! WARNING: some of it looks inconsistent with long-window 4dvar
!
! program history log:
! 2009-xx-xx derber - originally placed inside inquire
! 2009-01-05 todling - move time/type-check out of inquire
! 2010-09-13 pagowski - add anow bufr and one obs chem
! 2013-01-26 parrish - WCOSS debug compile fails with satid not initialized.
! Set satid=1 at start of subroutine to allow debug compile.
! 2013-02-13 eliu - add ssmis
! 2013-07-01 todling/guo - allow user to bypass this check (old bufr support)
! 2014-10-01 ejones - add gmi and amsr2
! 2015-01-16 ejones - add saphir
! 2016-09-19 guo - properly initialized nread, in case of for quick-return cases.
! 2017-11-16 dutta - adding KOMPSAT5 bufr i.d for reading the data.
! 2019-03-27 h. liu - add abi
!
!
! input argument list:
! lexist - file status
! filename - input filename
! jsatid - satellite id
! dtype - satellite type
!
! output argument list:
! lexist - file status
!
! attributes:
! language: f90
! machine: Linux-cluster
!
!$$$ end documentation block
use kinds, only: i_kind,i_llong,r_kind,r_double
use gsi_4dvar, only: iadatebgn,iadateend
use obsmod, only: offtime_data
use convinfo, only: nconvtype,ictype,ioctype,icuse
use chemmod, only : oneobtest_chem,oneob_type_chem,&
code_pm25_ncbufr,code_pm25_anowbufr,code_pm10_ncbufr,code_pm10_anowbufr
implicit none
logical ,intent(inout) :: lexist
character(len=*),intent(in) :: filename
character(len=*),intent(in) :: jsatid
character(len=*),intent(in) :: dtype
integer(i_kind) ,intent(in) :: minuse
integer(i_kind) ,intent(inout) :: nread
integer(i_kind) :: lnbufr,idate,idate2,iret,kidsat
integer(i_kind) :: ireadsb,ireadmg,kx,nc,said
real(r_double) :: satid,rtype
character(8) subset
satid=1 ! debug executable wants default value ???
idate=0
nread=0; if(lexist) nread=1 ! in case of a quick return
#ifdef _SKIP_READ_OBS_CHECK_
return
#endif
if(trim(dtype) == 'tcp' .or. trim(filename) == 'tldplrso')return
if(trim(filename) == 'mitmdat' .or. trim(filename) == 'mxtmdat')return
if(trim(filename) == 'satmar')return
if(trim(dtype) == 'dbz' )return
! Use routine as usual
if(lexist .and. trim(dtype) /= 'tcp' )then
lnbufr = 15
open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
call readmg(lnbufr,subset,idate,iret)
if(iret == 0)then
! Extract date and check for consistency with analysis date
if (idate<iadatebgn.or.idate>iadateend) then
if(offtime_data) then
write(6,*)'***read_obs_check analysis and data file date differ, but use anyway'
else
write(6,*)'***read_obs_check*** ',&
'incompatable analysis and observation date/time',trim(filename),trim(dtype)
lexist=.false.
end if
write(6,*)'Analysis start :',iadatebgn
write(6,*)'Analysis end :',iadateend
write(6,*)'Observation time:',idate
endif
else
write(6,*)'***read_obs_check*** iret/=0 for reading date for ',trim(filename),dtype,jsatid,iret
lexist=.false.
end if
if(lexist)then
if(jsatid == '')then
kidsat=0
else if(jsatid == 'metop-a')then
kidsat=4
else if(jsatid == 'metop-b')then
kidsat=3
else if(jsatid == 'metop-c')then
kidsat=5
else if(jsatid == 'm08')then
kidsat = 55
else if(jsatid == 'm09')then
kidsat = 56
else if(jsatid == 'm10')then
kidsat = 57
else if(jsatid == 'm11')then
kidsat = 70
else if(jsatid == 'n08')then
kidsat=200
else if(jsatid == 'n09')then
kidsat=201
else if(jsatid == 'n10')then
kidsat=202
else if(jsatid == 'n11')then
kidsat=203
else if(jsatid == 'n12')then
kidsat=204
else if(jsatid == 'n14')then
kidsat=205
else if(jsatid == 'n15')then
kidsat=206
else if(jsatid == 'n16')then
kidsat=207
else if(jsatid == 'n17')then
kidsat=208
else if(jsatid == 'n18')then
kidsat=209
else if(jsatid == 'n19')then
kidsat=223
else if(jsatid == 'npp')then
kidsat=224
else if(jsatid == 'n20')then
kidsat=225
else if(jsatid == 'n21')then
kidsat=226
else if(jsatid == 'f08')then
kidsat=241
else if(jsatid == 'f10')then
kidsat=243
else if(jsatid == 'f11')then
kidsat=244
else if(jsatid == 'f13')then
kidsat=246
else if(jsatid == 'f14')then
kidsat=247
else if(jsatid == 'f15')then
kidsat=248
else if(jsatid == 'f16')then
kidsat=249
else if(jsatid == 'trmm')then
kidsat=282
else if(jsatid == 'f17')then
kidsat=285
else if(jsatid == 'f18')then
kidsat=286
else if(jsatid == 'f19')then
kidsat=287
else if(jsatid == 'g08' .or. jsatid == 'g08_prep')then
kidsat=252
else if(jsatid == 'g09' .or. jsatid == 'g09_prep')then
kidsat=253
else if(jsatid == 'g10' .or. jsatid == 'g10_prep')then
kidsat=254
else if(jsatid == 'g11' .or. jsatid == 'g11_prep')then
kidsat=255
else if(jsatid == 'g12' .or. jsatid == 'g12_prep')then
kidsat=256
else if(jsatid == 'g13' .or. jsatid == 'g13_prep')then
kidsat=257
else if(jsatid == 'g14' .or. jsatid == 'g14_prep')then
kidsat=258
else if(jsatid == 'g15' .or. jsatid == 'g15_prep')then
kidsat=259
else if(jsatid == 'g16' .or. jsatid == 'g16_prep')then
kidsat=270
else if(jsatid == 'g17' .or. jsatid == 'g17_prep')then
kidsat=271
else if(jsatid == 'n05')then
kidsat=705
else if(jsatid == 'n06')then
kidsat=706
else if(jsatid == 'n07')then
kidsat=707
else if(jsatid == 'tirosn')then
kidsat=708
else if ( jsatid == 'terra' ) then
kidsat = 783
else if ( jsatid == 'aqua' ) then
kidsat = 784
else if ( jsatid == 'aura' ) then
kidsat = 785
else if ( jsatid == 'gcom-w1' ) then
kidsat = 122
! Temporary comment gpm out here; discrepancy between SAID in bufr file and
! kidsat.
! else if ( jsatid == 'gpm' ) then
! kidsat = 288
else if ( jsatid == 'meghat' ) then
kidsat = 440
else
kidsat = 0
end if
call closbf(lnbufr)
open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
if(kidsat /= 0)then
lexist = .false.
satloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(ireadsb(lnbufr)==0)then
call ufbint(lnbufr,satid,1,1,iret,'SAID')
end if
if(nint(satid) == kidsat) then
lexist=.true.
exit satloop
end if
nread = nread + 1
end do satloop
else if(trim(filename) == 'prepbufr')then ! RTod: wired-in filename is not a good idea
lexist = .false.
fileloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloop
end if
end do
end do
nread = nread + 1
end do fileloop
else if(trim(filename) == 'wcpbufr')then
lexist = .false.
file2loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit file2loop
end if
end do
end do
nread = nread + 1
end do file2loop
else if(trim(filename) == 'gps_ref' .or. trim(filename) == 'gps_bnd')then
lexist = .false.
gpsloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(ireadsb(lnbufr)==0)then
call ufbint(lnbufr,satid,1,1,iret,'SAID')
end if
said=nint(satid)
if(((said > 739) .and.(said < 746)).or.(said == 820) .or. (said == 825) .or. &
(said == 786).or. (said == 4) .or.(said == 3).or. &
(said == 421).or. (said == 440).or.(said == 821)) then
lexist=.true.
exit gpsloop
end if
nread = nread + 1
end do gpsloop
else if(trim(filename) == 'prepbufr_profl')then
lexist = .false.
airploop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
if (trim(dtype)=='uv') then
if (kx==330 .or. kx==430 .or. kx==530) kx=230
if (kx==331 .or. kx==431 .or. kx==531) kx=231
if (kx==332 .or. kx==432 .or. kx==532) kx=232
if (kx==333 .or. kx==433 .or. kx==533) kx=233
if (kx==334 .or. kx==434 .or. kx==534) kx=234
if (kx==335 .or. kx==435 .or. kx==535) kx=235
else
if (kx==330 .or. kx==430 .or. kx==530) kx=130
if (kx==331 .or. kx==431 .or. kx==531) kx=131
if (kx==332 .or. kx==432 .or. kx==532) kx=132
if (kx==333 .or. kx==433 .or. kx==533) kx=133
if (kx==334 .or. kx==434 .or. kx==534) kx=134
if (kx==335 .or. kx==435 .or. kx==535) kx=135
end if
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit airploop
end if
end do
end do
nread = nread + 1
end do airploop
else if(trim(filename) == 'satwndbufr')then
lexist = .false.
loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034 and NC005039)
! are added as the GOES-R bufr file provide do not contain other winds.
! May not be necessary with the operational satwnd BUFR
if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or.&
trim(subset) == 'NC005070' .or. trim(subset) == 'NC005071' .or.&
trim(subset) == 'NC005044' .or. trim(subset) == 'NC005045' .or.&
trim(subset) == 'NC005046' .or. trim(subset) == 'NC005064' .or.&
trim(subset) == 'NC005065' .or. trim(subset) == 'NC005066' .or.&
trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or.&
trim(subset) == 'NC005032' .or. trim(subset) == 'NC005034' .or.&
trim(subset) == 'NC005039' .or. trim(subset) == 'NC005067' .or.&
trim(subset) == 'NC005068' .or. trim(subset) == 'NC005069' ) then
lexist = .true.
exit loop
endif
nread = nread + 1
end do loop
else if(trim(filename) == 'oscatbufr')then
lexist = .false.
oscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC012255') then
lexist = .true.
exit oscatloop
endif
end do oscatloop
else if(trim(filename) == 'rapidscatbufr')then
lexist = .false.
rapidscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC012255') then
lexist = .true.
exit rapidscatloop
endif
nread = nread + 1
end do rapidscatloop
else if(trim(filename) == 'hdobbufr')then
lexist = .false.
loop_hdob: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC004015') then
lexist = .true.
exit loop_hdob
endif
nread = nread + 1
end do loop_hdob
else if(trim(dtype) == 'pm2_5')then
if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
lexist=.true.
else
lexist = .false.
fileloopanow_pm2_5:do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
if (subset == 'ANOWPM') then
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
else if ( (subset == 'NC008031') .or. &
(subset == 'NC008032' ) ) then
call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
kx=nint(rtype)
if (kx/=code_pm25_ncbufr) then
cycle
else
kx=code_pm25_anowbufr
endif
else
cycle
endif
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. &
kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloopanow_pm2_5
end if
end do
end do
nread = nread + 1
enddo fileloopanow_pm2_5
endif
if (lexist) then
write(6,*)'found pm2_5 in anow bufr'
else
write(6,*)'did not find pm2_5 in anow bufr'
endif
else if(trim(dtype) == 'pm10')then
lexist = .false.
fileloopanow_pm10:do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
if (subset == 'NC008033') then
call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
kx=nint(rtype)
IF (kx/=code_pm10_ncbufr) then
cycle
else
kx=code_pm10_anowbufr
endif
else
cycle
endif
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. &
kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloopanow_pm10
end if
end do
end do
nread = nread + 1
enddo fileloopanow_pm10
if (lexist) then
write(6,*)'found pm10 in anow bufr'
else
write(6,*)'did not find pm10 in anow bufr'
endif
end if
end if
call closbf(lnbufr)
end if
if(lexist)then
write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype,jsatid
else
write(6,*)'read_obs_check: bufr file ',dtype,jsatid,' not available ',trim(filename)
end if
return
end subroutine read_obs_check
subroutine read_obs(ndata,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_obs read, select, reformat obs data
! prgmmr: parrish org: np22 date: 1990-10-07
!
! abstract: This routine is a driver for routines which read different
! types of observational data.
!
! program history log:
! 1990-10-07 parrish
! 1998-05-15 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-06-16 treadon - update documentation
! 2004-07-23 derber - modify to include conventional sst
! 2004-07-29 treadon - add only to module use, add intent in/out
! 2005-01-20 okamoto - add calling read_ssmi
! 2005-06-14 wu - add OMI oz
! 2005-07-06 derber - add mhs, hirs/4 and ears data
! 2005-08-02 derber - modify to use convinfo file
! 2005-09-08 derber - modify to use input group time window
! 2005-09-20 xu & pawlak - modify calling read_ssmis and read_amsre
! 2005-09-28 derber - modify to simplify obs handling
! 2005-10-17 derber - pass obs_load1 into read_amsre and read_ssmis
! 2005-10-18 treadon - remove obs_load and obs_load1
! 2005-10-20 kazumori - modify to read real AMSR-E data
! 2005-11-28 derber move determination of which ob data sets to read inside read_obs
! 2005-11-14 li, xu - modify sst obs read and add avhrr gac 1b obs read
! 2006-01-25 treadon - remove read_ieeetovs
! 2006-02-01 parrish - add getsfc and destroy_sfc for full surface fields
! 2006-02-03 derber - modify for new obs control and obs count
! 2005-02-03 treadon - gather guess 3d pressure to full grid array
! 2006-02-08 derber - modify to use new convinfo module
! 2006-02-01 liu - add ssu
! 2006-02-24 derber - modify to take advantage of convinfo module
! 2006-03-07 derber - consolidate processing of 1x1 and 5x5 goessndr
! 2006-04-20 kistler - moved conv_read to gsisub
! 2006-05-25 treadon - rename goesimg and goes_img and pcp_ssm/i as
! pcp_ssmi to make consistent with other obstype
! 2006-09-20 treadon - add mpi_io for select data file (obstype)s
! 2006-10-12 treadon - remove tendsflag check for pcp data (now in gsimain)
! 2007-02-21 sienkiewicz - bring in changes for MLS ozone
! 2007-03-15 su - add reading conventional error table option
! 2007-06-05 treadon - restructure mpi_querybf section to improve efficiency
! 2007-10-03 todling - skip most of this in 4dvar inner loop
! 2008-03-28 wu - move random seed for perturb_obs from setuprhsall
! 2008-04-18 safford - rm unused vars and uses
! 2008-05-01 h.liu - add gome ozone
! 2008-06-20 derber - move destroy_sfc to this routine
! 2008-09-08 lueken - merged ed''s cahnges into q1fy09 code
! 2008-12-30 todling - handle inquire for diff versions of fortran
! 2009-01-05 todling - need tendency alloc in observer mode
! 2009-01-23 todling - echo surface state info
! 2009-03-18 meunier - add a if statement to read lagrangian data
! 2009-08-19 guo - moved destroy_sfc_grid() to observer_finalize().
! 2009-12-20 gayno - modify argument lists so that fov-based surface
! calculation may be used.
! 2010-03-29 hu - add code to read in cloud observations including:
! prepbufr (metar, nesdis cloud product)
! radar reflectivity, lightning, NASA LaRC cloud
! 2010-04-01 treadon - move strip and reorder to gridmod
! 2010-04-08 hliu - add seviri
!
! 2010-04-05 huang - add aero and modis for reading modis aod from satellite
! currently read BUFR file only
! 2010-04-22 tangborn - read for carbon monoxide
! 2010-08-23 tong - add calcuation of hgtl_full used in read_radar.f90
! 2011-04-02 li - add nst_gsi, getnst and destroy_nst
! 2011-05-20 mccarty - add cris/atms handling
! 2011-05-26 todling - add call to create_nst
! 2012-03-07 akella - set NST variables; init & destroy in guess_grids
! [see create_sfc_grids & destroy_sfc_grids]
! 2012-05-16 todling - protect call to prebufr when NST is on (should be mods)
! 2012-09-08 wargan - add OMI with efficiency factors
! 2013-01-26 parrish - WCOSS debug compile fails--extra arguments in call read_aerosol.
! Commented out extra line of arguments not used.
! 2013-02-13 eliu - turn off parallel I/O for SSMIS (due to the need to
! do spatial averaging for noise reduction)
! 2013-06-01 zhu - add mype_airobst to handle aircraft temperature bias correction
! 2013-08-08 s.liu - add read NASA_LaRC_cloud product
! 2013-10-25 todling - reposition ltosi and others to commvars
! 2014-01-01 xli - add option to read NSST marine BUFR data file nsstbufr (on the
! top of prepbufr and modsbufr)
! 2014-02-03 guo - Hid processing (read) of non-EMC ozone obstypes
! through module m_extOzone, separated from read_ozone.
! - Added some -do- and -if- construct names, for easier
! understanding of the code.
! 2014-06-19 carley/zhu - Add tcamt and lcbas
! 2014-11-12 carley - Add call to read goes imager sky cover data for tcamt
! 2014-12-03 derber - modify for no radiance cases and read processor for
! surface fields
! 2015-01-16 ejones - added saphir, gmi, and amsr2 handling
! 2015-03-23 zaizhong ma - add Himawari-8 ahi
! 2015-05-30 li - modify for no radiance cases but sst (nsstbufr) and read processor for
! surface fields (use_sfc = .true. for data type of sst),
! to use deter_sfc in read_nsstbufr.f90)
! 2015-07-10 pondeca - add cloud ceiling height (cldch)
! 2015-08-12 pondeca - add capability to read min/maxT obs from ascii file
! 2015-08-20 zhu - use centralized radiance info from radiance_mod to specify
! all-sky and aerosol usages in radiance assimilation
! 2015-09-04 J. Jung - Added mods for CrIS full spectral resolution (FSR)
! 2015-10-19 lippi - Added logic to read l2rw or rw radial winds.
! 2016-03-02 s.liu/carley - remove use_reflectivity and use i_gsdcldanal_type
! 2016-04-28 J. Jung - added logic for RARS and direct broadcast data from NESDIS/UW.
! 2016-05-05 pondeca - add 10-m u-wind and v-wind (uwnd10m, vwnd10m)
! 2016-09-19 Guo - replaced open(obs_input_common) with "call unformatted_open(obs_input_common)"
! 2017-05-12 Y. Wang and X. Wang - add multi-interface to read in dBZ (nc) and radial velocity (ascii)
! 2016-12-14 lippi - Fixed bug of using observations twice when both
! l2rwbufr and radarbufr are in the OBS_INPUT table.
! Changed the dsis entries for l2rwbufr and radarbufr to
! l2rw and l3rw respectively. Also make use of nml
! option vadwnd_l2rw_qc.
! 2017-08-31 Li - move gsi_nstcoupler_init & gsi_nstcoupler_read to getsfc in sathin.F90
! - move gsi_nstcoupler_final from create_sfc_grids to here
! 2018-01-23 Apodaca - add GOES/GLM lightning data
! 2019-01-15 Li - add to handle mbuoyb
! 2019-03-27 h. liu - add abi
!
!
! input argument list:
! mype - mpi task id
!
! output 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
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,i_llong,r_double
use gridmod, only: nlon,nlat,nsig,iglobal,ijn,itotsub,lat1,lon1,&
displs_g,strip,reorder
use general_commvars_mod, only: ltosi,ltosj
use obsmod, only: iadate,ndat,time_window,dplat,dsfcalc,dfile,dthin, &
dtype,dval,dmesh,obsfile_all,ref_obs,nprof_gps,dsis,ditype,&
perturb_obs,lobserver,lread_obs_save,obs_input_common, &
reduce_diag,nobs_sub,dval_use
use gsi_nstcouplermod, only: nst_gsi,gsi_nstcoupler_final
use qcmod, only: njqc,vadwnd_l2rw_qc
use gsi_4dvar, only: l4dvar
use satthin, only: super_val,super_val1,superp,makegvals,getsfc,destroy_sfc
use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_rtype,mpi_integer,npe,&
setcomm
use constants, only: one,zero
use converr, only: converr_read
use converr_ps, only: converr_ps_read
use converr_q, only: converr_q_read
use converr_t, only: converr_t_read
use converr_uv, only: converr_uv_read
use converr_pw, only: converr_pw_read
use convb_ps,only: convb_ps_read
use convb_q,only:convb_q_read
use convb_t,only:convb_t_read
use convb_uv,only:convb_uv_read
use guess_grids, only: ges_prsl,geop_hgtl,ntguessig
use radinfo, only: nusis,iuse_rad,jpch_rad,diag_rad
use insitu_info, only: mbuoy_info,mbuoyb_info,read_ship_info
use aeroinfo, only: nusis_aero,iuse_aero,jpch_aero,diag_aero
use ozinfo, only: nusis_oz,iuse_oz,jpch_oz,diag_ozone
use pcpinfo, only: npcptype,nupcp,iusep,diag_pcp
use convinfo, only: nconvtype,ioctype,icuse,diag_conv,ithin_conv
use chemmod, only : oneobtest_chem,oneob_type_chem,oneobschem
use lightinfo, only: nlighttype,iuse_light,diag_light
use aircraftinfo, only: aircraft_t_bc,aircraft_t_bc_pof,aircraft_t_bc_ext,mype_airobst
use gsi_io, only: mype_io
use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type
use radiance_mod, only: rad_obs_type,radiance_obstype_search
use m_extOzone, only: is_extOzone
use m_extOzone, only: extOzone_read
use mpeu_util, only: warn
use gsi_unformatted, only: unformatted_open
use mrmsmod,only: l_mrms_sparse_netcdf
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: mype
integer(i_kind),dimension(ndat,3),intent( out) :: ndata
! Declare local parameters
integer(i_llong),parameter:: lenbuf=8388608_i_llong ! lenbuf=8*1024*1024
! Declare local variables
logical :: lexist,ssmis,amsre,sndr,hirs,avhrr,lexistears,lexistdb,use_prsl_full,use_hgtl_full
logical :: use_sfc,nuse,use_prsl_full_proc,use_hgtl_full_proc,seviri,mls,abi
logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible
logical :: modis,use_sfc_any
logical :: acft_profl_file
character(10):: obstype,platid
character(22):: string
character(120):: infile
character(20):: sis
integer(i_kind) i,j,k,ii,nmind,lunout,isfcalc,ithinx,ithin,nread,npuse,nouse
integer(i_kind) nprof_gps1,npem1,krsize,len4file,npemax,ilarge,nlarge,npestart
integer(i_llong) :: lenbytes
integer(i_kind):: npetot,npeextra,mmdat,nodata
integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix
integer(i_kind):: mype_root
integer(i_kind):: minuse,lunsave,maxproc,minproc
integer(i_kind),dimension(ndat):: npe_sub,npe_sub3,mpi_comm_sub,mype_root_sub,npe_order
integer(i_kind),dimension(ndat):: ntasks1,ntasks
integer(i_kind),dimension(ndat):: read_rec1,read_rec
integer(i_kind),dimension(ndat):: read_ears_rec1,read_ears_rec
integer(i_kind),dimension(ndat):: read_db_rec1,read_db_rec
integer(i_kind),dimension(ndat,3):: ndata1
integer(i_kind),dimension(npe,ndat):: mype_work,nobs_sub1
integer(i_kind),dimension(npe,ndat):: mype_sub
integer(i_kind),allocatable,dimension(:):: nrnd
integer(i_kind):: nmls_type,mype_io_sfc
integer(i_kind):: iread,ipuse,iouse
real(r_kind) gstime,val_dat,rmesh,twind,rseed
real(r_kind),allocatable,dimension(:) :: prslsm,hgtlsm,work1
real(r_kind),allocatable,dimension(:,:,:):: prsl_full,hgtl_full
type(rad_obs_type) :: radmod
data lunout / 81 /
data lunsave / 82 /
!*****************************************************************************
! Set analysis time and allocate/initialize arrays and variables
call w3fs21(iadate,nmind)
gstime=real(nmind,r_kind)
allocate(nobs_sub(npe,ndat))
nobs_sub = 0
nobs_sub1 = 0
call makegvals
do ii=1,ndat
ndata1(ii,1)=0
ndata1(ii,2)=0
ndata1(ii,3)=0
ntasks1(ii) =0
parallel_read=.false.
end do
npem1=npe-1
nprof_gps1=0
if(njqc) then
call converr_ps_read(mype)
call converr_q_read(mype)
call converr_t_read(mype)
call converr_uv_read(mype)
call converr_pw_read(mype)
call convb_ps_read(mype)
call convb_q_read(mype)
call convb_t_read(mype)
call convb_uv_read(mype)
else
call converr_read(mype)
endif
! Optionally set random seed to perturb observations
if (perturb_obs) then
rseed=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000+mype
call random_seed(size=krsize)
allocate(nrnd(krsize))
do i=1,krsize
nrnd(i)=rseed
end do
call random_seed(put=nrnd)
deallocate(nrnd)
endif
! Set data class and number of reader tasks. Set logical flag to indicate
! type type of GPS data (if present)
ii=0
ref_obs = .false. !.false. = assimilate GPS bending angle
ears_possible = .false.
db_possible = .false.
nmls_type=0
read_rec1 = 0
read_ears_rec1=0
read_db_rec1=0
do i=1,ndat
obstype=dtype(i) ! obstype - observation types to process
amsre= index(obstype,'amsre') /= 0
ssmis= index(obstype,'ssmis') /= 0
sndr = index(obstype,'sndr') /= 0
hirs = index(obstype,'hirs') /= 0
avhrr = index(obstype,'avhrr') /= 0
modis = index(obstype,'modis') /= 0
seviri = index(obstype,'seviri') /= 0
abi = index(obstype,'abi') /= 0
mls = index(obstype,'mls') /= 0
if(obstype == 'mls20' ) nmls_type=nmls_type+1
if(obstype == 'mls22' ) nmls_type=nmls_type+1
if(obstype == 'mls30' ) nmls_type=nmls_type+1
if(nmls_type>1) then
write(6,*) '******ERROR***********: there is more than one MLS data type, not allowed, please check'
call stop2(339)
end if
if (obstype == 't' .or. obstype == 'uv' .or. &
obstype == 'q' .or. obstype == 'ps' .or. &
obstype == 'pw' .or. obstype == 'spd'.or. &
obstype == 'sst'.or. &
obstype == 'tcp'.or. obstype == "lag".or. &
obstype == 'dw' .or. obstype == 'rw' .or. &
obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
obstype == 'rad_ref' .or. obstype=='lghtn' .or. &
obstype == 'larccld' .or. obstype == 'pm2_5' .or. obstype == 'pm10' .or. &
obstype == 'gust' .or. obstype=='vis' .or. &
obstype == 'pblh' .or. obstype=='wspd10m' .or. &
obstype == 'td2m' .or. obstype=='mxtm' .or. &
obstype == 'mitm' .or. obstype=='pmsl' .or. &
obstype == 'howv' .or. obstype=='tcamt' .or. &
obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. &
obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then
ditype(i) = 'conv'
else if (obstype == 'swcp' .or. obstype == 'lwcp') then
ditype(i) = 'wcp'
else if( hirs .or. sndr .or. seviri .or. abi .or. &
obstype == 'airs' .or. obstype == 'amsua' .or. &
obstype == 'msu' .or. obstype == 'iasi' .or. &
obstype == 'amsub' .or. obstype == 'mhs' .or. &
obstype == 'hsb' .or. obstype == 'goes_img' .or. &
obstype == 'ahi' .or. avhrr .or. &
amsre .or. ssmis .or. obstype == 'ssmi' .or. &
obstype == 'ssu' .or. obstype == 'atms' .or. &
obstype == 'cris' .or. obstype == 'cris-fsr' .or. &
obstype == 'amsr2' .or. &
obstype == 'gmi' .or. obstype == 'saphir' ) then
ditype(i) = 'rad'
else if (is_extOzone(dfile(i),obstype,dplat(i))) then
ditype(i) = 'ozone'
else if (obstype == 'sbuv2' &
.or. obstype == 'omi' &
.or. obstype == 'ompstc8' &
.or. obstype == 'ompsnp' &
.or. obstype == 'gome' &
.or. mls &
) then
ditype(i) = 'ozone'
else if (obstype == 'mopitt') then
ditype(i) = 'co'
else if (index(obstype,'pcp')/=0 )then
ditype(i) = 'pcp'
else if (obstype == 'gps_ref' .or. obstype == 'gps_bnd') then
ditype(i) = 'gps'
else if ( index(obstype,'aod') /= 0 ) then
ditype(i) = 'aero'
else if (obstype == 'goes_glm') then
ditype(i) = 'light'
else
write(6,*)'READ_OBS: ***ERROR*** - unknown ob type ',trim(obstype)
end if
! Set data class and number of reader tasks. Set logical flag to indicate
! type type of GPS data (if present)
if (index(dtype(i),'gps_ref') /= 0) ref_obs = .true.
! Check info files to see if data is read.
nuse=.false.
minuse=-1
if(ditype(i) == 'conv')then
if(diag_conv .and. .not. reduce_diag)minuse=-2
do j=1,nconvtype
if(trim(dtype(i)) == trim(ioctype(j)) .and. icuse(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'rad')then
if(diag_rad .and. .not. reduce_diag)minuse=-2
do j=1,jpch_rad
if(trim(dsis(i)) == trim(nusis(j)) .and. iuse_rad(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'ozone')then
if(diag_ozone .and. .not. reduce_diag)minuse=-2
do j=1,jpch_oz
if(trim(dsis(i)) == trim(nusis_oz(j)) .and. iuse_oz(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'pcp')then
if(diag_pcp .and. .not. reduce_diag)minuse=-2
do j=1,npcptype
if(trim(dsis(i)) == trim(nupcp(j)) .and. iusep(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'aero')then
if(diag_aero .and. .not. reduce_diag)minuse=-2
do j=1,jpch_aero
if(trim(dsis(i)) == trim(nusis_aero(j)) .and. iuse_aero(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'light')then
if(diag_light .and. .not. reduce_diag)minuse=-2
do j=1,nlighttype
if(iuse_light(j) > minuse) nuse=.true.
end do
else
nuse=.true.
end if
if(nuse)then
! Control parallel read for each ob type (currently just rad obs).
! To remove parallel read comment out line.
ithin=dthin(i)
if(ithin > 0 )then
if(dmesh(ithin) > one)then
if(hirs)then
parallel_read(i)= .true.
else if(obstype == 'amsua')then
parallel_read(i)= .true.
else if(obstype == 'airs' )then
parallel_read(i)= .true.
else if(obstype == 'iasi')then
parallel_read(i)= .true.
else if(obstype == 'amsub')then
parallel_read(i)= .true.
else if(obstype == 'mhs' )then
parallel_read(i)= .true.
else if(sndr )then
parallel_read(i)= .true.
! N.B. ATMS must be run on one processor for the filtering code to work.
else if(obstype == 'atms')then
! parallel_read(i)= .true.
else if(ssmis)then
! parallel_read(i)= .true.
else if(seviri)then
parallel_read(i)= .true.
else if(abi)then
parallel_read(i)= .true.
else if(obstype == 'cris' .or. obstype == 'cris-fsr')then
parallel_read(i)= .true.
else if(avhrr)then
parallel_read(i)= .true.
else if(amsre)then
parallel_read(i)= .true.
else if(obstype == 'goes_img' )then
parallel_read(i)= .true.
else if(obstype == 'ahi' )then
parallel_read(i)= .true.
else if(obstype == 'hsb' )then
parallel_read(i)= .true.
else if(obstype == 'ssmi' )then
parallel_read(i)= .true.
else if(obstype == 'ssu' )then
parallel_read(i)= .true.
else if(obstype == 'amsr2')then
! parallel_read(i)= .true. ! turn parallel read off for spatial averaging
else if(obstype == 'gmi')then
! parallel_read(i)= .true. ! turn parallel read off for spatial averaging
! Parallel read for SAPHIR not currently working. Leave parallel read off.
! else if(obstype == 'saphir')then
! parallel_read(i)= .true.
end if
end if
end if
! direct broadcast from EUMETSAT (EARS)
ears_possible(i) = ditype(i) == 'rad' .and. &
(obstype == 'amsua' .or. obstype == 'amsub' .or. &
obstype == 'mhs' .or. obstype == 'hirs3' .or. &
obstype == 'cris' .or. obstype == 'cris-fsr' .or. &
obstype == 'iasi' .or. obstype == 'atms') .and. &
(dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. &
dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. &
dplat(i) == 'n20' .or. &
dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. &
dplat(i) == 'metop-c')
! direct broadcast from NESDIS/UW
db_possible(i) = ditype(i) == 'rad' .and. &
(obstype == 'amsua' .or. obstype == 'amsub' .or. &
obstype == 'mhs' .or. obstype == 'atms' .or. &
obstype == 'cris' .or. obstype == 'cris-fsr' .or. &
obstype == 'iasi') .and. &
(dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. &
dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. &
dplat(i) == 'n20' .or. &
dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. &
dplat(i) == 'metop-c')
! Inquire data set to deterimine if input data available and size of dataset
ii=ii+1
if (ii>npem1) ii=0
if(mype==ii)then
call gsi_inquire(lenbytes,lexist,trim(dfile(i)),mype)
call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i))
! If no data set starting record to be 999999. Note if this is not large
! enough code should still work - just does a bit more work.
if(.not. lexist)read_rec1(i) = 999999
len4file=lenbytes/4
if (ears_possible(i))then
call gsi_inquire(lenbytes,lexistears,trim(dfile(i))//'ears',mype)
call read_obs_check (lexistears,trim(dfile(i))//'ears',dplat(i),dtype(i),minuse, &
read_ears_rec1(i))
! If no data set starting record to be 999999. Note if this is not large
! enough code should still work - just does a bit more work.
if(.not. lexistears)read_ears_rec1(i) = 999999
lexist=lexist .or. lexistears
len4file=len4file+lenbytes/4
else
read_ears_rec1(i) = 999999
end if
if (db_possible(i))then
call gsi_inquire(lenbytes,lexistdb,trim(dfile(i))//'_db',mype)
call read_obs_check (lexistdb,trim(dfile(i))//'_db',dplat(i),dtype(i),minuse, &
read_db_rec1(i))
! If no data set starting record to be 999999. Note if this is not large
! enough code should still work - just does a bit more work.
if(.not. lexistdb)read_db_rec1(i) = 999999
lexist=lexist .or. lexistdb
len4file=len4file+lenbytes/4
else
read_db_rec1(i) = 999999
end if
if(lexist) then
! Initialize number of reader tasks to 1. For the time being
! only allow number of reader tasks >= 1 for select obstype.
ntasks1(i)=1
if(parallel_read(i)) then
! Allow up to 16 processors/file increase loop bounds to increase number of processors allowed
do j=1,4
if(len4file < lenbuf)exit
ntasks1(i)=2*ntasks1(i)
len4file=len4file/2
end do
! if(ntasks1(i)*lenbuf < len4file) ntasks1(i)=ntasks1(i)+1
end if
end if
end if
else
if(mype == 0)write(6,*) 'data type ',dsis(i), &
'not used in info file -- do not read file ',trim(dfile(i))
end if
end do
! Distribute optimal number of reader tasks to all mpi tasks
call mpi_allreduce(ntasks1,ntasks,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(read_rec1,read_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(read_ears_rec1,read_ears_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(read_db_rec1,read_db_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror)
! Limit number of requested tasks per type to be <= total available tasks
npemax=0
npetot=0
do i=1,ndat
if (ntasks(i)>npe) then
write(6,*)'read_obs: ***WARNING*** i=',i,' dtype=',dtype(i),' dsis=',dsis(i),&
' requested ntasks=',ntasks(i),' > npe=',npe,' reset ntasks=',npe
ntasks(i)=npe
endif
npe_sub(i)=ntasks(i)
npetot=npetot+npe_sub(i)
npemax=max(npemax,npe_sub(i))
end do
if(l4dvar.and.(.not.lobserver)) then
!_RTod use_sfc=.false.
!_RTod call getsfc(mype,use_sfc)
return
endif
npeextra=0
if(mod(npetot,npe) > 0) npeextra=npe-mod(npetot,npe)
maxproc=32
if(npeextra > 0)then
if(mype == 0)write(6,*) ' number of extra processors ',npeextra
extraloop: do j=1,npeextra
npe_sub3=ntasks
minproc=999999
do i=1,ndat
if(ntasks(i) > 0 .and. parallel_read(i))minproc=min(minproc,ntasks(i))
end do
if(npeextra < minproc) exit extraloop
iix=minproc
do ii=1,6
do i=1,ndat
if(iix == npe_sub3(i) .and. parallel_read(i))then
if(ntasks(i) <= npeextra .and. ntasks(i) < maxproc)then
npeextra=npeextra-ntasks(i)
ntasks(i)=2*ntasks(i)
if(npeextra < iix)cycle extraloop
end if
end if
end do
iix=2*iix
if(iix >= maxproc) cycle extraloop
end do
end do extraloop
end if
! Set up locations of first processor
ilarge=0
npestart=0
npe_sub3=npe_sub
mype_root_sub=0
mmdat=0
loopx: do j=1,ndat
nlarge=0
do i=1,ndat
if(npe_sub3(i) > nlarge .and. npe_sub3(i)+npestart <= npe)then
ilarge=i
nlarge=npe_sub3(i)
end if
end do
if(nlarge == 0)exit loopx
npe_order(j)=ilarge
mype_root_sub(ilarge)=npestart
npestart=npestart+npe_sub3(ilarge)
mmdat=mmdat+1
npe_sub3(ilarge)=0
if(npestart + minval(npe_sub3, mask=npe_sub3>0)>= npe) npestart=0
end do loopx
! Define sub-communicators for each data file
mm1=mype+1
belong=.false.
mype_sub=-999
mype_root=0
next_mype=0
do ii=1,mmdat
i=npe_order(ii)
if(npe_sub(i) > 0)then
next_mype=mype_root_sub(i)
do k=1,npe_sub(i)
mype_work(k,i) = next_mype
mype_sub(mype_work(k,i)+1,i)=k-1
if(next_mype == mype)belong(i) = .true.
next_mype = next_mype + 1
if (next_mype>npem1) next_mype=0
end do
call setcomm(iworld,iworld_group,npe_sub(i),mype_work(1,i),&
mpi_comm_sub(i),ierror)
end if
end do
mype_airobst = mype_root
do ii=1,mmdat
i=npe_order(ii)
if(mype == 0 .and. npe_sub(i) > 0) write(6,'(1x,a,i4,1x,a,1x,2a,2i4,1x,i6,1x,i6,1x,i6)') &
'READ_OBS: read ',i,dtype(i),dsis(i),' using ntasks=',ntasks(i),mype_root_sub(i), &
read_rec(i),read_ears_rec(i),read_db_rec(i)
acft_profl_file = index(dfile(i),'_profl')/=0
if ((aircraft_t_bc_pof .or. aircraft_t_bc_ext .or. &
(aircraft_t_bc .and. acft_profl_file)) .and. dtype(i) == 't') &
mype_airobst = mype_root_sub(i)
end do
use_prsl_full=.false.
use_hgtl_full=.false.
use_sfc=.false.
use_prsl_full_proc=.false.
use_hgtl_full_proc=.false.
mype_io_sfc=mype_io
do ii=1,mmdat
i=npe_order(ii)
if(ditype(i) =='conv')then
obstype=dtype(i)
if (obstype == 't' .or. obstype == 'q' .or. &
obstype == 'uv' .or. obstype == 'wspd10m' .or. &
obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then
use_prsl_full=.true.
if(belong(i))use_prsl_full_proc=.true.
else
do j=1,nconvtype
if(obstype == trim(ioctype(j)) .and. ithin_conv(j) > 0)then
use_prsl_full=.true.
if(belong(i))use_prsl_full_proc=.true.
end if
end do
end if
if(obstype == 'rw')then
use_hgtl_full=.true.
if(belong(i))use_hgtl_full_proc=.true.
else if(obstype == 'dbz')then
use_hgtl_full=.true.
if(belong(i))use_hgtl_full_proc=.true.
end if
if(obstype == 'sst')then
if(belong(i))use_sfc=.true.
endif
else if(ditype(i) == 'rad' )then
if(belong(i)) use_sfc=.true.
end if
end do
use_sfc_any=.false.
loop: do ii=1,mmdat
i=npe_order(ii)
if(ditype(i) == 'rad' .or. ditype(i) == 'sst')then
mype_io_sfc=mype_root_sub(i)
use_sfc_any=.true.
exit loop
end if
end do loop
if(use_sfc_any .and. mype == mype_io)use_sfc=.true.
! Create full horizontal surface fields from local fields in guess_grids
call getsfc(mype,mype_io_sfc,use_sfc,use_sfc_any)
if(mype == mype_io) call prt_guessfc2('sfcges2',use_sfc)
! Get guess 3d pressure on full grid
allocate(work1(max(iglobal,itotsub)),prslsm(lat1*lon1))
if(use_prsl_full)then
if(use_prsl_full_proc)then
allocate(prsl_full(nlat,nlon,nsig))
else
allocate(prsl_full(1,1,1))
end if
do k=1,nsig
call strip(ges_prsl(:,:,k,ntguessig),prslsm)
call mpi_allgatherv(prslsm,ijn(mype+1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,mpi_comm_world,ierror)
if(use_prsl_full_proc)then
call reorder(work1,1,1)
do ii=1,iglobal
i=ltosi(ii)
j=ltosj(ii)
prsl_full(i,j,k)=work1(ii)
end do
end if
end do
else
allocate(prsl_full(1,1,1))
end if
! Get guess 3d geopotential height on full grid
if(use_hgtl_full)then
allocate(hgtlsm(lat1*lon1))
if(use_hgtl_full_proc)then
allocate(hgtl_full(nlat,nlon,nsig))
else
allocate(hgtl_full(1,1,1))
end if
do k=1,nsig
call strip(geop_hgtl(:,:,k,ntguessig),hgtlsm)
call mpi_allgatherv(hgtlsm,ijn(mype+1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,mpi_comm_world,ierror)
if(use_hgtl_full_proc)then
call reorder(work1,1,1)
do ii=1,iglobal
i=ltosi(ii)
j=ltosj(ii)
hgtl_full(i,j,k)=work1(ii)
end do
end if
end do
deallocate(hgtlsm)
else
allocate(hgtl_full(1,1,1))
end if
deallocate(work1,prslsm)
! Create moored buoy station ID
call mbuoy_info(mype)
! Create moored buoy station ID for mbuoyb with 7-digit station ID
call mbuoyb_info(mype)
! Create ships info(ID, Depth & Instrument)
call read_ship_info(mype)
! Loop over data files. Each data file is read by a sub-communicator
loop_of_obsdata_files: &
do ii=1,mmdat
i=npe_order(ii)
task_belongs: &
if (i > 0 .and. belong(i)) then
platid=dplat(i) ! platid - satellites to read
obstype=dtype(i) ! obstype - observation types to process
infile=trim(dfile(i)) ! infile - units from which to read data
sis=dsis(i) ! sensor/instrument/satellite indicator
val_dat=dval(i) ! weighting factors applied to super obs
ithin=dthin(i) ! ithin - flags to thin data
ithinx=max(1,abs(ithin))
rmesh=dmesh(ithinx) ! rmesh - thinning mesh sizes (km)
twind=time_window(i) ! time window (hours) for input group
isfcalc=dsfcalc(i) ! method to calculate surface fields within fov
nread=0
nouse=0
npuse=0
if (mype_sub(mm1,i)==mype_root) then
open(lunout,file=obsfile_all(i),form='unformatted')
rewind(lunout)
endif
! Process conventional (prepbufr) data
string="--"//trim(ditype(i))
ditype_select: &
if(ditype(i) == 'conv')then
conv_obstype_select: &
if (obstype == 't' .or. obstype == 'q' .or. obstype == 'ps' .or. &
obstype == 'pw' .or. obstype == 'spd'.or. &
obstype == 'gust' .or. obstype == 'vis'.or. &
obstype == 'td2m' .or. &
! obstype=='mxtm' .or. obstype == 'mitm' .or. &
! obstype=='howv' .or. obstype=='pmsl' .or. & #howv #ww3
obstype=='pmsl' .or. &
obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
obstype == 'lcbas' .or. obstype == 'cldch' ) then
! Process flight-letel high-density data not included in prepbufr
if ( index(infile,'hdobbufr') /=0 ) then
call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_FL_HDOB'
else
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
else if(obstype == 'howv') then
if ( index(infile,'satmar') /=0) then
call read_satmar (nread, npuse, nouse &
, infile, obstype, lunout, gstime, twind, sis &
, nobs_sub1(1,i))
string='READ_SATMAR'
else
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
else if(obstype == 'mitm') then
if ( index(infile,'mitmdat') /=0) then
call read_mitm_mxtm(nread,npuse,nouse,infile,obstype,lunout,gstime,sis, &
nobs_sub1(1,i))
string='READ_ASCII_MITM'
else
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
else if(obstype == 'mxtm') then
if ( index(infile,'mxtmdat') /=0) then
call read_mitm_mxtm(nread,npuse,nouse,infile,obstype,lunout,gstime,sis, &
nobs_sub1(1,i))
string='READ_ASCII_MXTM'
else
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
! Process total cloud amount (tcamt) in prepbufr -or- from goes imager sky cover products
else if(obstype == 'tcamt') then
! Process GOES Imager Sky Cover product separately from prepbufr-based sky cover obs
if ( index(infile,'goessky') /=0 ) then
call read_goesimgr_skycover(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_GOESIMGR_SKYCOVER'
else
! else read from prepbufr
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,prsl_full, &
nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
end if
! Process winds in the prepbufr
else if(obstype == 'uv' .or. obstype == 'wspd10m' .or. &
obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then
! Process satellite winds which seperate from prepbufr
if ( index(infile,'satwnd') /=0 ) then
call read_satwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_SATWND'
! Process oscat winds which seperate from prepbufr
elseif ( index(infile,'oscatbufr') /=0 ) then
call read_sfcwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_SFCWND'
! Process rapidscat winds which seperate from prepbufr
elseif ( index(infile,'rapidscatbufr') /=0 ) then
call read_rapidscat(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_RAPIDSCAT'
else if ( index(infile,'hdobbufr') /=0 ) then
call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs_sub1(1,i))
string='READ_FL_HDOB'
else
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
! Process conventional SST (nsstbufr, at this moment) data
elseif ( obstype == 'sst' ) then
string="--"//trim(ditype(i))//":sst:"//trim(platid)
if ( platid == 'nsst') then
call read_nsstbufr(nread,npuse,nouse,gstime,infile,obstype, &
lunout,twind,sis,nobs_sub1(1,i))
string='READ_NSSTBUFR'
elseif ( platid == 'mods') then
call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, &
lunout,twind,sis,nobs_sub1(1,i))
string='READ_MODSBUFR'
else
if(nst_gsi>0)then
write(6,*)'read_obs: should not handle SST via read_prepbufr when NSST on'
call stop2(999)
endif
call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_PREPBUFR'
endif
! Process radar reflectivity Mosaic
else if (obstype == 'rad_ref' ) then
call read_radarref_mosaic(nread,npuse,infile,obstype,lunout,twind,sis, &
nobs_sub1(1,i))
string='READ_RADARREF_MOSAIC'
! Process lightning
else if (obstype == 'lghtn' ) then
if(i_gsdcldanal_type==2) then
call read_lightning(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
else if(i_gsdcldanal_type==1 .or. i_gsdcldanal_type==6 &
.or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7) then
call read_lightning_grid(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
endif
string='READ_LIGHTNING'
! Process NASA LaRC
! for regional obs that are already mapped to analysis grid
else if (obstype == 'larccld' ) then
if(i_gsdcldanal_type==2) then
call read_NASA_LaRC_cloud(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
else if(i_gsdcldanal_type==1 .or. i_gsdcldanal_type==6 &
.or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7 .or. i_gsdcldanal_type==99) then
call read_nasa_larc(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
end if
string='READ_NASA_LaRC'
! for global NASA LaRC obs
else if (obstype == 'larcglb' ) then
call read_NASA_LaRC_cloud(nread,npuse,nouse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
string='READ_NASA_LaRC'
! Process radar winds
else if (obstype == 'rw') then
if( trim(infile) == 'vr_vol' )then
call read_radar_wind_ascii(nread,npuse,nouse,infile,lunout,obstype,sis,&
hgtl_full,nobs_sub1(1,i))
else
if (vadwnd_l2rw_qc) then
write(6,*)'READ_OBS: radial wind,read_radar,dfile=',infile,',dsis=',sis
call read_radar(nread,npuse,nouse,infile,lunout,obstype,twind,sis,&
hgtl_full,nobs_sub1(1,i))
string='READ_RADAR'
else if (sis == 'l2rw') then
write(6,*)'READ_OBS: radial wind,read_radar_l2rw_novadqc,dfile=',infile,',dsis=',sis
call read_radar_l2rw_novadqc(npuse,nouse,lunout,obstype,sis,nobs_sub1(1,i))
string='READ_RADAR_L2RW_NOVADQC'
end if
end if
string='READ_RADAR_WIND'
! Process radar reflectivity from MRMS
else if (obstype == 'dbz' ) then
print *, "calling read_dbz"
if(trim(infile)=='dbzobs.nc')then
call read_dbz_nc(nread,npuse,nouse,infile,lunout,obstype,sis,hgtl_full,nobs_sub1(1,i))
string='READ_dBZ'
else
call read_dbz_mrms_detect_format(infile,l_mrms_sparse_netcdf)
if(l_mrms_sparse_netcdf) then
call read_dbz_mrms_sparse_netcdf(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
string='READ_dbz_mrms_sparse_netcdf'
else
call read_dbz_mrms_netcdf(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
string='READ_dbz_mrms_netcdf'
endif
end if
! Process lagrangian data
else if (obstype == 'lag') then
call read_lag(nread,npuse,nouse,infile,lunout,obstype,&
&twind,gstime,sis)
string='READ_LAG'
! Process lidar winds
else if (obstype == 'dw') then
call read_lidar(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
nobs_sub1(1,i))
string='READ_LIDAR'
! Process synthetic tc-mslp obs
else if (obstype == 'tcp') then
call read_tcps(nread,npuse,nouse,infile,obstype,lunout,sis, &
nobs_sub1(1,i))
string='READ_TCPS'
else if (obstype == 'pm2_5' .or. obstype == 'pm10') then
if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
call oneobschem(nread,npuse,nouse,gstime,&
&infile,obstype,lunout,sis,nobs_sub1(1,i))
string='ONEOBSCHEM'
else
call read_anowbufr(nread,npuse,nouse,gstime,&
&infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
string='READ_ANOWBUFR'
endif
! Process pblh
else if (obstype == 'pblh') then
call read_pblh(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
nobs_sub1(1,i))
string='READ_PBLH'
end if conv_obstype_select
! Process swcp and lwcp
else if (ditype(i) == 'wcp') then
if ( obstype == 'swcp' .or. obstype == 'lwcp' ) then
call read_wcpbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
prsl_full,nobs_sub1(1,i),read_rec(i))
string='READ_WCPBUFR'
end if
else if (ditype(i) == 'rad')then
call radiance_obstype_search(obstype,radmod)
! Process TOVS 1b data
rad_obstype_select: &
if (platid /= 'aqua' .and. (obstype == 'amsua' .or. &
obstype == 'amsub' .or. obstype == 'msu' .or. &
obstype == 'mhs' .or. obstype == 'hirs4' .or. &
obstype == 'hirs3' .or. obstype == 'hirs2' .or. &
obstype == 'ssu' )) then
call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), nobs_sub1(1,i), &
read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
string='READ_BUFRTOVS'
! Process atms data
else if (obstype == 'atms') then
call read_atms(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i),&
read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
string='READ_ATMS'
! Process saphir data
else if (obstype == 'saphir') then
call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),dval_use)
string='READ_SAPHIR'
! Process airs data
else if(platid == 'aqua' .and. (obstype == 'airs' .or. &
obstype == 'amsua' .or. obstype == 'hsb' ))then
call read_airs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_AIRS'
! Process iasi data
else if(obstype == 'iasi')then
call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
string='READ_IASI'
! Process cris data
else if(obstype == 'cris' .or. obstype =='cris-fsr' )then
call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
string='READ_CRIS'
! Process GOES sounder data
! Process raw or prepbufr files (1x1 or 5x5)
else if (obstype == 'sndr' .or. &
obstype == 'sndrd1' .or. obstype == 'sndrd2' .or. &
obstype == 'sndrd3' .or. obstype == 'sndrd4') then
call read_goesndr(mype,val_dat,ithin,rmesh,platid,&
infile,lunout,obstype,nread,npuse,nouse,twind,gstime,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_GOESNDR'
! Process ssmi data
else if (obstype == 'ssmi' ) then
call read_ssmi(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_SSMI'
! Process amsre data
else if ( obstype == 'amsre_low' .or. obstype == 'amsre_mid' .or. &
obstype == 'amsre_hig' ) then
call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_AMSRE'
! Process ssmis data
else if (obstype == 'ssmis' .or. &
obstype == 'ssmis_las' .or. obstype == 'ssmis_uas' .or. &
obstype == 'ssmis_img' .or. obstype == 'ssmis_env' ) then
call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_SSMIS'
! Process AMSR2 data
else if(obstype == 'amsr2')then
call read_amsr2(mype,val_dat,ithin,rmesh,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i))
string='READ_AMSR2'
! Process GOES IMAGER RADIANCE data
else if(obstype == 'goes_img') then
call read_goesimg(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_GOESMIMG'
! Process GMI data
else if (obstype == 'gmi') then
call read_gmi(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),dval_use)
string='READ_GMI'
! Process Meteosat SEVIRI RADIANCE data
else if(obstype == 'seviri') then
call read_seviri(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_SEVIRI'
! Process GOES-R ABI RADIANCE data
else if(obstype == 'abi') then
call read_abi(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_ABI'
! Process Himawari-8 AHI RADIANCE data
else if(obstype == 'ahi') then
call read_ahi(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i))
string='READ_AHI'
! Process NAVY AVHRR RADIANCE data
else if(obstype == 'avhrr_navy') then
call read_avhrr_navy(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_AVH_NAVY'
! Process NESDIS AVHRR RADIANCE data
else if(obstype == 'avhrr') then
call read_avhrr(mype,val_dat,ithin,rmesh,platid,gstime,&
infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i),read_rec(i),dval_use)
string='READ_AVHRR'
end if rad_obstype_select
! Process ozone data
else if (ditype(i) == 'ozone')then
ozone_obstype_select: &
if (is_extOzone(infile,obstype,dplat(i))) then
call extOzone_read(infile,obstype,dplat(i),dsis(i), &
iread,ipuse,iouse, platid,gstime,lunout,twind,ithin,rmesh, &
nobs_sub1(:,i))
string='extOzone_read'
nread=nread+iread
npuse=npuse+ipuse
nouse=nouse+iouse
else
call read_ozone(nread,npuse,nouse,&
platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
nobs_sub1(1,i))
string='READ_OZONE'
endif ozone_obstype_select
! Process co data
else if (ditype(i) =='co')then
call read_co(nread,npuse,nouse,infile,gstime,lunout,obstype,sis, &
nobs_sub1(1,i))
string='READ_CO'
! Process precipitation
else if (ditype(i) == 'pcp')then
call read_pcp(nread,npuse,nouse,gstime,infile,lunout,obstype,twind,sis, &
nobs_sub1(1,i))
string='READ_PCP'
! Process gps observations
else if (ditype(i) == 'gps')then
call read_gps(nread,npuse,nouse,infile,lunout,obstype,twind, &
nprof_gps1,sis,nobs_sub1(1,i))
string='READ_GPS'
! Process aerosol data
else if (ditype(i) == 'aero' )then
call read_aerosol(nread,npuse,nouse,&
platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
nobs_sub1(1,i))
string='READ_AEROSOL'
! Process satellite lightning observations (e.g. GOES/GLM)
else if(ditype(i) == 'light')then
if (obstype == 'goes_glm' ) then
call read_goesglm(nread,ndata,nodata,infile,obstype,lunout,sis)
string='READ_GOESGLM'
endif
end if ditype_select
! Close unit to data file
! Accumulate data counts on "root" task
if (mype_sub(mm1,i)==mype_root) then
close(lunout)
ndata1(i,1)=ndata1(i,1)+npuse
ndata1(i,2)=ndata1(i,2)+nread
ndata1(i,3)=ndata1(i,3)+nouse
if(string(1:2)=='--') then
call warn('read_obs','task without reader, i =',i)
call warn('read_obs',' infile =',trim(infile))
call warn('read_obs',' ditype(i) =',trim(ditype(i)))
call warn('read_obs',' obstype =',trim(obstype))
call warn('read_obs',' platid =',trim(platid))
call warn('read_obs',' string =',trim(string))
endif
write(6,8000) adjustl(string),infile,obstype,sis,nread,ithin,&
rmesh,isfcalc,nouse,npe_sub(i)
8000 format(1x,a22,': file=',a15,&
' type=',a10, ' sis=',a20, ' nread=',i10,&
' ithin=',i2, ' rmesh=',f11.6,' isfcalc=',i2,&
' nkeep=',i10,' ntask=',i3)
endif
endif task_belongs
end do loop_of_obsdata_files
deallocate(prsl_full)
deallocate(hgtl_full)
! Broadcast aircraft new tail numbers for aircraft
! temperature bias correction
! if (aircraft_t_bc) then
! call mpi_barrier(mpi_comm_world,ierror)
! call mpi_bcast(ntail_update,1,mpi_itype,mype_airobst,mpi_comm_world,ierror)
! call mpi_bcast(idx_tail,max_tail,mpi_itype,mype_airobst,mpi_comm_world,ierror)
! call mpi_bcast(taillist,max_tail,MPI_CHARACTER,mype_airobst,mpi_comm_world,ierror)
! end if
! Deallocate arrays containing full horizontal surface fields
call destroy_sfc
! Deallocate arrays containing full horizontal nsst fields
if (nst_gsi > 0) call gsi_nstcoupler_final()
! Sum and distribute number of obs read and used for each input ob group
call mpi_allreduce(ndata1,ndata,ndat*3,mpi_integer,mpi_sum,mpi_comm_world,&
ierror)
! Collect super obs factors
if(mype == 0)write(6,*) ' dval_use = ',dval_use
if(dval_use)then
call mpi_allreduce(super_val,super_val1,superp+1,mpi_rtype,&
mpi_sum,mpi_comm_world,ierror)
else
super_val1=zero
end if
super_val1(0)=one
deallocate(super_val)
! Collect number of gps profiles (needed later for qc)
call mpi_allreduce(nprof_gps1,nprof_gps,1,mpi_integer,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(nobs_sub1,nobs_sub,npe*ndat,mpi_integer,mpi_sum,mpi_comm_world,&
ierror)
! Write collective obs selection information to scratch file.
if (lread_obs_save .and. mype==0) then
write(6,*)'READ_OBS: write collective obs selection info to ',trim(obs_input_common)
call unformatted_open(lunsave,file=obs_input_common,class=".obs_input.")
write(lunsave) ndata,ndat,npe,superp,nprof_gps,ditype
write(lunsave) super_val1
write(lunsave) nobs_sub
close(lunsave)
endif
! End of routine
return
end subroutine read_obs
end module read_obsmod