diff --git a/src/gsi/convinfo.f90 b/src/gsi/convinfo.f90 index f6ba29c2e1..83c84d3f38 100755 --- a/src/gsi/convinfo.f90 +++ b/src/gsi/convinfo.f90 @@ -155,6 +155,8 @@ subroutine convinfo_read ! 2009-01-22 todling - protect against non-initialized destroy call ! 2010-05-29 todling - interface consistent w/ similar routines ! 2014-07-10 carley - add check to bypass blank lines in convinfo file +! 2023-05-09 s.vetra-carvalho - extended iotype to be 8 characters long to +! allow for more descriptive observations names ! ! input argument list: ! mype - mpi task id @@ -171,7 +173,7 @@ subroutine convinfo_read implicit none character(len=1)cflg - character(len=7) iotype + character(len=8) iotype character(len=140) crecord integer(i_kind) lunin,i,nc,ier,istat integer(i_kind) nlines @@ -188,13 +190,13 @@ subroutine convinfo_read nlines=0 read1: do cflg=' ' - iotype=' ' + iotype=' ' read(lunin,1030,iostat=istat,end=1130)cflg,iotype,crecord -1030 format(a1,a7,2x,a140) +1030 format(a1,a8,1x,a140) if (istat /= 0) exit nlines=nlines+1 if(cflg == '!')cycle - if (cflg==' '.and.iotype==' ') then + if (cflg==' '.and.iotype==' ') then if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',nlines,' skipping!' cycle end if @@ -242,9 +244,9 @@ subroutine convinfo_read do i=1,nlines cflg=' ' - iotype=' ' + iotype=' ' read(lunin,1030)cflg,iotype,crecord - if (cflg==' '.and.iotype==' ') then + if (cflg==' '.and.iotype==' ') then if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',i,' skipping!' cycle end if @@ -282,7 +284,7 @@ subroutine convinfo_read if(print_verbose .and. mype == 0)write(6,1031)ioctype(nc),ictype(nc),icsubtype(nc),icuse(nc),ctwind(nc),ncnumgrp(nc), & ncgroup(nc),ncmiter(nc),cgross(nc),cermax(nc),cermin(nc),cvar_b(nc),cvar_pg(nc), & ithin_conv(nc),rmesh_conv(nc),pmesh_conv(nc),idum,pmot_conv(nc),ptime_conv(nc),index_sub(nc),ibeta(nc),ikapa(nc) -1031 format('READ_CONVINFO: ',a7,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5) +1031 format('READ_CONVINFO: ',a8,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5) enddo close(lunin) diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 95d885e2ee..97dcc3101f 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -218,6 +218,7 @@ gsi_dbzOper.F90 gsi_dwOper.F90 gsi_enscouplermod.f90 gsi_fedOper.F90 +gsi_gnssrspdOper.F90 gsi_gpsbendOper.F90 gsi_gpsrefOper.F90 gsi_gustOper.F90 @@ -276,6 +277,7 @@ intco.f90 intdbz.f90 intfed.f90 intdw.f90 +intgnssrspd.f90 intgps.f90 intgust.f90 inthowv.f90 @@ -342,6 +344,7 @@ m_dwNode.F90 m_extOzone.F90 m_fedNode.F90 m_find.f90 +m_gnssrspdNode.F90 m_gpsNode.F90 m_gpsrhs.F90 m_gsiBiases.f90 @@ -486,6 +489,7 @@ read_files.f90 read_fl_hdob.f90 read_gfs_ozone_for_regional.f90 read_gmi.f90 +read_gnssrspd.f90 read_goesglm.f90 read_goesimg.f90 read_goesimgr_skycover.f90 @@ -537,6 +541,7 @@ setupdbz.f90 setupdbz_lib.f90 setupdw.f90 setupfed.f90 +setupgnssrspd.f90 setupgust.f90 setuphowv.f90 setuplag.f90 @@ -597,6 +602,7 @@ stpco.f90 stpdbz.f90 stpfed.f90 stpdw.f90 +stpgnssrspd.f90 stpgps.f90 stpgust.f90 stphowv.f90 diff --git a/src/gsi/gsi_gnssrspdOper.F90 b/src/gsi/gsi_gnssrspdOper.F90 new file mode 100644 index 0000000000..92f0d0a2a6 --- /dev/null +++ b/src/gsi/gsi_gnssrspdOper.F90 @@ -0,0 +1,161 @@ +module gsi_gnssrspdOper +!$$$ subprogram documentation block +! . . . . +! subprogram: module gsi_gnssrspdOper +! (based on gsi_*Oper.F90 routines by j guo ) +! prgmmr: K. Apodaca +! org: Spire Global, Inc. +! date: 2023-04-28 +! +! abstract: an obOper extension for gnssrspdNode type +! +! program history log: +! 2023-04-21 k apodaca - initial version +! 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 gsi_obOper, only: obOper + use m_gnssrspdNode , only: gnssrspdNode + implicit none + public:: gnssrspdOper ! data stracture + + type,extends(obOper):: gnssrspdOper + contains + procedure,nopass:: mytype + procedure,nopass:: nodeMold + procedure:: setup_ + procedure:: intjo1_ + procedure:: stpjo1_ + end type gnssrspdOper + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*),parameter :: myname='gsi_gnssrspdOper' + type(gnssrspdNode),save,target:: myNodeMold_ + +contains + function mytype(nodetype) + implicit none + character(len=:),allocatable:: mytype + logical,optional, intent(in):: nodetype + mytype="[gnssrspdOper]" + if(present(nodetype)) then + if(nodetype) mytype=myNodeMold_%mytype() + endif + end function mytype + + function nodeMold() + !> %nodeMold() returns a mold of its corresponding obsNode + use m_obsNode, only: obsNode + implicit none + class(obsNode),pointer:: nodeMold + nodeMold => myNodeMold_ + end function nodeMold + + subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) + use gnssrspd_setup, only: setup + use kinds, only: i_kind + use gsi_obOper, only: len_obstype + use gsi_obOper, only: len_isis + + use m_rhs , only: awork => rhs_awork + use m_rhs , only: bwork => rhs_bwork + use m_rhs , only: iwork => i_gnssrspd + + use obsmod , only: write_diag + use convinfo, only: diag_conv + use jfunc , only: jiter + + use mpeu_util, only: die + implicit none + class(gnssrspdOper ), intent(inout):: self + integer(i_kind), intent(in):: lunin + integer(i_kind), intent(in):: mype + integer(i_kind), intent(in):: is + integer(i_kind), intent(in):: nobs + logical , intent(in):: init_pass ! supporting multi-pass setup() + logical , intent(in):: last_pass ! with incremental backgrounds. + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::setup_" + + character(len=len_obstype):: obstype + character(len=len_isis ):: isis + integer(i_kind):: nreal,nchanl,ier,nele + logical:: diagsave + + if(nobs == 0) return + + read(lunin,iostat=ier) obstype,isis,nreal,nchanl + if(ier/=0) call die(myname_,'read(obstype,...), iostat =',ier) + nele = nreal+nchanl + + diagsave = write_diag(jiter) .and. diag_conv + + call setup(self%obsLL(:), self%odiagLL(:), & + lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave) + + end subroutine setup_ + + subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) + use intgnssrspdmod, only: intjo => intgnssrspd + use gsi_bundlemod , only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds , only: i_kind, r_quad + implicit none + class(gnssrspdOper ),intent(in ):: self + integer(i_kind ),intent(in ):: ibin + type(gsi_bundle),intent(inout):: rval ! (ibin) + type(gsi_bundle),intent(in ):: sval ! (ibin) + real(r_quad ),target,dimension(:),intent(inout):: qpred ! (ibin) + type(predictors),target, intent(in ):: sbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::intjo1_" + class(obsNode),pointer:: headNode + + headNode => obsLList_headNode(self%obsLL(ibin)) + call intjo(headNode, rval,sval) + headNode => null() + + end subroutine intjo1_ + + subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) + use stpgnssrspdmod, only: stpjo => stpgnssrspd + use gsi_bundlemod, only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds, only: r_quad,r_kind,i_kind + implicit none + class(gnssrspdOper ),intent(in):: self + integer(i_kind ),intent(in):: ibin + type(gsi_bundle),intent(in):: dval + type(gsi_bundle),intent(in):: xval + real(r_quad ),dimension(:),intent(inout):: pbcjo ! (1:4) + real(r_kind ),dimension(:),intent(in ):: sges + integer(i_kind),intent(in):: nstep + + type(predictors),target, intent(in):: dbias + type(predictors),target, intent(in):: xbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::stpjo1_" + class(obsNode),pointer:: headNode + + headNode => obsLList_headNode(self%obsLL(ibin)) + call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep) + headNode => null() + end subroutine stpjo1_ + +end module gsi_gnssrspdOper diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index 6db7921905..dfcde80a32 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -12,7 +12,8 @@ module gsi_obOperTypeManager ! 2018-07-12 j guo - a type-manager for all obOper extensions. ! - an enum mapping of obsinput::dtype(:) to obOper type ! extensions. -! +! 2022-03-15 k apodaca - add GNSS-R L2 ocean wind speed type-manager +! 2023-03-20 k apodaca - add GNSS-R DDM type-manager ! ! input argument list: see Fortran 90 style document below ! ! output argument list: see Fortran 90 style document below @@ -52,6 +53,7 @@ module gsi_obOperTypeManager use gsi_radOper , only: radOper use gsi_rwOper , only: rwOper use gsi_spdOper , only: spdOper + use gsi_gnssrspdOper , only: gnssrspdOper use gsi_sstOper , only: sstOper use gsi_swcpOper , only: swcpOper use gsi_tcamtOper , only: tcamtOper @@ -102,6 +104,7 @@ module gsi_obOperTypeManager public:: iobOper_w public:: iobOper_q public:: iobOper_spd + public:: iobOper_gnssrspd public:: iobOper_rw public:: iobOper_dw public:: iobOper_sst @@ -148,6 +151,7 @@ module gsi_obOperTypeManager enumerator:: iobOper_w enumerator:: iobOper_q enumerator:: iobOper_spd + enumerator:: iobOper_gnssrspd enumerator:: iobOper_rw enumerator:: iobOper_dw enumerator:: iobOper_sst @@ -210,6 +214,7 @@ module gsi_obOperTypeManager type( wOper), target, save:: wOper_mold type( qOper), target, save:: qOper_mold type( spdOper), target, save:: spdOper_mold + type( gnssrspdOper), target, save:: gnssrspdOper_mold type( rwOper), target, save:: rwOper_mold type( dwOper), target, save:: dwOper_mold type( sstOper), target, save:: sstOper_mold @@ -264,6 +269,7 @@ function dtype2index_(dtype) result(index_) case("q" ,"[qoper]" ); index_= iobOper_q case("spd" ,"[spdoper]" ); index_= iobOper_spd + case("gnssrspd" ,"[gnssrspdoper]" ); index_= iobOper_gnssrspd case("rw" ,"[rwoper]" ); index_= iobOper_rw case("dw" ,"[dwoper]" ); index_= iobOper_dw case("sst" ,"[sstoper]" ); index_= iobOper_sst @@ -293,7 +299,7 @@ function dtype2index_(dtype) result(index_) case("ompslp" ); index_= iobOper_o3l case("ompslpuv" ); index_= iobOper_o3l case("ompslpvis"); index_= iobOper_o3l - case("ompslpnc" ); index_= iobOper_o3l + case("ompslpnc" ); index_= iobOper_o3l case("gpsbend","[gpsbendoper]"); index_= iobOper_gpsbend case("gps_bnd"); index_= iobOper_gpsbend @@ -457,6 +463,7 @@ function index2vmold_(iobOper) result(vmold_) case(iobOper_w ); vmold_ => wOper_mold case(iobOper_q ); vmold_ => qOper_mold case(iobOper_spd ); vmold_ => spdOper_mold + case(iobOper_gnssrspd ); vmold_ => gnssrspdOper_mold case(iobOper_rw ); vmold_ => rwOper_mold case(iobOper_dw ); vmold_ => dwOper_mold case(iobOper_sst ); vmold_ => sstOper_mold @@ -573,6 +580,7 @@ subroutine cobstype_config_() cobstype(iobOper_w ) ="wind " ! w_ob_type cobstype(iobOper_q ) ="moisture " ! q_ob_type cobstype(iobOper_spd ) ="wind speed " ! spd_ob_type + cobstype(iobOper_gnssrspd ) ="gnss-r wind speed " ! gnssrspd_ob_type cobstype(iobOper_rw ) ="radial wind " ! rw_ob_type cobstype(iobOper_dw ) ="doppler wind " ! dw_ob_type cobstype(iobOper_sst ) ="sst " ! sst_ob_type diff --git a/src/gsi/gsimain.f90 b/src/gsi/gsimain.f90 index 5eaa38286a..5c56259bc7 100644 --- a/src/gsi/gsimain.f90 +++ b/src/gsi/gsimain.f90 @@ -131,6 +131,8 @@ program gsi ! related to the GOES/GLM lightnig assimilation ! 2019-07-09 todling - add initialization of abstract layer defining use of GFS ensemble ! 2019-08-04 guo - moved ensemble object configuration into module gsi_fixture. +! 2022-03-15 K Apodaca - add CYGNSS and Spire Ocean wind speed observations +! 2023-03-15 K Apodaca - add GNSS-R Doppler Delay Map observations ! ! usage: ! input files: @@ -194,7 +196,7 @@ program gsi ! read_goesglm, read_goesndr, read_gps_ref, read_guess, read_ieeetovs, ! read_lidar, read_obs, read_ozone, read_pcp, read_prepbufr, read_radar, ! read_superwinds, read_wrf_mass_files, read_wrf_mass_guess, -! read_wrf_nmm_files, read_wrf_nmm_guess, rfdpar, rsearch, satthin, +! read_wrf_nmm_files, read_wrf_nmm_guess, read_gnssrspd, rfdpar, rsearch, satthin, ! setupdw, setupoz, setuppcp, setupps, setuppw, setupq, setuprad, ! setupref, setupbend, setuplight, setuprhsall, setuprw, setupspd, setupsst, ! setupt, setupw, simpin1, simpin1_init, smooth121, smoothrf, diff --git a/src/gsi/intgnssrspd.f90 b/src/gsi/intgnssrspd.f90 new file mode 100644 index 0000000000..1310ec7e14 --- /dev/null +++ b/src/gsi/intgnssrspd.f90 @@ -0,0 +1,239 @@ +module intgnssrspdmod +!$$$ module documentation block +! . . . . +! module: intgnssrspdmod module for intgnssrspd and its tangent linear i ntgnssrspd_tl +! (Based on intspd.f90) +! prgmmr: Apodaca, K. - Add CYGNSS TL and AD modules +! email: Karina.Apodaca@Spire.com +! +! abstract: module for intgnssrspd and its tangent linear intgnssrspd_tl +! +! program history log: +! 2023-04-21 k apodaca - initial version +! +! subroutines included: +! sub intgnssrspd_ +! +! variable definitions: +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +use m_obsNode, only: obsNode +use m_gnssrspdNode, only: gnssrspdNode +use m_gnssrspdNode, only: gnssrspdNode_typecast +use m_gnssrspdNode, only: gnssrspdNode_nextcast +use m_obsdiagNode, only: obsdiagNode_set +implicit none + +PRIVATE +PUBLIC intgnssrspd + +interface intgnssrspd; module procedure & + intgnssrspd_ +end interface + +contains + +subroutine intgnssrspd_(gnssrspdhead,rval,sval) +!$$$ subprogram documentation block +! . . . . +! subprogram: intgnssrspd apply nonlin qc obs operator for wind speed +! prgmmr: K. Apodaca org: Spire Global, Inc. date: 2023-03-15 +! email: Karina.Apodaca@Spire.com +! Based on intspd.f90 by prgmmr: derber org: np23 date: 1991-02-26 +! +! abstract: apply nonlinear observation operator and adjoint for winds +! +! program history log: +! 2023-04-21 k apodaca - initial version` +! +! input argument list: +! gnssrspdhead - obs type pointer to obs structure +! su - u increment in grid space +! sv - v increment in grid space +! ru +! rv +! +! output argument list: +! gnssrspdhead - obs type pointer to obs structure +! ru - u results from observation operator +! rv - v results from observation operator +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,i_kind + use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag + use qcmod, only: nlnqc_iter,varqc_iter + use constants, only: zero, half, one, tiny_r_kind,cg_term,r3600 + use gsi_4dvar, only: ltlint + use jfunc, only: jiter + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use gsi_4dvar, only: ladtest_obs + implicit none + +! Declare passed variables + class(obsNode), pointer, intent(in ) :: gnssrspdhead + type(gsi_bundle), intent(in ) :: sval + type(gsi_bundle), intent(inout) :: rval + +! Declare local variables + integer(i_kind) ier,istatus + integer(i_kind) j1,j2,j3,j4 + real(r_kind) w1,w2,w3,w4,term +! real(r_kind) penalty + real(r_kind) uanl,vanl,gnssrspdanl,gnssrspd,valv,valu + real(r_kind) uatl,vatl,gnssrspdatl,gnssrspdtra,grad + real(r_kind) cg_gnssrspd,p0,wnotgross,wgross,pg_gnssrspd + real(r_kind),pointer,dimension(:) :: su,sv + real(r_kind),pointer,dimension(:) :: ru,rv + type(gnssrspdNode), pointer :: gnssrspdptr + logical :: ltlint_tmp + +! If no gnssrspd data return + if(.not. associated(gnssrspdhead))return + +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier + if(ier/=0)return + + if( ladtest_obs) then + ltlint_tmp = ltlint + ltlint = .true. + end if + !gnssrspdptr => gnssrspdhead + gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead) + do while (associated(gnssrspdptr)) + + j1 = gnssrspdptr%ij(1) + j2 = gnssrspdptr%ij(2) + j3 = gnssrspdptr%ij(3) + j4 = gnssrspdptr%ij(4) + w1 = gnssrspdptr%wij(1) + w2 = gnssrspdptr%wij(2) + w3 = gnssrspdptr%wij(3) + w4 = gnssrspdptr%wij(4) + + + valu=zero + valv=zero + gnssrspdtra=sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges) + + if (ltlint) then + + if (gnssrspdtra>EPSILON(gnssrspdtra)) then +! Forward model + uatl=w1*su(j1)+w2*su(j2)+w3*su(j3)+w4*su(j4) + vatl=w1*sv(j1)+w2*sv(j2)+w3*sv(j3)+w4*sv(j4) + gnssrspdatl=gnssrspdptr%uges*uatl+gnssrspdptr%vges*vatl + gnssrspdatl=gnssrspdatl/gnssrspdtra + + if(luse_obsdiag)then + if (lsaveobsens) then + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspdatl + !-- gnssrspdptr%diags%obssen(jiter)=grad + call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad) + else + !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdatl + if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdatl) + endif + endif + + if (l_do_adjoint) then + if (.not. lsaveobsens) then + if( ladtest_obs ) then + grad=gnssrspdatl + else + gnssrspd=gnssrspdatl-gnssrspdptr%res+sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges) + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd + end if + endif + +! Adjoint + valu=grad*gnssrspdptr%uges/gnssrspdtra + valv=grad*gnssrspdptr%vges/gnssrspdtra + endif + else + if(luse_obsdiag)then + !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=zero + if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=zero) + !-- if (lsaveobsens) gnssrspdptr%diags%obssen(jiter)=zero + if (lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=zero) + end if + endif + + + else ! < ltlint > + +! Forward model + uanl=gnssrspdptr%uges+w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4) + vanl=gnssrspdptr%vges+w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4) + gnssrspdanl=sqrt(uanl*uanl+vanl*vanl) + !-- if (luse_obsdiag .and. gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdanl-gnssrspdtra + if (luse_obsdiag .and. gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdanl-gnssrspdtra) + + if (l_do_adjoint) then + valu=zero + valv=zero + gnssrspd=gnssrspdanl-gnssrspdptr%res + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd + +! Adjoint +! if(gnssrspdanl > tiny_r_kind*100._r_kind) then + if (gnssrspdanl>EPSILON(gnssrspdanl)) then + !-- if (luse_obsdiag .and. lsaveobsens) gnssrspdptr%diags%obssen(jiter)=grad + if (luse_obsdiag .and. lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad) + valu=uanl/gnssrspdanl + valv=vanl/gnssrspdanl + if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and. & + gnssrspdptr%b > tiny_r_kind) then + pg_gnssrspd=gnssrspdptr%pg*varqc_iter + cg_gnssrspd=cg_term/gnssrspdptr%b + wnotgross= one-pg_gnssrspd + wgross = pg_gnssrspd*cg_gnssrspd/wnotgross + p0 = wgross/(wgross+exp(-half*gnssrspdptr%err2*gnssrspd**2)) + term = (one-p0) + grad = grad*term + endif + end if + endif ! < l_do_adjoint > + + valu=valu*grad + valv=valv*grad + + endif ! < ltlint > + + + if (l_do_adjoint) then + ru(j1)=ru(j1)+w1*valu + ru(j2)=ru(j2)+w2*valu + ru(j3)=ru(j3)+w3*valu + ru(j4)=ru(j4)+w4*valu + rv(j1)=rv(j1)+w1*valv + rv(j2)=rv(j2)+w2*valv + rv(j3)=rv(j3)+w3*valv + rv(j4)=rv(j4)+w4*valv + + endif + + !gnssrspdptr => gnssrspdptr%llpoint + gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr) + + end do + if( ladtest_obs) ltlint = ltlint_tmp + return +end subroutine intgnssrspd_ + +end module intgnssrspdmod diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90 index a68355471b..4ce2952062 100644 --- a/src/gsi/intjo.f90 +++ b/src/gsi/intjo.f90 @@ -14,6 +14,7 @@ module intjomod ! 2016-08-29 J Guo - Separated calls to intozlay() and intozlev() ! 2018-08-10 guo - a new generic intjo() implementation replaced type ! specific intXYZ() calls with polymorphic %intjo(). +! 2022-03-15 K Apodaca - add CYGNSS and Spire ocean wind speed ! ! subroutines included: ! sub intjo_ @@ -31,7 +32,7 @@ module intjomod use gsi_obOperTypeManager, only: & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, iobOper_fed, & + iobOper_rw, iobOper_dbz, iobOper_fed, iobOper_gnssrspd, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & @@ -60,7 +61,7 @@ module intjomod integer(i_kind),parameter,dimension(obOper_count):: ix_obtype = (/ & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, iobOper_fed, & + iobOper_rw, iobOper_dbz, iobOper_fed, iobOper_gnssrspd, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & diff --git a/src/gsi/m_gnssrspdNode.F90 b/src/gsi/m_gnssrspdNode.F90 new file mode 100644 index 0000000000..1a8400440f --- /dev/null +++ b/src/gsi/m_gnssrspdNode.F90 @@ -0,0 +1,257 @@ +module m_gnssrspdNode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_gnssrspdNode +! prgmmr: K. Apodaca +! org: Spire Global, Inc. +! date: 2023-03-05 +! (based on m_*Node.F90 modules by: +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type gnssrspdNode (wind speed) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! 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_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsNode, only: obsNode + implicit none + private + + public:: gnssrspdNode + + type,extends(obsNode):: gnssrspdNode + !type(gnssrspd_ob_type),pointer :: llpoint => NULL() + type(obs_diag), pointer :: diags => NULL() + real(r_kind) :: res ! speed observation + real(r_kind) :: err2 ! speed error squared + real(r_kind) :: raterr2 ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b ! variational quality control parameter + real(r_kind) :: pg ! variational quality control parameter + real(r_kind) :: wij(4) ! horizontal interpolation weights + real(r_kind) :: uges ! guess u value + real(r_kind) :: vges ! guess v value + integer(i_kind) :: ij(4) ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !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 + real (r_kind) :: factw ! factor of 10m wind + contains + procedure,nopass:: mytype + procedure:: setHop => obsNode_setHop_ + procedure:: xread => obsNode_xread_ + procedure:: xwrite => obsNode_xwrite_ + procedure:: isvalid => obsNode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerRead => obsHeader_read_ + ! procedure, nopass:: headerWrite => obsHeader_write_ + ! procedure:: init => obsNode_init_ + ! procedure:: clean => obsNode_clean_ + end type gnssrspdNode + + public:: gnssrspdNode_typecast + public:: gnssrspdNode_nextcast + interface gnssrspdNode_typecast; module procedure typecast_ ; end interface + interface gnssrspdNode_nextcast; module procedure nextcast_ ; end interface + + public:: gnssrspdNode_appendto + interface gnssrspdNode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: MYNAME="m_gnssrspdNode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(aNode) result(ptr_) +!-- cast a class(obsNode) to a type(gnssrspdNode) + use m_obsNode, only: obsNode + implicit none + type(gnssrspdNode ),pointer:: ptr_ + class(obsNode),pointer,intent(in):: aNode + ptr_ => null() + if(.not.associated(aNode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(aNode) + type is(gnssrspdNode) + ptr_ => aNode + end select +return +end function typecast_ + +function nextcast_(aNode) result(ptr_) +!-- cast an obsNode_next(obsNode) to a type(gnssrspdNode) + use m_obsNode, only: obsNode,obsNode_next + implicit none + type(gnssrspdNode ),pointer:: ptr_ + class(obsNode),target ,intent(in):: aNode + + class(obsNode),pointer:: inode_ + inode_ => obsNode_next(aNode) + ptr_ => typecast_(inode_) +return +end function nextcast_ + +subroutine appendto_(aNode,oll) +!-- append aNode to linked-list oLL + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList,obsLList_appendNode + implicit none + type(gnssrspdNode),pointer,intent(in):: aNode + type(obsLList),intent(inout):: oLL + + class(obsNode),pointer:: inode_ + inode_ => aNode + call obsLList_appendNode(oLL,inode_) + inode_ => null() +end subroutine appendto_ + +! obsNode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[gnssrspdNode]" +end function mytype + +subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) + use m_obsdiagNode, only: obsdiagLookup_locate + implicit none + class(gnssrspdNode),intent(inout):: aNode + integer(i_kind),intent(in ):: iunit + integer(i_kind),intent( out):: istat + type(obs_diags),intent(in ):: diagLookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xread_' + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) aNode%res , & + aNode%err2 , & + aNode%raterr2, & + aNode%b , & + aNode%pg , & + aNode%uges , & + aNode%vges , & + aNode%factw , & + aNode%wij , & + aNode%ij + if (istat/=0) then + call perr(myname_,'read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + end if + + aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) + if(.not.associated(aNode%diags)) then + call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) + call perr(myname_,' %iob =',aNode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) +return +end subroutine obsNode_xread_ + +subroutine obsNode_xwrite_(aNode,junit,jstat) + implicit none + class(gnssrspdNode),intent(in):: aNode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xwrite_' +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) aNode%res , & + aNode%err2 , & + aNode%raterr2, & + aNode%b , & + aNode%pg , & + aNode%uges , & + aNode%vges , & + aNode%factw , & + aNode%wij , & + aNode%ij + if (jstat/=0) then + call perr(myname_,'write(%(res,err2,...)), iostat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) +return +end subroutine obsNode_xwrite_ + +subroutine obsNode_setHop_(aNode) + use m_cvgridLookup, only: cvgridLookup_getiw + implicit none + class(gnssrspdNode),intent(inout):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' +_ENTRY_(myname_) + call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,aNode%wij) + aNode%wij(1:4) = aNode%wij(1:4)*aNode%factw +_EXIT_(myname_) +return +end subroutine obsNode_setHop_ + +function obsNode_isvalid_(aNode) result(isvalid_) + implicit none + logical:: isvalid_ + class(gnssrspdNode),intent(in):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' +_ENTRY_(myname_) + isvalid_=associated(aNode%diags) +_EXIT_(myname_) +return +end function obsNode_isvalid_ + +pure subroutine gettlddp_(aNode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(gnssrspdNode), intent(in):: aNode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 +return +end subroutine gettlddp_ + +end module m_gnssrspdNode diff --git a/src/gsi/m_obsNodeTypeManager.F90 b/src/gsi/m_obsNodeTypeManager.F90 index 43b42e4bf2..9da356c808 100644 --- a/src/gsi/m_obsNodeTypeManager.F90 +++ b/src/gsi/m_obsNodeTypeManager.F90 @@ -15,6 +15,8 @@ module m_obsNodeTypeManager ! 2018-01-23 k apodaca - add a new observation type i.e. lightning (light) ! suitable for the GOES/GLM instrument ! +! 2024-04-07 k apodaca - add gnssrspdnode +! ! input argument list: see Fortran 90 style document below ! ! output argument list: see Fortran 90 style document below @@ -54,6 +56,8 @@ module m_obsNodeTypeManager use m_wspd10mNode, only: wspd10mNode use m_uwnd10mNode, only: uwnd10mNode use m_vwnd10mNode, only: vwnd10mNode + + use m_gnssrspdNode, only: gnssrspdNode use m_td2mNode , only: td2mNode use m_mxtmNode , only: mxtmNode @@ -122,10 +126,10 @@ module m_obsNodeTypeManager public:: iobsNode_cldch public:: iobsNode_swcp public:: iobsNode_lwcp - public:: iobsNode_light public:: iobsNode_dbz public:: iobsNode_fed + public:: iobsNode_gnssrspd public :: obsNode_typeMold public :: obsNode_typeIndex @@ -166,6 +170,8 @@ module m_obsNodeTypeManager type(wspd10mNode), target, save:: wspd10m_mold type(uwnd10mNode), target, save:: uwnd10m_mold type(vwnd10mNode), target, save:: vwnd10m_mold + + type(gnssrspdNode), target, save:: gnssrspd_mold type( td2mNode), target, save:: td2m_mold type( mxtmNode), target, save:: mxtm_mold @@ -214,6 +220,7 @@ module m_obsNodeTypeManager enumerator:: iobsNode_w enumerator:: iobsNode_q enumerator:: iobsNode_spd + enumerator:: iobsNode_gnssrspd enumerator:: iobsNode_rw enumerator:: iobsNode_dw enumerator:: iobsNode_sst @@ -301,6 +308,8 @@ function vname2index_(vname) result(index_) "[uwnd10mnode]"); index_ = iobsNode_uwnd10m case("vwnd10m", & "[vwnd10mnode]"); index_ = iobsNode_vwnd10m + case("gnssrspd", & + "[gnssrspdnode]"); index_ = iobsNode_gnssrspd case("td2m" , "[td2mnode]"); index_ = iobsNode_td2m case("mxtm" , "[mxtmnode]"); index_ = iobsNode_mxtm @@ -365,6 +374,8 @@ function vmold2index_select_(mold) result(index_) type is(wspd10mNode); index_ = iobsNode_wspd10m type is(uwnd10mNode); index_ = iobsNode_uwnd10m type is(vwnd10mNode); index_ = iobsNode_vwnd10m + + type is(gnssrspdNode); index_ = iobsNode_gnssrspd type is( td2mNode); index_ = iobsNode_td2m type is( mxtmNode); index_ = iobsNode_mxtm @@ -423,6 +434,8 @@ function index2vmold_(i_obType) result(obsmold_) case(iobsNode_wspd10m); obsmold_ => wspd10m_mold case(iobsNode_uwnd10m); obsmold_ => uwnd10m_mold case(iobsNode_vwnd10m); obsmold_ => vwnd10m_mold + + case(iobsNode_gnssrspd ); obsmold_ => gnssrspd_mold case(iobsNode_td2m ); obsmold_ => td2m_mold case(iobsNode_mxtm ); obsmold_ => mxtm_mold diff --git a/src/gsi/m_rhs.F90 b/src/gsi/m_rhs.F90 index aea417fe27..723b0a6918 100644 --- a/src/gsi/m_rhs.F90 +++ b/src/gsi/m_rhs.F90 @@ -17,6 +17,7 @@ module m_rhs ! through an enum block. ! - removed external dimension argument aworkdim2 of ! rhs_alloc(). +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed index (i_gnssrspd) ! ! input argument list: see Fortran 90 style document below ! @@ -66,6 +67,7 @@ module m_rhs public:: i_vis public:: i_pblh public:: i_wspd10m + public:: i_gnssrspd public:: i_td2m public:: i_mxtm public:: i_mitm @@ -133,6 +135,7 @@ module m_rhs enumerator:: i_vis enumerator:: i_pblh enumerator:: i_wspd10m + enumerator:: i_gnssrspd enumerator:: i_td2m enumerator:: i_mxtm enumerator:: i_mitm @@ -207,7 +210,6 @@ subroutine rhs_alloc() allocate(rhs_stats_oz(9,jpch_oz)) allocate(rhs_toss_gps(max(1,nprof_gps))) - rhs_awork =zero rhs_bwork =zero rhs_aivals =zero diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 1c45c62bc8..19f6a072fb 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -160,6 +160,7 @@ module obsmod ! 2021-11-16 Zhao - add option l_obsprvdiag (if true) to trigger the output of ! observation provider and sub-provider information into ! obsdiags files (used for AutoObsQC) +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed observations (CYGNSS, Spire) ! 2023-07-10 Y. Wang, D. Dowell - add variables for flash extent density ! 2023-10-10 H. Wang (GSL) - add variables for flash extent density EnVar DA ! @@ -463,8 +464,8 @@ module obsmod public :: iout_lag,iout_uv,iout_gps,iout_ps,iout_light,mype_light public :: mype_gust,mype_vis,mype_pblh,iout_gust,iout_vis,iout_pblh public :: mype_tcamt,mype_lcbas,iout_tcamt,iout_lcbas - public :: mype_wspd10m,mype_td2m,iout_wspd10m,iout_td2m - public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m + public :: mype_wspd10m,mype_gnssrspd,mype_td2m,iout_wspd10m,iout_gnssrspd,iout_td2m + public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m public :: mype_mxtm,mype_mitm,iout_mxtm,iout_mitm public :: mype_pmsl,mype_howv,iout_pmsl,iout_howv public :: mype_swcp,mype_lwcp,iout_swcp,iout_lwcp @@ -600,12 +601,12 @@ module obsmod integer(i_kind) iout_dw,iout_gps,iout_sst,iout_tcp,iout_lag integer(i_kind) iout_co,iout_gust,iout_vis,iout_pblh,iout_tcamt,iout_lcbas integer(i_kind) iout_cldch - integer(i_kind) iout_wspd10m,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv + integer(i_kind) iout_wspd10m,iout_gnssrspd,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv integer(i_kind) iout_uwnd10m,iout_vwnd10m,iout_fed integer(i_kind) mype_t,mype_q,mype_uv,mype_ps,mype_pw, & mype_rw,mype_dw,mype_gps,mype_sst, & mype_tcp,mype_lag,mype_co,mype_gust,mype_vis,mype_pblh, & - mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,& + mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,& mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz, mype_fed integer(i_kind) mype_cldch integer(i_kind) iout_swcp, iout_lwcp @@ -887,6 +888,7 @@ subroutine init_obsmod_dflts iout_light=237 ! lightning iout_dbz=238 ! radar reflectivity iout_fed=239 ! flash extent density + iout_gnssrspd=240 ! GNSS-R wind speed mype_ps = npe-1 ! surface pressure mype_t = max(0,npe-2) ! temperature @@ -922,6 +924,7 @@ subroutine init_obsmod_dflts mype_light=max(0,npe-32)! GOES/GLM lightning mype_dbz=max(0,npe-33) ! radar reflectivity mype_fed= max(0,npe-34) ! flash extent density + mype_gnssrspd= max(0,npe-35) ! surface GNSS-R speed ! Initialize arrays used in namelist obs_input diff --git a/src/gsi/read_gnssrspd.f90 b/src/gsi/read_gnssrspd.f90 new file mode 100644 index 0000000000..859345a9cc --- /dev/null +++ b/src/gsi/read_gnssrspd.f90 @@ -0,0 +1,463 @@ +subroutine read_gnssrspd(nread,ndata,nodata,infile,obstype,lunout,twind,sis,& + nobs) + +!$$$ subprogram documentation block +! . . . . +! subprogram: read_gnssrspd read obs from gnssrspd bufr file +! prgmmr: kapodaca org: Spire Global, Inc. date: 2022-03-12 +! Largely based on other read_* routines. For the !Winds --- surface wind speed section +! ~L438, we used the exixting read_fl_hdob.f90 routine +! +! abstract: This routine reads GNSSRSPD L2 wind speed observations +! +! When running the gsi in regional mode, the code only +! retains those observations that fall within the regional +! domain + +! program history log: +! 2022-03-12 k apodaca- initial coding +! +! input argument list: +! infile - unit from which to read BUFR data +! obstype - observation type to process +! lunout - unit to which to write data for further processing +! twind - input group time window (hours) +! sis - satellite/instrument/sensor indicator +! +! output argument list: +! nread - number of type "obstype" observations read +! nodata - number of individual "obstype" observations read +! ndata - number of type "obstype" observations retained for further processing +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: +! +!$$$ + use kinds, only: r_single,r_kind,r_double,i_kind + use constants, only: zero,one_tenth,one,two,ten,deg2rad,t0c,half,& + three,four,rad2deg,r0_01,& + r60inv,r10,r100,r2000,hvap + use gridmod, only: diagnostic_reg,regional,nlon,nlat,& + tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& + rlats,rlons,twodvar_regional + use convinfo, only: nconvtype, & + icuse,ictype,ioctype + use obsmod, only: ran01dom + use obsmod, only: iadate,bmiss,offtime_data + use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen + use qcmod, only: errormod + use convthin, only: make3grids,del3grids + use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error + use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 + use mpimod, only: npe + + implicit none + +! Declare passed variables + character(len=*), intent(in ) :: infile,obstype + character(len=20),intent(in ) :: sis + integer(i_kind) , intent(in ) :: lunout + integer(i_kind) , dimension(npe), intent(inout) :: nobs + integer(i_kind) , intent(inout) :: nread,ndata,nodata + real(r_kind) , intent(in ) :: twind + +! Declare local variables +! Logical variables + logical :: outside + logical :: inflate_error + logical :: lspdob + +! Character variables + character(40) :: timestr,locstr,wndstr,oestr + character( 8) :: subset + character( 8) :: c_prvstg,c_sprvstg + character( 8) :: c_station_id + character(10) date +! Integer variables + integer(i_kind), parameter :: mxib = 31 + + integer(i_kind) :: i,k + integer(i_kind) :: ihh,idd,idate,im,iy + integer(i_kind) :: lunin + integer(i_kind) :: ireadmg,ireadsb + integer(i_kind) :: ilat,ilon + integer(i_kind) :: nlv + integer(i_kind) :: nreal,nchanl + integer(i_kind) :: idomsfc,isflg + integer(i_kind) :: iout + integer(i_kind) :: nc,ncsave + integer(i_kind) :: ntmatch,ntb + integer(i_kind) :: nmsg + integer(i_kind) :: maxobs + integer(i_kind) :: itype,iecol + integer(i_kind) :: qcm,lim_qm + integer(i_kind) :: wspd_qm + integer(i_kind) :: ntest,nvtest + integer(i_kind) :: igood + integer(i_kind) :: iuse + + integer(i_kind) :: idate5(5) + integer(i_kind) :: minobs,minan + + integer(i_kind), allocatable,dimension(:) :: isort + +! Real variables + real(r_kind), parameter :: r0_001 = 0.001_r_kind + real(r_kind), parameter :: r1_2 = 1.2_r_kind + real(r_kind), parameter :: r3_0 = 3.0_r_kind + real(r_kind), parameter :: r0_7 = 0.7_r_kind + real(r_kind), parameter :: r6 = 6.0_r_kind + real(r_kind), parameter :: r50 = 50.0_r_kind + real(r_kind), parameter :: r1200 = 1200.0_r_kind + real(r_kind), parameter :: emerr = 0.2_r_kind ! RH + real(r_kind), parameter :: missing = 1.0e+11_r_kind + real(r_kind), parameter :: r180 = 180.0_r_kind + real(r_kind), parameter :: r360 = 360.0_r_kind + + real(r_kind) :: toff,t4dv + real(r_kind) :: usage + real(r_kind) :: woe,obserr + real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth + real(r_kind) :: dlat_earth_deg,dlon_earth_deg + real(r_kind) :: cdist,disterr,disterrmax,rlon00,rlat00 + real(r_kind) :: vdisterrmax + real(r_kind) :: spdob + real(r_kind) :: gob + real(r_kind) :: dlnpsob + real(r_kind) :: tdiff + real(r_kind) :: tsavg,ff10,sfcr,zz + real(r_kind) :: errmin + real(r_kind) :: log100 + + real(r_kind) :: obstime(6,1) + real(r_kind) :: obsloc(2,1) + real(r_kind) :: gnssrw(2,1) + + real(r_double) :: rstation_id + real(r_double) :: r_prvstg(1,1),r_sprvstg(1,1) + + real(r_kind), allocatable,dimension(:,:) :: cdata_all,cdata_out + +! Equivalence to handle character names + equivalence(r_prvstg(1,1),c_prvstg) + equivalence(r_sprvstg(1,1),c_sprvstg) + equivalence(rstation_id,c_station_id) + +! Data + !data timestr / 'YEAR MNTH DAYS HOUR MINU SECO' / + data timestr / 'DHR RPT' / + !data locstr / 'CLAT CLON' / + data locstr / 'XOB YOB' / + !data wndstr / 'WSPD' / !GNSSRSPD Wind speed + data wndstr / 'SOB' / !GNSSRSPD Wind speed + data oestr / 'WSU' / !GNSSRSPD Wind speed uncertainty/error + data lunin / 13 / + +!------------------------------------------------------------------------------------------------ + + write(6,*)'READ_GNSSRSPD: begin to read gnssrspd satellite data ...' + +! Initialize parameters + +! Set common variables + lspdob = obstype == 'gnssrspd' + + nreal = 0 + iecol = 0 + log100=log(100._r_kind) + + + lim_qm = 4 + iecol=0 + if (lspdob) then + nreal = 23 + iecol = 4 + errmin = one + else + write(6,*) ' illegal obs type in read_gnssrspd ' + call stop2(94) + end if + + inflate_error = .true. + +! Check if the obs type specified in the convinfo is in the fl hdob bufr file +! If found, get the index (nc) from the convinfo for the specified type +ntmatch = 0 +ncsave = 0 +do nc = 1, nconvtype + if (trim(ioctype(nc)) == trim(obstype)) then + if (trim(ioctype(nc)) == 'gnssrspd' .and. ictype(nc) == 600 ) then + ntmatch = ntmatch + 1 + ncsave = nc + itype = ictype(nc) + end if + end if +end do +if (ntmatch == 0) then ! Return if not specified in convinfo + write(6,*) ' READ_GNSSRSPD: No matching obstype found in obsinfo ', obstype + return +else + nc = ncsave + write(6,*) ' READ_GNSSRSPD: Processing GNSSRSPD data : ', ntmatch, nc, ioctype(nc), ictype(nc), itype +end if + + +!------------------------------------------------------------------------------------------------ + +! Go through the bufr file to find out how mant subsets to process + nmsg = 0 + maxobs = 0 + call closbf(lunin) + open(lunin,file=trim(infile),form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + + loop_msg1: do while(ireadmg(lunin,subset,idate) >= 0) + if(nmsg == 0) call time_4dvar(idate,toff) ! time offset (hour) + + nmsg = nmsg+1 + loop_readsb1: do while(ireadsb(lunin) == 0) + maxobs = maxobs+1 + end do loop_readsb1 + end do loop_msg1 + call closbf(lunin) + write(6,*) 'READ_GNSSRSPD: total number of data found in the bufr file ',maxobs,obstype + write(6,*) 'READ_GNSSRSPD: time offset is ',toff,' hours' + +!--------------------------------------------------------------------------------------------------- + +! Allocate array to hold data + allocate(cdata_all(nreal,maxobs)) + allocate(isort(maxobs)) + +! Initialize + cdata_all = zero + isort = 0 + nread = 0 + nchanl = 0 + ntest = 0 + nvtest = 0 + ilon = 2 + ilat = 3 + +! Open bufr file again for reading + call closbf(lunin) + open(lunin,file=trim(infile),form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + ntb = 0 + igood = 0 +! Loop through BUFR file + loop_msg2: do while(ireadmg(lunin,subset,idate) >= 0) + loop_readsb2: do while(ireadsb(lunin) == 0) + + ntb = ntb+1 + + c_station_id = subset + +! QC mark 9: will be monitored but not assimilated +! QC mark 4: reject - will not be monitored nor assimilated +! QC mark 3: suspect +! QC mark 2: neutral or not checked +! QC mark 1: good +! QC mark 0: keep - will be always assimilated + qcm = 0 + wspd_qm = 0 + + +! Read observation time + call ufbint(lunin,obstime,2,1,nlv,timestr) + +! If date in gnssrspd file does not agree with analysis date, +! print message and stop program execution. + write(date,'( i10)') idate + read (date,'(i4,3i2)') iy,im,idd,ihh + if(offtime_data) then + +! in time correction for observations to account for analysis + idate5(1)=iy + idate5(2)=im + idate5(3)=idd + idate5(4)=ihh + idate5(5)=0 + call w3fs21(idate5,minobs) ! obs ref time in minutes relative to historic date + idate5(1)=iadate(1) + idate5(2)=iadate(2) + idate5(3)=iadate(3) + idate5(4)=iadate(4) + idate5(5)=0 + call w3fs21(idate5,minan) ! analysis ref time in minutes relative to historic date +! Add obs reference time, then subtract analysis time to get obs time relative to analysis + + tdiff=float(minobs-minan)*r60inv + + else + tdiff=zero + end if + + t4dv = toff+obstime(1,1) + + if (l4dvar.or.l4densvar) then + if (t4dv < zero .OR. t4dv > winlen) cycle loop_readsb2 + else + if (abs(tdiff)>twind) cycle loop_readsb2 + endif + nread = nread+1 + + usage = zero ! will be considered for assimilation + ! subject to further QC in setupt subroutine + iuse = icuse(nc) ! assimilation flag + if (iuse <=0) usage = r100 ! will be monitored but not assimilated + +! Read observation location (lat/lon degree) + call ufbint(lunin,obsloc,2,1,nlv,locstr) + + if (obsloc(1,1) == missing .or. abs(obsloc(1,1)) < -180.0_r_kind .or. & + obsloc(1,1) == missing .or. abs(obsloc(1,1)) > 180.0_r_kind .or. & + obsloc(2,1) == missing .or. abs(obsloc(2,1)) < -90.0_r_kind .or. & + obsloc(2,1) == missing .or. abs(obsloc(2,1)) > 90.0_r_kind ) then + write(6,*) 'READ_GNSSRSPD: bad lon/lat values: ', obsloc(1,1),obsloc(2,1) + cycle loop_readsb2 + endif +! GNSSRSPD BUFR longitudes are in the +/- 180 range, need to convert to 0 to 360 +! deg + if (obsloc(1,1) < 0.0_r_kind) obsloc(1,1) = obsloc(1,1) + 360.0_r_kind + !if (obsloc(1,1) >= 0.0_r_kind .or. obsloc(1,1) <= 180.0_r_kind) obsloc(1,1) = obsloc(1,1) + 180.0_r_kind + + dlon_earth_deg = obsloc(1,1) + dlat_earth_deg = obsloc(2,1) + dlon_earth = obsloc(1,1)*deg2rad ! degree to radian + dlat_earth = obsloc(2,1)*deg2rad ! degree to radian + +! Convert obs lat/lon to rotated coordinate and check +! if the obs is outside of the regional domain + if (regional) then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if (diagnostic_reg) then + call txy2ll(dlon,dlat,rlon00,rlat00) + ntest = ntest+1 + cdist = sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & + (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) + cdist = max(-one,min(cdist,one)) + disterr = acos(cdist)*rad2deg + disterrmax = max(disterrmax,disterr) + end if + if(outside) cycle loop_readsb2 + else + dlon = dlon_earth + dlat = dlat_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + + +! Read surface wind speed [m/s] from GNSSRSPD + if (lspdob) then + !usage=r100 + !usage = zero ! will be considered for assimilation + + ! Get Wind Speed observations from bufr file + call ufbint(lunin,gnssrw,1,1,nlv,wndstr) + spdob = gnssrw(1,1) ! surface wind speed + + ! Don't permit observations with ws <= 1 m/s + if (spdob <= one) cycle loop_readsb2 + if (spdob >= missing) cycle loop_readsb2 + + ! Get observation error from bufr file + call ufbint(lunin,gnssrw,1,1,nlv,oestr) + obserr = max(gnssrw(1,1),1.5_r_kind) ! surface wind speed error + endif + + + if ( .not. twodvar_regional) then + call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg) + endif + + ! Get information from surface file necessary for conventional data + call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) + ! Process data passed quality control + igood = igood + 1 + ndata = ndata + 1 + nodata = nodata + 1 + iout = ndata + +! Read extrapolated surface pressure [pa] and convert to [cb] + dlnpsob = log100 ! default (1000mb) + +!------------------------------------------------------------------------------------------------- + + ! Winds --- surface wind speed + if (lspdob) then + woe = obserr + !if (inflate_error) woe = woe*r3_0 + !if (inflate_error) woe = woe*r1_2 + !if (qcm > lim_qm ) woe = woe*1.0e6_r_kind + cdata_all( 1,iout)=woe ! wind error + cdata_all( 2,iout)=dlon ! grid relative longitude + cdata_all( 3,iout)=dlat ! grid relative latitude + cdata_all( 4,iout)=dlnpsob ! ln(surface pressure in cb) + cdata_all( 5,iout)=spdob*sqrt(two)*half ! u obs + cdata_all( 6,iout)=spdob*sqrt(two)*half ! v obs + cdata_all( 7,iout)=rstation_id ! station id + cdata_all( 8,iout)=t4dv ! time + cdata_all( 9,iout)=nc ! type + cdata_all(10,iout)=r10 ! elevation of observation ! 10-m wind + cdata_all(11,iout)=qcm ! quality mark + cdata_all(12,iout)=obserr ! original obs error + cdata_all(13,iout)=usage ! usage parameter + cdata_all(14,iout)=idomsfc ! dominate surface type + cdata_all(15,iout)=tsavg ! skin temperature + cdata_all(16,iout)=ff10 ! 10 meter wind factor + cdata_all(17,iout)=sfcr ! surface roughness + cdata_all(18,iout)=dlon_earth_deg ! earth relative longitude (degree) + cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degree) + cdata_all(20,iout)=gob ! station elevation (m) + cdata_all(21,iout)=zz ! terrain height at ob location + cdata_all(22,iout)=r_prvstg(1,1) ! provider name + cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name + endif + + end do loop_readsb2 + end do loop_msg2 + +! Close unit to bufr file + call closbf(lunin) + +! Write header record and data to output file for further processing + write(6,*) "READ_GNSSRSPD: nreal=",nreal," ndata=",ndata + allocate(cdata_out(nreal,ndata)) + do i=1,ndata + do k=1,nreal + cdata_out(k,i)=cdata_all(k,i) + end do + end do + deallocate(cdata_all) + + call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon + write(lunout) cdata_out + deallocate(cdata_out) +900 continue + if(diagnostic_reg .and. ntest>0) write(6,*)'READ_GNSSRSPD: ',& + 'ntest, disterrmax=', ntest,disterrmax + if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GNSSRSPD: ',& + 'nvtest,vdisterrmax=',ntest,vdisterrmax + + if (ndata == 0) then + call closbf(lunin) + write(6,*)'READ_GNSSRSPD: no data to process' + endif + write(6,*)'READ_GNSSRSPD: nreal=',nreal + write(6,*)'READ_GNSSRSPD: ntb,nread,ndata,nodata=',ntb,nread,ndata,nodata + + call closbf(lunin) + close(lunin) + +! End of routine + return + +end subroutine read_gnssrspd + diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index aa0f11b4e3..6f0e32a4b6 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -141,6 +141,9 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) ! 2019-08-21 H. Shao - add METOPC-C, COSMIC-2 and PAZ to the GPS check list ! 2020-05-21 H. Shao - add commercial GNSSRO (Spire, PlanetIQ, GeoOptics) and other existing missions to the check list ! 2021-02-20 X.Li - add viirs-m and get_hsst +! 2022-03-12 K. Apodaca - Enable CYGNSS ocean wind speed reading capability +! 2023-03-12 K. Apodaca - Enable GNSS-R L2 ocean wind speed reading (CYGNSS, Spire) +! 2023-03-12 K. Apodaca - Enable GNSS-R DDM reading capability ! ! input argument list: ! lexist - file status @@ -487,6 +490,16 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) endif nread = nread + 1 end do loop_hdob + else if(trim(filename) == 'gnssrwndbufr')then + lexist = .false. + gnssrwndloop: do while(ireadmg(lnbufr,subset,idate2) >= 0) + if (trim(dtype)=='gnssrspd') then + lexist = .true. + exit gnssrwndloop + endif + nread = nread + 1 + end do gnssrwndloop + else if(trim(dtype) == 'pm2_5')then if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then lexist=.true. @@ -901,7 +914,7 @@ subroutine read_obs(ndata,mype) if (obstype == 't' .or. obstype == 'uv' .or. & obstype == 'q' .or. obstype == 'ps' .or. & obstype == 'pw' .or. obstype == 'spd'.or. & - obstype == 'sst'.or. & + obstype == 'sst'.or. obstype == 'gnssrspd'.or. & obstype == 'tcp'.or. obstype == "lag".or. & obstype == 'dw' .or. obstype == 'rw' .or. & obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. & @@ -1449,8 +1462,13 @@ subroutine read_obs(ndata,mype) 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 == 'gnssrspd' .and. index(infile,'gnssrwndbufr') /=0 ) then + call read_gnssrspd(nread,npuse,nouse,infile,obstype,lunout,twind,sis, & + nobs_sub1(1,i)) + string='READ_GNSSRSPD' + else if(obstype == 'howv') then if ( index(infile,'satmar') /=0) then diff --git a/src/gsi/setupgnssrspd.f90 b/src/gsi/setupgnssrspd.f90 new file mode 100644 index 0000000000..e90b94a1bd --- /dev/null +++ b/src/gsi/setupgnssrspd.f90 @@ -0,0 +1,940 @@ +module gnssrspd_setup + implicit none + private + public:: setup + interface setup; module procedure setupgnssrspd; end interface + +contains +subroutine setupgnssrspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsave) +!$$$ subprogram documentation block +! . . . . +! subprogram: setupgnssrspd compute rhs of oi for wind speed obs +! prgmmr: Karina Apodaca org: Spire Global, Inc. date: 2023-03-15 +! Based on: setup10mspd.f90 by: +! prgmmr: parrish org: np22 date: 1990-10-06 +! +! abstract: For wind speed observations, this routine +! a) reads obs assigned to given mpi task (geographic region), +! b) simulates obs from guess, +! c) apply some quality control to obs, +! d) load weight and innovation arrays used in minimization +! e) collects statistics for runtime diagnostic output +! f) writes additional diagnostic information to output file +! +! program history log: +! 2023-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed +! +! input argument list: +! lunin - unit from which to read observations +! mype - mpi task id +! nele - number of data elements per observation +! nobs - number of observations +! +! 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 mpeu_util, only: die,perr,getindex + use kinds, only: r_kind,r_single,r_double,i_kind + use m_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use m_obsdiagNode, only: obsdiagLList_nextNode + use m_obsdiagNode, only: obsdiagNode_set + use m_obsdiagNode, only: obsdiagNode_get + use m_obsdiagNode, only: obsdiagNode_assert + + use obsmod, only: rmiss_single,& + lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,& + lobsdiag_forenkf + use obsmod, only: netcdf_diag, binary_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 m_obsNode, only: obsNode + use m_gnssrspdNode, only: gnssrspdNode + use m_gnssrspdNode, only: gnssrspdNode_appendto + use m_obsLList, only: obsLList + use obsmod, only: luse_obsdiag + use gsi_4dvar, only: nobs_bins,hr_obsbin + use guess_grids, only: nfldsig,hrdifsig,ges_lnprsl, & + comp_fact10,sfcmod_gfs,sfcmod_mm5 + use guess_grids, only: geop_hgtl + use gridmod, only: nsig,get_ij,twodvar_regional + use qcmod, only: npres_print,ptop,pbot + use constants, only: one,grav,rd,zero,four,tiny_r_kind, & + half,two,cg_term,huge_single,r1000,wgtlim,ten + use jfunc, only: jiter,last,miter,jiterstart + use state_vectors, only: svars3d, levels + use qcmod, only: dfact,dfact1 + use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype + use convinfo, only: icsubtype + use gsi_bundlemod, only : gsi_bundlegetpointer + use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle + use m_dtime, only: dtime_setup, dtime_check + use sparsearr, only: sparr2, new, size, writearray, fullarray + + ! The following variables are the coefficients that describe the + ! linear regression fits that are used to define the dynamic + ! observation error (DOE) specifications for all reconnissance + ! observations collected within hurricanes/tropical cyclones; these + ! apply only to the regional forecast models (e.g., HWRF); Henry + ! R. Winterbottom (henry.winterbottom@noaa.gov). + + + implicit none + +! Declare passed variables + type(obsLList ),target,dimension(:),intent(in):: obsLL + type(obs_diags),target,dimension(:),intent(in):: odiagLL + + logical ,intent(in ) :: conv_diagsave + integer(i_kind) ,intent(in ) :: lunin,mype,nele,nobs + real(r_kind),dimension(100+7*nsig) ,intent(inout) :: awork + real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout) :: bwork + integer(i_kind) ,intent(in ) :: is ! ndat index + +! Declare local variables + character(len=*),parameter:: myname='setupgnssrspd' + +! Declare external calls for code analysis + external:: tintrp2a1,tintrp2a11 + external:: tintrp31 + external:: grdcrd1 + external:: stop2 + +! Declare local variables + + real(r_double) rstation_id + + real(r_kind) uob,vob,gnssrspdges,gnssrspdob,gnssrspdob0,goverrd,ratio_errors + real(r_kind) presw,factw,dpres,ugesin,vgesin,sfcr,skint + real(r_kind) scale + real(r_kind) val2,ressw,ress,error,ddiff,dx10,rhgh,prsfc,r0_001 + real(r_kind) sfcchk,prsln2,rwgt,tfact + real(r_kind) thirty,rsig,ratio,residual,obserrlm,obserror + real(r_kind) val,valqc,psges,drpx,dlat,dlon,dtime,rlow + real(r_kind) cg_gnssrspd,wgross,wnotgross,wgt,arg,exp_arg,term,rat_err2 + real(r_kind) errinv_input,errinv_adjst,errinv_final + real(r_kind) err_input,err_adjst,err_final + real(r_kind),dimension(nele,nobs):: data + real(r_kind),dimension(nobs):: dup + real(r_kind),dimension(nsig)::prsltmp,tges + real(r_single),allocatable,dimension(:,:)::rdiagbuf + + integer(i_kind) mm1,ibin,ioff,ioff0 + integer(i_kind) ii,jj,i,nchar,nreal,k,j,l,nty,nn,ikxx + integer(i_kind) ier,ilon,ilat,ipres,iuob,ivob,id,itime,ikx + integer(i_kind) ihgt,iqc,ier2,iuse,ilate,ilone,istnelv,izz,iprvd,isprvd + integer(i_kind) idomsfc,iskint,iff10,isfcr,isli + + type(sparr2) :: dhx_dx + integer(i_kind) :: iz, u_ind, v_ind, nnz, nind + real(r_kind) :: delz + + logical,dimension(nobs):: luse,muse + integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID + logical proceed + + character(8) station_id + character(8),allocatable,dimension(:):: cdiagbuf + character(8),allocatable,dimension(:):: cprvstg,csprvstg + character(8) c_prvstg,c_sprvstg + real(r_double) r_prvstg,r_sprvstg + + logical:: in_curbin, in_anybin, save_jacobian + type(gnssrspdNode),pointer:: my_head + type(obs_diag),pointer:: my_diag + type(obs_diags),pointer:: my_diagLL + + logical z_height + real(r_kind) zsges,dstn + real(r_kind),dimension(nsig):: zges + real(r_kind) dz,zob,z1,z2,p1,p2,dz21,dlnp21,pobl + integer(i_kind) k1,k2 + + equivalence(rstation_id,station_id) + equivalence(r_prvstg,c_prvstg) + equivalence(r_sprvstg,c_sprvstg) + + real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps + real(r_kind),allocatable,dimension(:,:,: ) :: ges_z + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv + + type(obsLList),pointer,dimension(:):: gnssrspdhead + gnssrspdhead => obsLL(:) + +!****************************************************************************** +! Read and reformat observations in work arrays. + read(lunin)data,luse,ioid + +! index information for data array (see reading routine) + ier=1 ! index of obs error + ilon=2 ! index of grid relative obs location (x) + ilat=3 ! index of grid relative obs location (y) + ipres=4 ! index of pressure + iuob=5 ! index of u observation + ivob=6 ! index of v observation + id=7 ! index of station id + itime=8 ! index of observation time in data array + ikxx=9 ! index of ob type + ihgt=10 ! index of observation elevation + iqc=11 ! index of quality mark + ier2=12 ! index of original-original obs error ratio + iuse=13 ! index of use parameter + idomsfc=14 ! index of dominant surface type + iskint=15 ! index of surface skin temperature + iff10=16 ! index of 10 meter wind factor + isfcr=17 ! index of surface roughness + ilone=18 ! index of longitude (degrees) + ilate=19 ! index of latitude (degrees) + istnelv=20 ! index of station elevation (m) + izz=21 ! index of surface height + iprvd=22 ! index of observation provider + isprvd=23 ! index of observation subprovider + + mm1=mype+1 + scale=one + rsig=nsig + thirty = 30.0_r_kind + r0_001=0.001_r_kind + goverrd=grav/rd + + save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf + +! Check to see if required guess fields are available + call check_vars_(proceed) + if(.not.proceed) return ! not all vars available, simply return + +! If require guess vars available, extract from bundle ... + call init_vars_ + +! If requested, save select data for output to diagnostic file + if(conv_diagsave)then + ii=0 + nchar=1 + ioff0=21 + nreal=ioff0 + if (lobsdiagsave) nreal=nreal+4*miter+1 + if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif + if (save_jacobian) then + nnz = 4 ! number of non-zero elements in dH(x)/dx profile + nind = 2 + call new(dhx_dx, nnz, nind) + nreal = nreal + size(dhx_dx) + endif + allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs)) + if(netcdf_diag) call init_netcdf_diag_ + end if + + !!awork(:) = zero + + do i=1,nobs + muse(i)=nint(data(iuse,i)) <= jiter + end do + + dup=one + do k=1,nobs + do l=k+1,nobs + if(data(ilat,k) == data(ilat,l) .and. & + data(ilon,k) == data(ilon,l) .and. & + data(ipres,k)== data(ipres,l) .and. & + data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & + muse(l) .and. muse(k))then + + tfact=min(one,abs(data(itime,k)-data(itime,l))/dfact1) + dup(k)=dup(k)+one-tfact*tfact*(one-dfact) + dup(l)=dup(l)+one-tfact*tfact*(one-dfact) + end if + end do + end do + + call dtime_setup() + do i=1,nobs + dtime=data(itime,i) + call dtime_check(dtime, in_curbin, in_anybin) + if(.not.in_anybin) cycle + + if(in_curbin) then + dlat=data(ilat,i) + dlon=data(ilon,i) + + dpres=data(ipres,i) + error=data(ier2,i) + ikx=nint(data(ikxx,i)) + endif + +! Link observation to appropriate observation bin + if (nobs_bins>1) then + ibin = NINT( dtime/hr_obsbin ) + 1 + else + ibin = 1 + endif + IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin + + if (luse_obsdiag) my_diagLL => odiagLL(ibin) + +! Link obs to diagnostics structure + if (luse_obsdiag) then + my_diag => obsdiagLList_nextNode(my_diagLL ,& + create = .not.lobsdiag_allocated ,& + idv = is ,& + iob = ioid(i) ,& + ich = 1 ,& + elat = data(ilate,i) ,& + elon = data(ilone,i) ,& + luse = luse(i) ,& + miter = miter ) + + if(.not.associated(my_diag)) call die(myname, & + 'obsdiagLList_nextNode(), create =', .not.lobsdiag_allocated) + endif + + if(.not.in_curbin) cycle + +! Load obs error and u,v obs + obserror = max(cermin(ikx),min(cermax(ikx),data(ier,i))) + uob = data(iuob,i) + vob = data(ivob,i) + + + gnssrspdob=sqrt(uob*uob+vob*vob) + call tintrp2a1(ges_tv,tges,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + + factw = data(iff10,i) + if(sfcmod_gfs .or. sfcmod_mm5)then + sfcr = data(isfcr,i) + skint = data(iskint,i) + isli=data(idomsfc,i) + call comp_fact10(dlat,dlon,dtime,skint,sfcr,isli,mype,factw) + end if + + nty=ictype(ikx) + +! Initialize logical + z_height = .false. + +! Process observations reported with height differently than those +! reported with pressure. Type 260=nacelle 261=tower wind spd are +! encoded in NCEP prepbufr files with geometric height above +! sea level. + +! Process cygnss observations at mslp + if ( nty == 600 ) then + z_height = .true. + data(ihgt,i) = ten + endif + + if (z_height) then + + drpx = zero + dpres = data(ihgt,i) + dstn = data(istnelv,i) + call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + +! Get guess surface elevation and geopotential height profile +! at observation location. + call tintrp2a1(geop_hgtl,zges,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + +! Convert observation height (in dpres) from meters to grid relative +! units. Save the observation height in zob for later use. + zob = dpres + call grdcrd1(dpres,zges,nsig,1) + factw=one + rlow=zero + rhgh=zero + +! Compute observation pressure (only used for diagnostics) +! Set indices of model levels below (k1) and above (k2) observation. + if (dpresnsig) then + z1=zges(nsig-1); p1=prsltmp(nsig-1) + z2=zges(nsig); p2=prsltmp(nsig) + drpx = 1.e6_r_kind + else + k=dpres + k1=min(max(1,k),nsig) + k2=max(1,min(k+1,nsig)) + z1=zges(k1); p1=prsltmp(k1) + z2=zges(k2); p2=prsltmp(k2) + endif + + dz21 = z2-z1 + dlnp21 = p2-p1 + dz = zob-z1 + pobl = p1 + (dlnp21/dz21)*dz + presw = ten*exp(pobl) + + prsfc=psges + prsln2=log(exp(prsltmp(1))/prsfc) + sfcchk=log(psges) + if(dpres <= prsln2)then + factw=one + else + dx10=-goverrd*ten/tges(1) + if (dpres < dx10)then + factw=(dpres-dx10+factw*(prsln2-dpres))/(prsln2-dx10) + end if + end if + +! Put obs pressure in correct units to get grid coord. number + dpres=log(exp(dpres)*prsfc) + call grdcrd1(dpres,prsltmp(1),nsig,-1) + +! Get approx k value of sfc by using surface pressure of 1st ob + call grdcrd1(sfcchk,prsltmp(1),nsig,-1) + + +! Check to see if observations is below what is seen to be the surface + rlow=max(sfcchk-dpres,zero) + + rhgh=max(dpres-r0_001-rsig,zero) + + endif ! end of process observations with reported pressure + + if(luse(i))then + awork(1) = awork(1) + one + if(rlow/=zero) awork(2) = awork(2) + one + if(rhgh/=zero) awork(3) = awork(3) + one + end if + + ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+four*rlow) + + +! Interpolate guess u and v to observation location and time. + call tintrp31(ges_u,ugesin,dlat,dlon,dpres,dtime, & + hrdifsig,mype,nfldsig) + call tintrp31(ges_v,vgesin,dlat,dlon,dpres,dtime, & + hrdifsig,mype,nfldsig) + + +! Apply 10-meter wind reduction factor to guess winds. Compute +! guess wind speed. + ugesin=factw*ugesin + vgesin=factw*vgesin + gnssrspdges=sqrt(ugesin*ugesin+vgesin*vgesin) + + iz = max(1, min( int(dpres), nsig)) + delz = max(zero, min(dpres - float(iz), one)) + + if (save_jacobian) then + + u_ind = getindex(svars3d, 'u') + if (u_ind < 0) then + print *, 'Error: no variable u in state vector. Exiting.' + call stop2(1300) + endif + v_ind = getindex(svars3d, 'v') + if (v_ind < 0) then + print *, 'Error: no variable v in state vector. Exiting.' + call stop2(1300) + endif + + dhx_dx%st_ind(1) = iz + sum(levels(1:u_ind-1)) + dhx_dx%end_ind(1) = min(iz + 1,nsig) + sum(levels(1:u_ind-1)) + + dhx_dx%val(1) = (one - delz) * two * ugesin + dhx_dx%val(2) = delz * two * ugesin + + dhx_dx%st_ind(2) = iz + sum(levels(1:v_ind-1)) + dhx_dx%end_ind(2) = min(iz + 1,nsig) + sum(levels(1:v_ind-1)) + + dhx_dx%val(3) = (one - delz) * two * vgesin + dhx_dx%val(4) = delz * two * ugesin + endif + + + ddiff = gnssrspdob-gnssrspdges + +! write(6,*) "cygnss gnssrspdob=",gnssrspdob," gnssrspdges=",gnssrspdges," ddiff=",ddiff + + error=one/error + +! Check to see if observations is above the top of the model (regional mode) + if (dpres>rsig) ratio_errors=zero + + +! Gross error checks + obserror = one/max(ratio_errors*error,tiny_r_kind) + obserrlm = max(cermin(ikx),min(cermax(ikx),obserror)) + residual = abs(ddiff) + ratio = residual/obserrlm + if (ratio>cgross(ikx) .or. ratio_errors < tiny_r_kind) then + if (luse(i)) awork(4) = awork(4)+one + error = zero + ratio_errors = zero + muse(i)=.false. + else + ratio_errors=ratio_errors/sqrt(dup(i)) + end if + + if (ratio_errors*error <=tiny_r_kind) muse(i)=.false. + if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i)) + +! Compute penalty terms (linear & nonlinear qc). + val = error*ddiff + val2 = val*val + exp_arg = -half*val2 + rat_err2 = ratio_errors**2 + if (cvar_pg(ikx) > tiny_r_kind .and. error > tiny_r_kind) then + arg = exp(exp_arg) + wnotgross= one-cvar_pg(ikx) + cg_gnssrspd=cvar_b(ikx) + wgross = cg_term*cvar_pg(ikx)/(cg_gnssrspd*wnotgross) + term = log((arg+wgross)/(one+wgross)) + wgt = one-wgross/(arg+wgross) + wgt = wgt/wgtlim + else + term = exp_arg + wgt = wgtlim + rwgt = wgt/wgtlim + endif + valqc = -two*rat_err2*term + + +! Accumulate statistics for obs belonging to this task + if (luse(i) .and. muse(i)) then + if(rwgt < one) awork(61) = awork(61)+one + awork(5)=awork(5) + val2*rat_err2 + awork(6)=awork(6) + one + awork(22)=awork(22) + valqc + end if + +! Loop over pressure level groupings and obs to accumulate statistics +! as a function of observation type. + do k = 1,npres_print + if(luse(i) .and.presw >ptop(k) .and. presw<=pbot(k))then + ress = scale*ddiff + ressw = ress*ress + nn=1 + if (.not. muse(i)) then + nn=2 + if(ratio_errors*error >=tiny_r_kind)nn=3 + end if + bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count + bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+ddiff ! bias + bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw ! (o-g)**2 + bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2 ! penalty + bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty + + end if + end do + + !write(6,*) "cygnss: i,luse(i),muse(i)=",i,luse(i),muse(i) + + if (luse_obsdiag) then + call obsdiagNode_set(my_diag, wgtjo=(error*ratio_errors)**2, & + jiter=jiter, muse=muse(i), nldepart=gnssrspdob-sqrt(ugesin*ugesin+vgesin*vgesin)) + endif + +! If obs is "acceptable", load array with obs info for use +! in inner loop minimization (int* and stp* routines) + if (.not. last .and. muse(i)) then + + allocate(my_head) + call gnssrspdNode_appendto(my_head,gnssrspdhead(ibin)) + + my_head%idv = is + my_head%iob = ioid(i) + my_head%elat= data(ilate,i) + my_head%elon= data(ilone,i) + +! Set (i,j) indices of guess gridpoint that bound obs location + call get_ij(mm1,dlat,dlon,my_head%ij,my_head%wij) + + my_head%factw=factw + do j=1,4 + my_head%wij(j)=factw*my_head%wij(j) + end do + my_head%raterr2= ratio_errors**2 + my_head%res = gnssrspdob + my_head%uges = ugesin + my_head%vges = vgesin + my_head%err2 = error**2 + my_head%time = dtime + my_head%luse = luse(i) + my_head%b = cvar_b(ikx) + my_head%pg = cvar_pg(ikx) + + if (luse_obsdiag) then + call obsdiagNode_assert(my_diag, my_head%idv,my_head%iob,1,myname,'my_diag:my_head') + my_head%diags => my_diag + endif + + my_head => null() + end if +! Save select output for diagnostic file + if(conv_diagsave .and. luse(i))then + ii=ii+1 + gnssrspdob0 = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i)) + rstation_id = data(id,i) + err_input = data(ier2,i) + err_adjst = data(ier,i) + if (ratio_errors*error>tiny_r_kind) then + err_final = one/(ratio_errors*error) + else + err_final = huge_single + endif + + errinv_input = huge_single + errinv_adjst = huge_single + errinv_final = huge_single + if (err_input>tiny_r_kind) errinv_input = one/err_input + if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst + if (err_final>tiny_r_kind) errinv_final = one/err_final + + if(binary_diag) call contents_binary_diag_(my_diag) + if(netcdf_diag) call contents_netcdf_diag_(my_diag) + + end if + +! End of loop over observations + end do + +! Release memory of local guess arrays + call final_vars_ + +! Write information to diagnostic file + if(conv_diagsave) then + if(netcdf_diag) call nc_diag_write + if(binary_diag .and. ii>0)then + write(7)'gnssrspd',nchar,nreal,ii,mype,ioff0 + write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + deallocate(cdiagbuf,rdiagbuf) + + if (twodvar_regional) then + write(7)cprvstg(1:ii),csprvstg(1:ii) + deallocate(cprvstg,csprvstg) + endif + end if + end if + +! End of routine + + return + contains + + subroutine check_vars_ (proceed) + logical,intent(inout) :: proceed + integer(i_kind) ivar, istatus +! Check to see if required guess fields are available + call gsi_metguess_get ('var::ps', ivar, istatus ) + proceed=ivar>0 + call gsi_metguess_get ('var::z' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::u' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::v' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::tv', ivar, istatus ) + proceed=proceed.and.ivar>0 + end subroutine check_vars_ + + subroutine init_vars_ + + real(r_kind),dimension(:,: ),pointer:: rank2=>NULL() + real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL() + character(len=5) :: varname + integer(i_kind) ifld, istatus + +! If require guess vars available, extract from bundle ... + if(size(gsi_metguess_bundle)==nfldsig) then +! get ps ... + varname='ps' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_ps))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig)) + ges_ps(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_ps(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get z ... + varname='z' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_z))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) + ges_z(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_z(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get u ... + varname='u' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_u))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_u(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_u(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get v ... + varname='v' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_v))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_v(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_v(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get tv ... + varname='tv' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_tv))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_tv(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_tv(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + else + write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& + nfldsig,size(gsi_metguess_bundle) + call stop2(999) + endif + end subroutine init_vars_ + + subroutine init_netcdf_diag_ + character(len=80) string + character(len=128) diag_conv_file + integer(i_kind) ncd_fileid,ncd_nobs + logical append_diag + logical,parameter::verbose=.false. + write(string,900) jiter +900 format('conv_gnssrspd_',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 ) + if (save_jacobian) then + call nc_diag_header("jac_nnz", nnz) + call nc_diag_header("jac_nind", nind) + endif + endif + end subroutine init_netcdf_diag_ + subroutine contents_binary_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag + cdiagbuf(ii) = station_id ! station id + + rdiagbuf(1,ii) = ictype(ikx) ! observation type + rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype + + rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) + rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) + rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) + rdiagbuf(6,ii) = presw ! observation pressure (hPa) + rdiagbuf(7,ii) = data(ihgt,i) ! observation height (meters) + rdiagbuf(8,ii) = dtime-time_offset ! obs time (hours relative to analysis time) + + rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark + rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark + rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag + if(muse(i)) then + rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) + else + rdiagbuf(12,ii) = -one + endif + + rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight + rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (m/s)**-1 + rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (m/s)**-1 + rdiagbuf(16,ii) = errinv_final ! final inverse observation error (m/s)**-1 + + rdiagbuf(17,ii) = gnssrspdob ! wind speed observation (m/s) + rdiagbuf(18,ii) = ddiff ! obs-ges used in analysis (m/s) + rdiagbuf(19,ii) = gnssrspdob0-gnssrspdges ! obs-ges w/o bias correction (m/s) (future slot) + + rdiagbuf(20,ii) = factw ! 10m wind reduction factor + + rdiagbuf(21,ii) = 1.e+10_r_single ! ges ensemble spread (filled in by EnKF) + + ioff=ioff0 + if (lobsdiagsave) then + do jj=1,miter + ioff=ioff+1 + if (odiag%muse(jj)) then + rdiagbuf(ioff,ii) = one + else + rdiagbuf(ioff,ii) = -one + endif + enddo + do jj=1,miter+1 + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%nldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%tldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%obssen(jj) + enddo + endif + + if (twodvar_regional) then + ioff = ioff + 1 + rdiagbuf(ioff,ii) = data(idomsfc,i) ! dominate surface type + ioff = ioff + 1 + rdiagbuf(ioff,ii) = data(izz,i) ! model terrain at ob location + r_prvstg = data(iprvd,i) + cprvstg(ii) = c_prvstg ! provider name + r_sprvstg = data(isprvd,i) + csprvstg(ii) = c_sprvstg ! subprovider name + endif + + if (save_jacobian) then + call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,ii)) + ioff = ioff + size(dhx_dx) + endif + + end subroutine contents_binary_diag_ + subroutine contents_netcdf_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag +! Observation class + character(8),parameter :: obsclass = 'gnssrspd' + real(r_kind),dimension(miter) :: obsdiag_iuse + call nc_diag_metadata("Station_ID", station_id ) + call nc_diag_metadata("Observation_Class", obsclass ) + call nc_diag_metadata("Observation_Type", ictype(ikx) ) + call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) + call nc_diag_metadata("Latitude", sngl(data(ilate,i)) ) + call nc_diag_metadata("Longitude", sngl(data(ilone,i)) ) + call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) ) + call nc_diag_metadata("Pressure", sngl(presw) ) + call nc_diag_metadata("Height", sngl(data(ihgt,i)) ) + call nc_diag_metadata("Time", sngl(dtime-time_offset)) + call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) + call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) +! call nc_diag_metadata("Nonlinear_QC_Var_Jb", var_jb ) + call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", sngl(rwgt) ) + if(muse(i)) then + call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) + else + call nc_diag_metadata("Analysis_Use_Flag", sngl(-one) ) + endif + + call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) + call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) + call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) + + call nc_diag_metadata("Observation", sngl(gnssrspdob) ) + call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ddiff) ) + call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(gnssrspdob0-gnssrspdges) ) + + if (lobsdiagsave) then + do jj=1,miter + if (odiag%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", odiag%nldepart ) + call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart ) + call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) + endif + + if (twodvar_regional) then + call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) ) + call nc_diag_metadata("Model_Terrain", data(izz,i) ) + r_prvstg = data(iprvd,i) + call nc_diag_metadata("Provider_Name", c_prvstg ) + r_sprvstg = data(isprvd,i) + call nc_diag_metadata("Subprovider_Name", c_sprvstg ) + endif + if (save_jacobian) then + call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind) + call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind) + call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single)) + endif + + + end subroutine contents_netcdf_diag_ + + subroutine final_vars_ + if(allocated(ges_tv)) deallocate(ges_tv) + if(allocated(ges_v )) deallocate(ges_v ) + if(allocated(ges_u )) deallocate(ges_u ) + if(allocated(ges_z )) deallocate(ges_z ) + if(allocated(ges_ps)) deallocate(ges_ps) + end subroutine final_vars_ + +end subroutine setupgnssrspd +end module gnssrspd_setup diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 8075956431..6c00de5803 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -101,6 +101,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! 2019-03-15 Ladwig - add option for cloud analysis in observer ! 2019-03-28 Ladwig - add metar cloud obs as pseudo water vapor in var analysis ! 2020-09-08 CAPS(G. Zhao) - add 'l_use_dbz_directDA' flag not to sort obsdiag +! 2023-03-20 K Apodaca - add GNSS-R L2 Ocean Wind Speed ! ! input argument list: ! ndata(*,1)- number of prefiles retained for further processing @@ -165,7 +166,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) use m_rhs, only: toss_gps_sub => rhs_toss_gps use m_rhs, only: i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp use m_rhs, only: i_dbz use m_rhs, only: i_fed @@ -625,7 +626,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! Compute and print statistics for "conventional" data call statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz, & size(awork1,2),bwork1,awork1,ndata) diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90 index fc105515ff..b6bc849a37 100644 --- a/src/gsi/statsconv.f90 +++ b/src/gsi/statsconv.f90 @@ -1,6 +1,6 @@ subroutine statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,& i_swcp,i_lwcp,i_fed,i_dbz,i_ref,bwork,awork,ndata) !$$$ subprogram documentation block @@ -44,6 +44,7 @@ subroutine statsconv(mype,& ! 2015-07-10 pondeca - add cldch ! 2016-05-05 pondeca - add uwnd10m, vwnd10m ! 2017-05-12 Y. Wang and X. Wang - add dbz, POC: xuguang.wang@ou.edu +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed ! ! input argument list: ! mype - mpi task number @@ -62,6 +63,7 @@ subroutine statsconv(mype,& ! i_vis - index in awork array holding vis info ! i_pblh - index in awork array holding pblh info ! i_wspd10m- index in awork array holding wspd10m info +! i_gnssrspd - index in awork array holding gnssrspd info ! i_td2m - index in awork array holding td2m info ! i_mxtm - index in awork array holding mxtm info ! i_mitm - index in awork array holding mitm info @@ -94,13 +96,13 @@ subroutine statsconv(mype,& use constants, only: zero,three,five use obsmod, only: iout_sst,iout_pw,iout_t,iout_rw,iout_dw,& iout_uv,iout_gps,iout_ps,iout_q,iout_tcp,iout_lag,& - iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_td2m,& + iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_gnssrspd,iout_td2m,& iout_mxtm,iout_mitm,iout_pmsl,iout_howv,iout_tcamt,iout_lcbas,iout_cldch,& iout_uwnd10m,iout_vwnd10m,& iout_fed,iout_dbz,iout_swcp,iout_lwcp,& mype_dw,mype_rw,mype_sst,mype_gps,mype_uv,mype_ps,& mype_t,mype_pw,mype_q,mype_tcp,ndat,dtype,mype_lag,mype_gust,& - mype_vis,mype_pblh,mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,& + mype_vis,mype_pblh,mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,& mype_pmsl,mype_howv,mype_tcamt,mype_lcbas,mype_cldch,mype_uwnd10m,mype_vwnd10m,& mype_fed,mype_dbz,mype_swcp,mype_lwcp use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq @@ -112,7 +114,7 @@ subroutine statsconv(mype,& ! Declare passed variables integer(i_kind) ,intent(in ) :: mype,i_ps,i_uv,& i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag,i_gust,i_vis,i_pblh,& - i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,& + i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,& i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz,i_ref real(r_kind),dimension(7*nsig+100,i_ref) ,intent(in ) :: awork real(r_kind),dimension(npres_print,nconvtype,5,3),intent(in ) :: bwork @@ -122,7 +124,7 @@ subroutine statsconv(mype,& character(100) mesage integer(i_kind) numgrspw,numsst,nsuperp,nump,nhitopo,ntoodif - integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numuwnd10m,numvwnd10m + integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numgnssrspd,numuwnd10m,numvwnd10m integer(i_kind) numtd2m,nummxtm,nummitm,numpmsl,numhowv,numtcamt,numlcbas,numcldch integer(i_kind) numgrsswcp,numgrslwcp integer(i_kind) ntot,numlow,k,numssm,i,j @@ -677,6 +679,49 @@ subroutine statsconv(mype,& end if end if +! Summary report for conventional gnssrspd +if (mype == mype_gnssrspd) then + nread = 0 + nkeep = 0 + do i = 1, ndat + if (dtype(i) == 'gnssrspd') then + nread = nread + ndata(i, 2) + nkeep = nkeep + ndata(i, 3) + end if + end do + if (nread > 0) then + if (first) then + open(iout_gnssrspd) + else + open(iout_gnssrspd, position = 'append') + end if + + numgnssrspd = nint(awork(5, i_gnssrspd)) + pw = zero + pw3 = zero + if (nkeep > 0) then + mesage = 'current fit of conventional gnssrspd data, ranges in m/s' + do j = 1, nconvtype + pflag(j) = trim(ioctype(j)) == 'gnssrspd' + end do + call dtast(bwork, 1, pbotall, ptopall, mesage, jiter, iout_gnssrspd, pflag) + + numgross = nint(awork(6, i_gnssrspd)) + numfailqc = nint(awork(21, i_gnssrspd)) + if (numgnssrspd > 0) then + pw = awork(4, i_gnssrspd) / numgnssrspd + pw3 = awork(22, i_gnssrspd) / numgnssrspd + end if + write(iout_gnssrspd, 925) 'gnssrspd', numgross, numfailqc + end if + write(iout_gnssrspd, 950) 'gnssrspd', jiter, nread, nkeep, numgnssrspd + write(iout_gnssrspd, 951) 'gnssrspd', awork(4, i_gnssrspd), awork(22, i_gnssrspd), pw, pw3 + + close(iout_gnssrspd) + end if +end if + + ! Summary report for conventional td2m if(mype==mype_td2m) then nread=0 diff --git a/src/gsi/stpgnssrspd.f90 b/src/gsi/stpgnssrspd.f90 new file mode 100644 index 0000000000..3909c40534 --- /dev/null +++ b/src/gsi/stpgnssrspd.f90 @@ -0,0 +1,175 @@ +module stpgnssrspdmod + + +!$$$ module documentation block +! . . . . +! module: stpgnssrspdmod module for stpgnssrspd and its tangent linear stpgnssrspd_tl +! prgmmr: K. Apodaca org: Spire Global, Inc. date: 2022-03-12 +! Largely based on other stp_* routines +! +! abstract: module for stpgnssrspd and its tangent linear stpgnssrspd_tl +! +! program history log: +! 2023-09-21 K. Apodaca - add documentation +! subroutine included: +! sub stpgnssrspd +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +implicit none + +PRIVATE +PUBLIC stpgnssrspd + +contains + +subroutine stpgnssrspd(gnssrspdhead,rval,sval,out,sges,nstep) +!$$$ subprogram documentation block +! . . . . +! subprogram: stpgnssrspd calculate penalty and stepsize terms +! for wind speed, with nonlinear qc. +! 2023-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed +! +! abstract: calculate penalty and stepsize terms for wind speed +! +! program history log: + +! +! input argument list: +! gnssrspdhead +! ru - search direction for u +! rv - search direction for v +! su - analysis increment for u +! sv - analysis increment for v +! sges - step size estimates (nstep) +! nstep - number of stepsizes (==0 means use outer iteration values) +! +! output argument list +! out(1:nstep) - contribution to penalty from wind speed sges(1:nstep) +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,i_kind,r_quad + use qcmod, only: nlnqc_iter,varqc_iter + use constants, only: zero,half,one,two,tiny_r_kind,cg_term,zero_quad,r3600 + use gsi_4dvar, only: ltlint + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use m_obsNode, only: obsNode + use m_gnssrspdNode, only: gnssrspdNode + use m_gnssrspdNode, only: gnssrspdNode_typecast + use m_gnssrspdNode, only: gnssrspdNode_nextcast + implicit none + +! Declare passed variables + class(obsNode), pointer ,intent(in ) :: gnssrspdhead + integer(i_kind) ,intent(in ) :: nstep + real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges + real(r_quad),dimension(max(1,nstep)),intent(inout) :: out + type(gsi_bundle) ,intent(in ) :: rval,sval + +! Declare local variables + integer(i_kind) j1,j2,j3,j4,kk,ier,istatus + real(r_kind) w1,w2,w3,w4,time_gnssrspd + real(r_kind) valu,valv,ucur,vcur,gnssrspdnl,gnssrspdtl,uu,vv,gnssrspd + real(r_kind),dimension(max(1,nstep)):: pen + real(r_kind) cg_gnssrspd,pencur,wgross,wnotgross + real(r_kind) pg_gnssrspd,pentl + real(r_kind),pointer,dimension(:) :: su,sv + real(r_kind),pointer,dimension(:) :: ru,rv + type(gnssrspdNode), pointer :: gnssrspdptr + + out=zero_quad + +! If no gnssrspd data return + if(.not. associated(gnssrspdhead))return + + time_gnssrspd=zero +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier + if(ier/=0)return + + gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead) + do while (associated(gnssrspdptr)) + + if(gnssrspdptr%luse)then + if(nstep > 0)then + j1 = gnssrspdptr%ij(1) + j2 = gnssrspdptr%ij(2) + j3 = gnssrspdptr%ij(3) + j4 = gnssrspdptr%ij(4) + w1 = gnssrspdptr%wij(1) + w2 = gnssrspdptr%wij(2) + w3 = gnssrspdptr%wij(3) + w4 = gnssrspdptr%wij(4) + + valu=w1* ru(j1)+w2* ru(j2)+w3* ru(j3)+w4* ru(j4) + valv=w1* rv(j1)+w2* rv(j2)+w3* rv(j3)+w4* rv(j4) + ucur=w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4)+gnssrspdptr%uges + vcur=w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4)+gnssrspdptr%vges + + if (ltlint) then + gnssrspd=sqrt(ucur*ucur+vcur*vcur)-gnssrspdptr%res + pencur=gnssrspd*gnssrspd*gnssrspdptr%err2 + do kk=1,nstep + gnssrspdnl=sqrt(ucur*ucur+vcur*vcur) + gnssrspdtl=ucur*valu+vcur*valv + if (gnssrspdnl>tiny_r_kind*100._r_kind) then + gnssrspdtl=gnssrspdtl/gnssrspdnl + else + gnssrspdtl=zero + endif + pentl =two*gnssrspdtl*gnssrspd*gnssrspdptr%err2 + pen(kk)=pencur+sges(kk)*pentl + end do + else + do kk=1,nstep + uu=ucur+sges(kk)*valu + vv=vcur+sges(kk)*valv + gnssrspd=sqrt(uu*uu+vv*vv)-gnssrspdptr%res + pen(kk)=gnssrspd*gnssrspd*gnssrspdptr%err2 + end do + end if + else + pen(1)=gnssrspdptr%res*gnssrspdptr%res*gnssrspdptr%err2 + end if + +! Modify penalty term if nonlinear QC + if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and. & + gnssrspdptr%b > tiny_r_kind) then + pg_gnssrspd=gnssrspdptr%pg*varqc_iter + cg_gnssrspd=cg_term/gnssrspdptr%b + wnotgross= one-pg_gnssrspd + wgross = pg_gnssrspd*cg_gnssrspd/wnotgross + pencur = -two*log((exp(-half*pencur) + wgross)/(one+wgross)) + do kk=1,max(1,nstep) + pen(kk) = -two*log((exp(-half*pen(kk) ) + wgross)/(one+wgross)) + enddo + endif + + out(1) = out(1)+pen(1)*gnssrspdptr%raterr2 + do kk=2,nstep + out(kk) = out(kk)+(pen(kk)-pen(1))*gnssrspdptr%raterr2 + end do + + end if + + gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr) + + end do + return +end subroutine stpgnssrspd + +end module stpgnssrspdmod