module m_gpsStats
!$$$ subprogram documentation block
! . . . .
! subprogram: module m_gpsStats
! prgmmr: j guo <jguo@nasa.gov>
! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3
! date: 2015-08-28
!
! abstract: a modular wrapper of (gsp_allhead, gps_alltail), and genstats_gps()
!
! program history log:
! 2007-06-22 cucurull - modify gps_all_ob_type structure
! 2015-08-28 j guo - created this module on top of genstats_gps();
! . completed with type/data components from obsmod;
! . changed code where this module needs to be used;
! . added this document block;
! . for earlier history log, see the history log section
! inside ::genstats_gps() below.
! 2016-05-18 j guo - Made the type private, since this is only a single
! instance module object, with its components defined
! as module variables (gps_allhead, and gsp_alltail).
! . Removed old interface names, which are no longer used
! in this version of GSI.
! . Edited the in-file documentation.
!
! input argument list: see Fortran 90 style document below
!
! output argument list: see Fortran 90 style document below
!
! attributes:
! language: Fortran 90 and/or above
! machine:
!
!$$$ end subprogram documentation block
! module interface:
use m_gpsNode, only: gps_ob_type => gpsNode
use kinds , only: r_kind,i_kind
implicit none
private ! except
! Data Structure:
!public:: gps_all_ob_type ! currently not required to be public
type gps_all_ob_type
type(gps_all_ob_type),pointer :: llpoint => NULL()
type(gps_ob_type),pointer :: mmpoint => NULL()
real(r_kind) :: ratio_err
real(r_kind) :: obserr
real(r_kind) :: dataerr
real(r_kind) :: pg
real(r_kind) :: b
real(r_kind) :: loc
real(r_kind) :: type
real(r_kind),dimension(:),pointer :: rdiag => NULL()
integer(i_kind) :: kprof
logical :: luse ! flag indicating if ob is used in pen.
logical :: muse ! flag indicating if ob is used in pen.
character(8) :: cdiag
integer(i_kind) :: idv,iob ! device id and obs index for sorting
real (r_kind) :: elat, elon ! earth lat-lon for redistribution
!real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution
end type gps_all_ob_type
type gps_all_ob_head
integer(i_kind):: n_alloc=0
type(gps_all_ob_type),pointer :: head => NULL()
end type gps_all_ob_head
! Data objects:
public:: gps_allhead
public:: gps_alltail
! interfaces:
public:: gpsStats_create
public:: gpsStats_destroy
interface gpsStats_create ; module procedure create_; end interface
interface gpsStats_destroy; module procedure destroy_genstats_gps; end interface
public:: gpsStats_genStats
interface gpsStats_genStats; module procedure genstats_gps; end interface
! Synopsis:
! - []_create: allocated in gsimod::gsimain_initialize(). It was
! done through ::create_obsmod_vars(), and now next to it.
!
! - externally built, node-by-node, in setupbend() or setupref().
! This is the reason why the ADT can not be defined as "private".
!
! - []_genStats: used to update m_rhs::[ab]work, in setuprhsall(),
! through ::genstats_gps();
!
! - []_destroy: deallocated within genstats_gps(), when it it
! finished.
!
! The use of []_create()/[]_destroy() pair is obviously not symmetric. It
! would cleaner if they are be moved to the level of setuprhsall(),
! where m_rhs::[ab]work are computed. e.g.,
!
! if(init_pass) call gpsStats_create()
! ...
! if(last_pass) then
! call gpsStats_genstats(bwork,awork,...)
! call gpsStats_destroy()
! endif
!
! As it is now, the second half has been done, but the first half
! stayed at where it has been, for later.
! Most implementations of this module, are snap-shots of gps_all_ob_type, from
! obsmod, its original home. These implementations include, the type
! definition (gps_all_ob_type and gsp_app_ob_head), instantiation of the type
! (gps_allhead and gsp_alltail), and interfaces for the operations of this
! object ([]_create, []_genstats, []_destroy). The part of implementation for
! the node-growing, is remained in setupbend() and setupref() as is.
type(gps_all_ob_head),dimension(:),pointer :: gps_allhead => null()
type(gps_all_ob_head),dimension(:),pointer :: gps_alltail => null()
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
character(len=*),parameter :: myname='genstats_gps'
contains
subroutine create_()
use gsi_4dvar, only: nobs_bins
implicit none
ALLOCATE(gps_allhead(nobs_bins))
ALLOCATE(gps_alltail(nobs_bins))
end subroutine create_
subroutine destroy_genstats_gps()
!$$$ subprogram documentation block
! . . . .
! subprogram: destroy_genstats_gps
! prgmmr: treadon org: np20 date: 2005-12-21
!
! abstract: deallocate arrays holding gps information
!
! program history log:
! 2005-12-21 treadon
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
use gsi_4dvar, only: nobs_bins
use kinds, only: i_kind
implicit none
integer(i_kind):: istatus,ii
do ii=1,nobs_bins
gps_alltail(ii)%head => gps_allhead(ii)%head
do while (associated(gps_alltail(ii)%head))
gps_allhead(ii)%head => gps_alltail(ii)%head%llpoint
deallocate(gps_alltail(ii)%head,stat=istatus)
if (istatus/=0) write(6,*)'DESTROY_GENSTATS_GPS: deallocate error for gps_all, istatus=',istatus
gps_alltail(ii)%head => gps_allhead(ii)%head
end do
end do
return
end subroutine destroy_genstats_gps
subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: genstats_gps generate statistics for gps observations
! prgmmr: treadon org: np20 date: 2005-12-21
!
! abstract: For gps observations, this routine
! a) collects statistics for runtime diagnostic output
! f) adjusts observation error ratio based on superobs factor
!
! program history log:
! 2005-12-21 treadon
! 2006-02-24 derber - modify to take advantage of convinfo module
! 2006-07-28 derber - modify to use new inner loop obs data structure
! 2006-09-20 cucurull - replace superobs factor for obs in a top (non-full) layer
! 2007-03-01 treadon - add array toss_gps
! 2007-03-19 tremolet - binning of observations
! 2007-06-21 cucurull - add conv_diagsave and mype in argument list;
! modify qc and output for diagnostic file based on toss_gps
! print out diagnostic files if requested
! add wgtlim and huge_single in constants module
! 2008-02-27 cucurull - modify diagnostics output file
! 2008-04-14 treadon - compute super_gps within this routine
! 2008-06-04 safford - rm unused vars and uses
! 2008-09-05 lueken - merged ed's changes into q1fy09 code
! 2008-25-08 todling - adapt obs-binning change to GSI-May2008
! 2009-02-05 cucurull - modify latitude range four statistics output
! 2009-10-22 shen - add high_gps
! 2010-04-09 cucurull - fix several bugs for high_gps (diag information, counters,
! - consider failure of gross check, obs-binning structures, QC for CL profiles)
! - reorganize high_gps structure
! - modify dimension of diagnostic structure
! 2010-07-23 treadon - add ratio_error=zero to reqional QC block, replace (izero,ione) with (0,1),
! remove _i_kind suffix from integer constants, clean up use statements
! 2010-08-17 treadon - convert high_gps from m to km one time only; break out regional
! QC as separate if/then block (global will bypass); replace
! ratio_errors_reg with logical toss
! 2010-10-25 cucurull - add quality control options for C/NOFS satellite
! 2011-01-18 cucurull - increase the size of nreal and mreal by one element to
! add gps_dtype information
! 2012-10-16 cucurull - increase the size of nreal and mreal by one element to
! add qrefges information, replace qcfail=5 by 4, add regional QC for MetOpB
! add dtype, dobs to distinguish use of toss_gps between ref/bending, add SR QC for obs
! 2014-01-28 todling - write sensitivity slot indicator (ioff) to header of diagfile
! 2014-12-13 derber - minor optimization modifications
! 2015-07-28 cucurull - add QC for regional bending angle assimilation
! 2015-08-28 guo - wrapped as a module (m_gpsStats)
! moved the call to obsmod::destroy_genstats_gps() to
! where this routine was used (setuprhsall()), with its
! new module interface name. gpsStats_destroy().
! 2016-11-29 shlyaeva - increase the size of nreal for saving linearized Hx for EnKF
! 2016-12-09 mccarty - add ncdiag writing support
!
! input argument list:
! toss_gps_sub - array of qc'd profile heights
! conv_diagsave - logical to save innovation diagnostics
! mype - mpi task id
!
! output argument list:
! bwork - array containing information about obs-ges statistics
! awork - array containing information for data counts and gross checks
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use obsmod, only: nprof_gps,lobsdiag_forenkf
use obsmod, only: obs_diag,lobsdiagsave,luse_obsdiag
use obsmod, only: binary_diag,netcdf_diag,dirname,ianldate
use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
nc_diag_write, nc_diag_data2d
use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close
use gridmod, only: nsig,regional
use constants, only: tiny_r_kind,half,wgtlim,one,two,zero,five,four
use qcmod, only: npres_print,ptop,pbot
use mpimod, only: ierror,mpi_comm_world,mpi_rtype,mpi_sum,mpi_max
use jfunc, only: jiter,miter,jiterstart
use gsi_4dvar, only: nobs_bins
use convinfo, only: nconvtype
use state_vectors, only: nsdim
implicit none
! Declare passed variables
logical ,intent(in):: conv_diagsave
integer(i_kind) ,intent(in) :: mype
real(r_kind),dimension(100+7*nsig) ,intent(inout):: awork
real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout):: bwork
real(r_kind),dimension(max(1,nprof_gps)) ,intent(in):: toss_gps_sub
! Declare local parameters
real(r_kind),parameter:: ten = 10.0_r_kind
real(r_kind),parameter:: six = 6.0_r_kind
real(r_kind),parameter:: r1em3 = 1.0e-3_r_kind
real(r_kind),parameter:: r20 = 20.0_r_kind
real(r_kind),parameter:: scale = 100.0_r_kind
! Declare local variables
logical:: luse,muse,toss,save_jacobian
integer(i_kind):: k,jsig,icnt,khgt,kprof,ikx,nn,j,nchar,nreal,mreal,ii,ioff
real(r_kind):: pressure,arg,wgross,wgt,term,cg_gps,valqc,elev,satid,dtype,dobs
real(r_kind):: ress,val,ratio_errors,val2
real(r_kind):: exp_arg,data_ikx,data_rinc,cg_term,rat_err2,elat
real(r_kind):: wnotgross,data_ipg,data_ier,data_ib,factor,super_gps_up,rhgt
real(r_kind),dimension(nsig,max(1,nprof_gps)):: super_gps_sub,super_gps
real(r_kind),dimension(max(1,nprof_gps)):: toss_gps
real(r_kind),dimension(max(1,nprof_gps)):: high_gps,high_gps_sub
real(r_kind),dimension(max(1,nprof_gps)):: dobs_height,dobs_height_sub
real(r_single),allocatable,dimension(:,:)::sdiag
character(8),allocatable,dimension(:):: cdiag
real(r_single), dimension(nsdim) :: dhx_dx_array
type(obs_diag), pointer :: obsptr => NULL()
integer(i_kind) :: nnz, nind
type(gps_ob_type), pointer:: gpsptr
type(gps_all_ob_type), pointer:: gps_allptr
save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf
!*******************************************************************************
! Check to see if there are any profiles to process. If none, return.
if (nprof_gps==0) then
if (mype==0) write(6,*)'GENSTATS_GPS: no profiles to process (nprof_gfs=',nprof_gps,'), EXIT routine'
return
endif
! Reduce sub-domain specific QC'd profile height cutoff values to
! maximum global value for each profile
toss_gps=zero
call mpi_allreduce(toss_gps_sub,toss_gps,nprof_gps,mpi_rtype,mpi_max,&
mpi_comm_world,ierror)
! If netcdf diag, initialize it
if (conv_diagsave.and.netcdf_diag) call init_netcdf_diag_
! Get height of maximum bending angle
dobs_height_sub = zero
DO ii=1,nobs_bins
gps_allptr => gps_allhead(ii)%head
do while (associated(gps_allptr))
! Load local work variables
kprof = gps_allptr%kprof
dtype = gps_allptr%rdiag(20)
dobs = gps_allptr%rdiag(17)
if (dtype == one .and. toss_gps(kprof) > zero .and. dobs == toss_gps(kprof)) then
dobs_height_sub(kprof) = gps_allptr%rdiag(7)
endif
gps_allptr => gps_allptr%llpoint
! End loop over observations
end do
! End of loop over time bins
END DO
! Reduce sub-domain specific QC'd profile height to maximum global value for each profile
dobs_height=zero
call mpi_allreduce(dobs_height_sub,dobs_height,nprof_gps,mpi_rtype,mpi_max,&
mpi_comm_world,ierror)
! Compute superobs factor on sub-domains using global QC'd profile height
super_gps_sub=zero
high_gps_sub = zero
DO ii=1,nobs_bins
gps_allptr => gps_allhead(ii)%head
do while (associated(gps_allptr))
! Load local work variables
ratio_errors = gps_allptr%ratio_err
data_ier = gps_allptr%obserr
luse = gps_allptr%luse
kprof = gps_allptr%kprof
dtype = gps_allptr%rdiag(20)
! Accumulate superobs factors and get highest good gps obs within a profile
if (dtype == zero) then ! refractivity
rhgt = gps_allptr%loc
if (rhgt >toss_gps(kprof)) then
if(ratio_errors*data_ier>tiny_r_kind) then
elev = gps_allptr%rdiag(7)
high_gps_sub(kprof)=max(high_gps_sub(kprof),elev)
if(luse) then
khgt = gps_allptr%loc
k=min(max(1,khgt),nsig)
super_gps_sub(k,kprof)=super_gps_sub(k,kprof)+one
endif
endif
endif
else ! bending angle
dobs = gps_allptr%rdiag(17)
if(toss_gps(kprof) == zero .or. (toss_gps(kprof) > zero .and. dobs < toss_gps(kprof))) then ! will not fail SR from obs qc
elev = gps_allptr%rdiag(7)
if(elev > dobs_height(kprof)) then
if(ratio_errors*data_ier>tiny_r_kind) then
high_gps_sub(kprof)=max(high_gps_sub(kprof),elev)
if(luse) then
khgt = gps_allptr%loc
k=min(max(1,khgt),nsig)
super_gps_sub(k,kprof)=super_gps_sub(k,kprof)+one
endif
endif
endif
endif
endif
gps_allptr => gps_allptr%llpoint
! End loop over observations
end do
! End of loop over time bins
END DO
super_gps = zero
high_gps = zero
! Reduce sub-domain specifc superobs factors to global values for each profile
call mpi_allreduce(super_gps_sub,super_gps,nsig*nprof_gps,mpi_rtype,mpi_sum,&
mpi_comm_world,ierror)
! Reduce sub-domain specific high_gps values to global values for each profile
call mpi_allreduce(high_gps_sub,high_gps,nprof_gps,mpi_rtype,mpi_max,&
mpi_comm_world,ierror)
! Convert high_gps from meters to kilometers
high_gps = r1em3*high_gps
! If generating diagnostic output, need to determine dimension of output arrays.
nreal=0
ioff =nreal
if (conv_diagsave) then
icnt = zero
DO ii=1,nobs_bins
gps_allptr => gps_allhead(ii)%head
do while (associated(gps_allptr))
luse = gps_allptr%luse
if(luse)icnt=icnt+1
gps_allptr => gps_allptr%llpoint
end do
END DO
if(icnt > 0)then
nreal =22
ioff =nreal
if (lobsdiagsave) nreal=nreal+4*miter+1
if (save_jacobian) then
nnz = 3*nsig
nind = 3
nreal = nreal + 2*nind + nnz + 2
endif
allocate(cdiag(icnt),sdiag(nreal,icnt))
end if
endif
! Loop over data to apply final qc, superobs factors, accumulate
! statistics and (optionally) load diagnostic output arrays
icnt=0
DO ii=1,nobs_bins
gps_allptr => gps_allhead(ii)%head
do while (associated(gps_allptr))
! Load local work variables
ratio_errors = gps_allptr%ratio_err
data_ier = gps_allptr%obserr
luse = gps_allptr%luse
muse = gps_allptr%muse
khgt = gps_allptr%loc
kprof = gps_allptr%kprof
dtype = gps_allptr%rdiag(20)
gpsptr => gps_allptr%mmpoint
if(muse .and. associated(gpsptr) .and. luse_obsdiag)then
obsptr => gpsptr%diags
endif
! Determine model level to which observation is mapped to
k=min(max(1,khgt),nsig)
! Normalize ratio_errors by superobs factor. Update ratio_error
! term used in minimization
super_gps_up=zero
if (super_gps(k,kprof)>tiny_r_kind) then
do j=min(k+1,nsig),nsig
super_gps_up = max(super_gps_up,super_gps(j,kprof))
enddo
if (super_gps_up >tiny_r_kind) then
factor = one / sqrt(super_gps(k,kprof))
else
factor = one / sqrt(max(super_gps(k-1,kprof),super_gps(k,kprof)))
endif
ratio_errors = ratio_errors * factor
if(conv_diagsave .and. luse) then
if(gps_allptr%rdiag(16) >tiny_r_kind) gps_allptr%rdiag(16)=ratio_errors*data_ier
endif
! Adjust error ratio for observations used in inner loop
if (associated(gpsptr)) then
gpsptr%raterr2 = ratio_errors **2
if(associated(obsptr) .and. luse_obsdiag)then
obsptr%wgtjo=(ratio_errors*data_ier)**2
end if
endif
endif
! For given profile, check if observation level is below level at
! which profile data is tossed. If so, set error parameter to
! zero (effectively tossing the obs).
rhgt = gps_allptr%loc
mreal = 22
if(dtype == zero) then !refractivity
if (rhgt<=toss_gps(kprof)) then
if(ratio_errors*data_ier > tiny_r_kind) then ! obs was good
if (luse) then
if(conv_diagsave) then
gps_allptr%rdiag(10) = four
gps_allptr%rdiag(12) = -one
gps_allptr%rdiag(16) = zero
if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
endif
elat = gps_allptr%rdiag(3)
if(elat > r20) then
awork(22) = awork(22)+one
else if(elat< -r20)then
awork(23) = awork(23)+one
else
awork(24) = awork(24)+one
end if
endif
endif
ratio_errors = zero
if (associated(gpsptr)) then
gpsptr%raterr2 = ratio_errors **2
if(associated(obsptr) .and. luse_obsdiag)then
obsptr%wgtjo=zero
obsptr%muse(jiter)=.false.
end if
endif
endif
else
elev = gps_allptr%rdiag(7)
dobs = gps_allptr%rdiag(17)
if (toss_gps(kprof) > zero .and. (dobs == toss_gps(kprof) .or. elev < dobs_height(kprof))) then ! SR from obs
if(ratio_errors*data_ier > tiny_r_kind) then ! obs was good
if (luse) then
if(conv_diagsave) then
gps_allptr%rdiag(10) = four
gps_allptr%rdiag(12) = -one
gps_allptr%rdiag(16) = zero
if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
endif
elat = gps_allptr%rdiag(3)
if(elat > r20) then
awork(22) = awork(22)+one
else if(elat< -r20)then
awork(23) = awork(23)+one
else
awork(24) = awork(24)+one
end if
endif
endif
ratio_errors = zero
if (associated(gpsptr)) then
gpsptr%raterr2 = ratio_errors **2
if(associated(obsptr) .and. luse_obsdiag)then
obsptr%wgtjo=zero
obsptr%muse(jiter)=.false.
end if
endif
endif
endif
! Regional QC. Remove obs if highest good obs in
! profile is below platform specific threshold height.
if(regional) then
toss=.false.
if(ratio_errors*data_ier > tiny_r_kind) then
if(dtype==zero) then !refractivity
satid = gps_allptr%rdiag(1)
if((satid==41).or.(satid==722).or.(satid==723).or.(satid==4).or.(satid==786).or.(satid==3)) then
if ((high_gps(kprof)) < ten) toss=.true.
else ! OL
if ((high_gps(kprof)) < five) toss=.true.
endif
else !bending angle
if ((high_gps(kprof)) <= six) toss=.true.
endif
endif
if (toss) then
if (luse) then
if(conv_diagsave) then
gps_allptr%rdiag(10) = four
gps_allptr%rdiag(12) = -one
gps_allptr%rdiag(16) = zero
if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
endif
elat = gps_allptr%rdiag(3)
if(elat > r20) then
awork(22) = awork(22)+one
else if(elat< -r20)then
awork(23) = awork(23)+one
else
awork(24) = awork(24)+one
end if
end if
ratio_errors = zero
if (associated(gpsptr)) then
gpsptr%raterr2 = ratio_errors **2
if(associated(obsptr) .and. luse_obsdiag)then
obsptr%wgtjo=zero
obsptr%muse(jiter)=.false.
end if
endif
endif
endif ! regional
! Compute penalty terms
if (ratio_errors*data_ier <= tiny_r_kind) muse = .false.
if(luse)then
val = gps_allptr%dataerr
val2 = val*val
exp_arg = -half*val2
rat_err2 = ratio_errors**2
data_ipg = gps_allptr%pg
if (data_ipg > tiny_r_kind) then
data_ib = gps_allptr%b
cg_gps=cg_term/data_ib
wnotgross= one-data_ipg
wgross = data_ipg*cg_gps
arg = exp(exp_arg)
term = log(wnotgross*arg+wgross)
wgt = wnotgross*arg/(wnotgross*arg+wgross)
else
term = exp_arg
wgt = one
endif
if(conv_diagsave) gps_allptr%rdiag(13) = wgt/wgtlim
valqc = -two*rat_err2*term
! Accumulate statistics for obs belonging to this task
! based on interface (not mid-point) level
val2=val2*rat_err2
if(muse)then
if(wgt < wgtlim) awork(21) = awork(21)+one
! Accumulate values for penalty and data count
jsig=max(1,khgt)
awork(jsig+3*nsig+100)=awork(jsig+3*nsig+100)+valqc
awork(jsig+5*nsig+100)=awork(jsig+5*nsig+100)+one
awork(jsig+6*nsig+100)=awork(jsig+6*nsig+100)+val2
nn=1
else
nn=2
if(ratio_errors*data_ier >=tiny_r_kind)nn=3
endif
data_ikx = gps_allptr%type
ikx = nint(data_ikx)
pressure = gps_allptr%rdiag(6)
data_rinc = gps_allptr%rdiag(5)*scale
! Loop over pressure level groupings and obs to accumulate
! statistics as a function of observation type.
do k = 1,npres_print
if(pressure>ptop(k) .and. pressure<=pbot(k))then
ress=data_rinc
bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count
bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+ress ! (o-g)
bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ress*ress ! (o-g)**2
bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2 ! penalty
bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty
end if
end do
end if
! Transfer diagnostic information to output arrays
if(conv_diagsave .and. luse) then
icnt=icnt+1
cdiag(icnt) = gps_allptr%cdiag
do j=1,nreal
sdiag(j,icnt)= gps_allptr%rdiag(j)
enddo
endif
if (conv_diagsave .and. netcdf_diag .and. luse) call contents_netcdf_diag_
gps_allptr => gps_allptr%llpoint
! End loop over observations
end do
! End of loop over time bins
END DO
! If requested, write information to diagnostic file
if(conv_diagsave) then
if (netcdf_diag) call nc_diag_write
if (binary_diag .and. icnt > 0)then
nchar = 1
write(7)'gps',nchar,nreal,icnt,mype,ioff
write(7)cdiag,sdiag
deallocate(cdiag,sdiag)
endif
endif
! Destroy arrays holding gps data
call destroy_genstats_gps
contains
subroutine init_netcdf_diag_
integer(i_kind) ncd_fileid, ncd_nobs
character(len=80) string
character(len=128) diag_conv_file
logical append_diag
logical,parameter::verbose=.false.
write(string,900) jiter
900 format('conv_gps_',i2.2,'.nc4')
diag_conv_file=trim(dirname) // trim(string)
inquire(file=diag_conv_file, exist=append_diag)
if (append_diag) then
call nc_diag_read_init(diag_conv_file,ncd_fileid)
ncd_nobs = nc_diag_read_get_dim(ncd_fileid,'nobs')
call nc_diag_read_close(diag_conv_file)
if (ncd_nobs > 0) then
if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists. Appending. nobs,mype=',ncd_nobs,mype
else
if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs. Not appending. nobs,mype=',ncd_nobs,mype
append_diag = .false. ! if there are no obs in existing file, then do not try to append
endif
end if
call nc_diag_init(diag_conv_file, append=append_diag)
if (.not. append_diag) then ! don't write headers on append - the module will break?
call nc_diag_header("date_time",ianldate )
call nc_diag_header("Number_of_state_vars", nsdim )
endif
end subroutine init_netcdf_diag_
subroutine contents_binary_diag_
end subroutine contents_binary_diag_
subroutine contents_netcdf_diag_
use sparsearr, only: sparr2, readarray, fullarray
integer(i_kind),dimension(miter) :: obsdiag_iuse
integer(i_kind) :: obstype, obssubtype
type(sparr2) :: dhx_dx
! Observation class
character(7),parameter :: obsclass = ' gps'
call nc_diag_metadata("Station_ID", gps_allptr%cdiag )
call nc_diag_metadata("Observation_Class", obsclass )
obstype = gps_allptr%rdiag(1)
obssubtype = gps_allptr%rdiag(2)
call nc_diag_metadata("Observation_Type", obstype )
call nc_diag_metadata("Observation_Subtype", obssubtype )
call nc_diag_metadata("Latitude", sngl(gps_allptr%rdiag(3)) )
call nc_diag_metadata("Longitude", sngl(gps_allptr%rdiag(4)) )
call nc_diag_metadata("Incremental_Bending_Angle", sngl(gps_allptr%rdiag(5)) )
call nc_diag_metadata("Pressure", sngl(gps_allptr%rdiag(6)) )
call nc_diag_metadata("Height", sngl(gps_allptr%rdiag(7)) )
call nc_diag_metadata("Time", sngl(gps_allptr%rdiag(8)) )
call nc_diag_metadata("Model_Elevation", sngl(gps_allptr%rdiag(9)) )
call nc_diag_metadata("Setup_QC_Mark", sngl(gps_allptr%rdiag(10)) )
call nc_diag_metadata("Prep_Use_Flag", sngl(gps_allptr%rdiag(11)) )
call nc_diag_metadata("Analysis_Use_Flag", sngl(gps_allptr%rdiag(12)) )
call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", sngl(gps_allptr%rdiag(13)) )
call nc_diag_metadata("Errinv_Input", sngl(gps_allptr%rdiag(14)) )
call nc_diag_metadata("Errinv_Adjust", sngl(gps_allptr%rdiag(15)) )
call nc_diag_metadata("Errinv_Final", sngl(gps_allptr%rdiag(16)) )
call nc_diag_metadata("Observation", sngl(gps_allptr%rdiag(17)) )
call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(gps_allptr%rdiag(17))*sngl(gps_allptr%rdiag(5)) )
call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(gps_allptr%rdiag(17))*sngl(gps_allptr%rdiag(5)) )
call nc_diag_metadata("GPS_Type", sngl(gps_allptr%rdiag(20)) )
call nc_diag_metadata("Temperature_at_Obs_Location", sngl(gps_allptr%rdiag(18)) )
call nc_diag_metadata("Specific_Humidity_at_Obs_Location", sngl(gps_allptr%rdiag(21)) )
if (save_jacobian) then
call readarray(dhx_dx, gps_allptr%rdiag(ioff+1:nreal))
call fullarray(dhx_dx, dhx_dx_array)
call nc_diag_data2d("Observation_Operator_Jacobian", dhx_dx_array)
endif
! call nc_diag_data2d("T_Jacobian", gps_allptr%mmpoint%jac_t )
if (lobsdiagsave) then
print *,'ERROR: OBSDIAGSAVE SKIPPED IN NCDIAG DEVELOPMENT. STOPPING.'
call stop2(55)
! do jj=1,miter
! if (gps_allptr%diags%muse(jj)) then
! obsdiag_iuse(jj) = one
! else
! obsdiag_iuse(jj) = -one
! endif
! enddo
!
! call nc_diag_data2d("ObsDiagSave_iuse", obsdiag_iuse )
! call nc_diag_data2d("ObsDiagSave_nldepart", gps_allptr%diags%nldepart )
! call nc_diag_data2d("ObsDiagSave_tldepart", gps_allptr%diags%tldepart )
! call nc_diag_data2d("ObsDiagSave_obssen", gps_allptr%diags%obssen )
endif
end subroutine contents_netcdf_diag_
end subroutine genstats_gps
end module m_gpsStats