diff --git a/src/gsi/gsi_fedOper.F90 b/src/gsi/gsi_fedOper.F90 new file mode 100644 index 0000000000..b2b2400ff0 --- /dev/null +++ b/src/gsi/gsi_fedOper.F90 @@ -0,0 +1,174 @@ +module gsi_fedOper +!$$$ subprogram documentation block +! . . . . +! subprogram: module gsi_fedOper +! +! abstract: an obOper extension for fedNode type +! +! program history log: +! 2023-07-10 D. Dowell - created new module for FED (flash extent +! density); gsi_dbzOper.F90 code used as a +! starting point for developing this new module +! +! 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_fedNode , only: fedNode + implicit none + public:: fedOper ! data structure + public:: diag_fed + + type,extends(obOper):: fedOper + contains + procedure,nopass:: mytype + procedure,nopass:: nodeMold + procedure:: setup_ + procedure:: intjo1_ + procedure:: stpjo1_ + end type fedOper + +! def diag_fed- namelist logical to compute/write (=true) FED diag files + logical,save:: diag_fed=.false. + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*),parameter :: myname='gsi_fedOper' + type(fedNode),save,target:: myNodeMold_ + +contains + function mytype(nodetype) + implicit none + character(len=:),allocatable:: mytype + logical,optional, intent(in):: nodetype + mytype="[fedOper]" + 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 fed_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_fed + + use obsmod , only: write_diag + use jfunc , only: jiter + + use mpeu_util, only: die + + use obsmod, only: dirname, ianldate + + implicit none + class(fedOper ), 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 + integer(i_kind):: lu_diag + character(128):: diag_file + character(80):: string + + if(nobs == 0) then + + if( (mype == 0) .and. init_pass ) then + write(string,600) jiter +600 format('fed_',i2.2) + diag_file=trim(dirname) // trim(string) + write(6,*) 'write ianldate to ', diag_file + open(newunit=lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + write(lu_diag) ianldate + close(lu_diag) + endif + + return + + endif + + 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_fed + + call setup(self%obsLL(:), self%odiagLL(:), & + lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave,init_pass) + + end subroutine setup_ + + subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) + 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(fedOper ),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 + + end subroutine intjo1_ + + subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) + 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(fedOper ),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 + + end subroutine stpjo1_ + +end module gsi_fedOper diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index b98cd2d0da..b514e11c1e 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -217,6 +217,7 @@ gsi_colvkOper.F90 gsi_dbzOper.F90 gsi_dwOper.F90 gsi_enscouplermod.f90 +gsi_fedOper.F90 gsi_gpsbendOper.F90 gsi_gpsrefOper.F90 gsi_gustOper.F90 @@ -338,6 +339,7 @@ m_distance.f90 m_dtime.F90 m_dwNode.F90 m_extOzone.F90 +m_fedNode.F90 m_find.f90 m_gpsNode.F90 m_gpsrhs.F90 @@ -478,6 +480,7 @@ read_cris.f90 read_dbz_nc.f90 read_dbz_netcdf.f90 read_diag.f90 +read_fed.f90 read_files.f90 read_fl_hdob.f90 read_gfs_ozone_for_regional.f90 @@ -532,6 +535,7 @@ setupco.f90 setupdbz.f90 setupdbz_lib.f90 setupdw.f90 +setupfed.f90 setupgust.f90 setuphowv.f90 setuplag.f90 diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index 5df899825a..6db7921905 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -66,6 +66,7 @@ module gsi_obOperTypeManager use gsi_lightOper , only: lightOper use gsi_dbzOper , only: dbzOper + use gsi_fedOper , only: fedOper use gsi_cldtotOper , only: cldtotOper use kinds , only: i_kind @@ -136,6 +137,7 @@ module gsi_obOperTypeManager public:: iobOper_lwcp public:: iobOper_light public:: iobOper_dbz + public:: iobOper_fed public:: iobOper_cldtot enum, bind(C) @@ -181,6 +183,7 @@ module gsi_obOperTypeManager enumerator:: iobOper_lwcp enumerator:: iobOper_light enumerator:: iobOper_dbz + enumerator:: iobOper_fed enumerator:: iobOper_cldtot enumerator:: iobOper_extra_ @@ -242,6 +245,7 @@ module gsi_obOperTypeManager type( lwcpOper), target, save:: lwcpOper_mold type( lightOper), target, save:: lightOper_mold type( dbzOper), target, save:: dbzOper_mold + type( fedOper), target, save:: fedOper_mold type( cldtotOper), target, save:: cldtotOper_mold contains @@ -390,6 +394,7 @@ function dtype2index_(dtype) result(index_) case("goes_glm" ); index_= iobOper_light case("dbz" ,"[dbzoper]" ); index_= iobOper_dbz + case("fed" ,"[fedoper]" ); index_= iobOper_fed case("cldtot" ,"[cldtotoper]" ); index_= iobOper_cldtot case("mta_cld" ); index_= iobOper_cldtot @@ -487,6 +492,7 @@ function index2vmold_(iobOper) result(vmold_) case(iobOper_lwcp ); vmold_ => lwcpOper_mold case(iobOper_light ); vmold_ => lightOper_mold case(iobOper_dbz ); vmold_ => dbzOper_mold + case(iobOper_fed ); vmold_ => fedOper_mold case(iobOper_cldtot ); vmold_ => cldtotOper_mold case( obOper_undef ); vmold_ => null() @@ -602,6 +608,7 @@ subroutine cobstype_config_() cobstype(iobOper_lwcp ) ="lwcp " ! lwcp_ob_type cobstype(iobOper_light ) ="light " ! light_ob_type cobstype(iobOper_dbz ) ="dbz " ! dbz_ob_type + cobstype(iobOper_fed ) ="fed " ! fed_ob_type cobstype(iobOper_cldtot ) ="cldtot " ! using q_ob_type cobstype_configured_=.true. diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index cf885c2b64..de19c85fab 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -21,6 +21,7 @@ module gsimod lread_obs_save,lread_obs_skip,time_window_rad,tcp_posmatch,tcp_box, & neutral_stability_windfact_2dvar,use_similarity_2dvar,ta2tb use gsi_dbzOper, only: diag_radardbz + use gsi_fedOper, only: diag_fed use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,& @@ -560,6 +561,7 @@ module gsimod ! diag_co - logical to turn off or on the diagnostic carbon monoxide file (true=on) ! diag_light - logical to turn off or on the diagnostic lightning file (true=on) ! diag_radardbz - logical to turn off or on the diagnostic radar reflectivity file (true=on) +! diag_fed - logical to turn off or on the diagnostic flash extent density file (true=on) ! write_diag - logical to write out diagnostic files on outer iteration ! lobsdiagsave - write out additional observation diagnostics ! ltlint - linearize inner loop @@ -738,8 +740,8 @@ module gsimod min_offset,pseudo_q2,& iout_iter,npredp,retrieval,& tzr_qc,tzr_bufrsave,& - diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,iguess, & - write_diag,reduce_diag, & + diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,diag_fed, & + iguess,write_diag,reduce_diag, & oneobtest,sfcmodel,dtbduv_on,ifact10,l_foto,offtime_data,& use_pbl,use_compress,nsig_ext,gpstop,commgpstop, commgpserrinf, & perturb_obs,perturb_fact,oberror_tune,preserve_restart_date, & @@ -1977,6 +1979,7 @@ subroutine gsimain_initialize diag_pcp=.false. diag_light=.false. diag_radardbz=.false. + diag_fed=.false. use_limit = 0 end if if(reduce_diag) use_limit = 0 diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90 index 91b811147e..a68355471b 100644 --- a/src/gsi/intjo.f90 +++ b/src/gsi/intjo.f90 @@ -31,7 +31,7 @@ module intjomod use gsi_obOperTypeManager, only: & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, & + iobOper_rw, iobOper_dbz, iobOper_fed, & 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 +60,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_rw, iobOper_dbz, iobOper_fed, & 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_fedNode.F90 b/src/gsi/m_fedNode.F90 new file mode 100644 index 0000000000..84a319cd12 --- /dev/null +++ b/src/gsi/m_fedNode.F90 @@ -0,0 +1,248 @@ +module m_fedNode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_fedNode +! prgmmr: YPW +! org: CIMMS +! date: 2019-09-24 +! +! abstract: class-module of obs-type fedNode (GLM flash extent density) +! Modified based on m_tdNode.f90 +! +! program history log: +! 2019-09-24 YPW - 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:: fedNode + + type,extends(obsNode):: fedNode + type(obs_diag), pointer :: diags => NULL() + real(r_kind) :: res ! flash extent density residual + real(r_kind) :: err2 ! flash extent density error squared + real(r_kind) :: raterr2 ! square of ratio of final obs error + ! to original obs error + real(r_kind) :: b ! variational quality control parameter + real(r_kind) :: pg ! variational quality control parameter + real(r_kind) :: jb ! variational quality control parameter + real(r_kind) :: wij(8) ! horizontal interpolation weights + real(r_kind) :: fedpertb ! random number adding to the obs + integer(i_kind) :: k1 ! level of errtable 1-33 + integer(i_kind) :: kx ! ob type + integer(i_kind) :: ij(8) ! horizontal locations + + real (r_kind) :: dlev ! reference to the vertical grid + contains + procedure,nopass:: mytype + procedure:: setHop => obsNode_setHop_ + procedure:: xread => obsNode_xread_ + procedure:: xwrite => obsNode_xwrite_ + procedure:: isvalid => obsNode_isvalid_ + procedure:: gettlddp => gettlddp_ + + end type fedNode + + public:: fedNode_typecast + public:: fedNode_nextcast + interface fedNode_typecast; module procedure typecast_ ; end interface + interface fedNode_nextcast; module procedure nextcast_ ; end interface + + public:: fedNode_appendto + interface fedNode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: MYNAME="m_fedNode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(aNode) result(ptr_) +!-- cast a class(obsNode) to a type(tdNode) + use m_obsNode, only: obsNode + implicit none + type(fedNode),pointer:: ptr_ + class(obsNode),pointer,intent(in):: aNode + ptr_ => null() + if(.not.associated(aNode)) return + select type(aNode) + type is(fedNode) + ptr_ => aNode + end select +return +end function typecast_ + +function nextcast_(aNode) result(ptr_) +!-- cast an obsNode_next(obsNode) to a type(fedNode) + use m_obsNode, only: obsNode,obsNode_next + implicit none + type(fedNode),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(fedNode),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="[fedNode]" +end function mytype + + +subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) + use m_obsdiagNode, only: obsdiagLookup_locate + implicit none + class(fedNode),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%jb , & + aNode%k1 , & + aNode%kx , & + aNode%dlev , & + aNode%wij , & + aNode%ij + if(istat/=0) then + call perr(myname_,'read(%(res,err2,...), iostat =',istat) + _EXIT_(myname_) + return + endif + + 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(fedNode),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%jb , & + aNode%k1 , & + aNode%kx , & + aNode%dlev , & + aNode%wij , & + aNode%ij + if(jstat/=0) then + call perr(myname_,'write(%res,err2,...), iostat =',jstat) + _EXIT_(myname_) + return + endif +_EXIT_(myname_) +return +end subroutine obsNode_xwrite_ + +subroutine obsNode_setHop_(aNode) + use m_cvgridLookup, only: cvgridLookup_getiw + implicit none + class(fedNode),intent(inout):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' +_ENTRY_(myname_) + call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%dlev,aNode%ij,aNode%wij) +_EXIT_(myname_) +return +end subroutine obsNode_setHop_ + +function obsNode_isvalid_(aNode) result(isvalid_) + implicit none + logical:: isvalid_ + class(fedNode),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(fedNode), 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_fedNode diff --git a/src/gsi/m_obsNodeTypeManager.F90 b/src/gsi/m_obsNodeTypeManager.F90 index b5ecc6e1ba..43b42e4bf2 100644 --- a/src/gsi/m_obsNodeTypeManager.F90 +++ b/src/gsi/m_obsNodeTypeManager.F90 @@ -70,6 +70,7 @@ module m_obsNodeTypeManager use m_lightNode, only: lightNode use m_dbzNode , only: dbzNode + use m_fedNode, only: fedNode use kinds, only: i_kind use m_obsNode, only: obsNode @@ -124,6 +125,7 @@ module m_obsNodeTypeManager public:: iobsNode_light public:: iobsNode_dbz + public:: iobsNode_fed public :: obsNode_typeMold public :: obsNode_typeIndex @@ -179,6 +181,7 @@ module m_obsNodeTypeManager type( lwcpNode), target, save:: lwcp_mold type( lightNode), target, save:: light_mold type( dbzNode), target, save:: dbz_mold + type( fedNode), target, save:: fed_mold !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ character(len=*),parameter :: myname='m_obsNodeTypeManager' @@ -245,6 +248,7 @@ module m_obsNodeTypeManager enumerator:: iobsNode_lwcp enumerator:: iobsNode_light enumerator:: iobsNode_dbz + enumerator:: iobsNode_fed enumerator:: iobsNode_extra_ end enum @@ -314,6 +318,7 @@ function vname2index_(vname) result(index_) case("light","[lightnode]"); index_ = iobsNode_light case("dbz" , "[dbznode]"); index_ = iobsNode_dbz + case("fed" , "[fednode]"); index_ = iobsNode_fed end select end function vname2index_ @@ -377,6 +382,7 @@ function vmold2index_select_(mold) result(index_) type is(lightNode); index_ = iobsNode_light type is( dbzNode); index_ = iobsNode_dbz + type is( fedNode); index_ = iobsNode_fed end select end function vmold2index_select_ @@ -434,6 +440,7 @@ function index2vmold_(i_obType) result(obsmold_) case(iobsNode_light); obsmold_ => light_mold case(iobsNode_dbz); obsmold_ => dbz_mold + case(iobsNode_fed); obsmold_ => fed_mold end select end function index2vmold_ diff --git a/src/gsi/m_rhs.F90 b/src/gsi/m_rhs.F90 index baee074688..aea417fe27 100644 --- a/src/gsi/m_rhs.F90 +++ b/src/gsi/m_rhs.F90 @@ -80,6 +80,7 @@ module m_rhs public:: i_lwcp public:: i_light public:: i_dbz + public:: i_fed public:: i_cldtot public:: awork_size @@ -146,6 +147,7 @@ module m_rhs enumerator:: i_lwcp enumerator:: i_light enumerator:: i_dbz + enumerator:: i_fed enumerator:: i_cldtot enumerator:: i_outbound diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 3066cdb5ca..a059586e67 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) +! 2023-07-10 Y. Wang, D. Dowell - add variables for flash extent density ! ! Subroutines Included: ! sub init_obsmod_dflts - initialize obs related variables to default values @@ -186,6 +187,7 @@ module obsmod ! def write_diag - namelist logical array to compute/write (=true) diag files ! def diag_radardbz- namelist logical to compute/write (=true) radar ! reflectiivty diag files +! def diag_fed - namelist logical to compute/write (=true) flash extent density diag files ! def reduce_diag - namelist logical to produce reduced radiance diagnostic files ! def use_limit - parameter set equal to -1 if diag files produced or 0 if not diag files or reduce_diag ! def obs_setup - prefix for files passing pe relative obs data to setup routines @@ -434,6 +436,7 @@ module obsmod public :: ran01dom,dval_use public :: iout_pcp,iout_rad,iadate,iadatemn,write_diag,reduce_diag,oberrflg,bflag,ndat,dthin,dmesh,l_do_adjoint public :: diag_radardbz + public :: diag_fed public :: lsaveobsens public :: iout_cldch, mype_cldch public :: nprof_gps,time_offset,ianldate,tcp_box @@ -483,7 +486,9 @@ module obsmod public :: iout_dbz, mype_dbz ! --- DBZ DA --- - + + public :: iout_fed, mype_fed + public :: obsmod_init_instr_table public :: obsmod_final_instr_table public :: nobs_sub @@ -583,12 +588,12 @@ module obsmod 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_uwnd10m,iout_vwnd10m + 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_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz + mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz, mype_fed integer(i_kind) mype_cldch integer(i_kind) iout_swcp, iout_lwcp integer(i_kind) mype_swcp, mype_lwcp @@ -638,6 +643,7 @@ module obsmod logical lobserver,l_do_adjoint, lobsdiag_forenkf logical,dimension(0:50):: write_diag logical diag_radardbz + logical diag_fed logical reduce_diag logical offtime_data logical hilbert_curve @@ -789,6 +795,7 @@ subroutine init_obsmod_dflts end do write_diag(1)=.true. diag_radardbz = .false. + diag_fed = .false. reduce_diag = .false. use_limit = -1 lobsdiagsave=.false. @@ -853,6 +860,7 @@ subroutine init_obsmod_dflts iout_lwcp=236 ! liquid-water content path iout_light=237 ! lightning iout_dbz=238 ! radar reflectivity + iout_fed=239 ! flash extent density mype_ps = npe-1 ! surface pressure mype_t = max(0,npe-2) ! temperature @@ -887,6 +895,7 @@ subroutine init_obsmod_dflts mype_lwcp=max(0,npe-31) ! liquid-water content path 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 ! Initialize arrays used in namelist obs_input diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 new file mode 100644 index 0000000000..c478b3d93f --- /dev/null +++ b/src/gsi/read_fed.f90 @@ -0,0 +1,525 @@ +subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) +!$$$ subprogram documentation block +! . . . . +! ABSTRACT: +! This routine reads in netcdf or prepbufr flash-extent density (FED) data. +! +! PROGRAM HISTORY LOG: +! 2018-07-25 Rong Kong (CAPS/OU) - modified based on read_radarref_mosaic.f90 +! 2019-09-20 Yaping Wang (CIMMS/OU) +! 2021-07-01 David Dowell (DCD; NOAA GSL) - added maximum flashes/min for observed FED +! +! input argument list: +! infile - unit from which to read observation information file +! obstype - observation type to process +! lunout - unit to which to write data for further processing +! twind - input group time window (hours) +! sis - observation variable name +! +! output argument list: +! nread - number of type "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: FORTRAN 90 +! +!$$$ +! +!_____________________________________________________________________ +! + use kinds, only: r_kind,r_double,i_kind + use constants, only: zero,one,deg2rad + use convinfo, only: nconvtype,ctwind,icuse,ioctype + use gsi_4dvar, only: l4dvar,l4densvar,winlen + use gridmod, only: tll2xy + use mod_wrfmass_to_a, only: wrfmass_obs_to_a8 + use mpimod, only: npe + use obsmod, only: perturb_obs,iadatemn + + use netcdf + implicit none + + include 'netcdf.inc' +! + character(len=*), intent(in) :: infile,obstype + integer(i_kind), intent(in) :: lunout + integer(i_kind), intent(inout) :: nread,ndata + integer(i_kind), intent(inout) :: nodata + integer(i_kind), dimension(npe) ,intent(inout) :: nobs + real(r_kind), intent(in ) :: twind + character(len=*), intent(in) :: sis + +! Declare local parameters + real(r_kind),parameter:: r90 = 90.0_r_kind + real(r_kind),parameter:: r360 = 360.0_r_kind + real(r_kind),parameter:: oe_fed = 1.0_r_kind + real(r_kind),parameter:: fed_lowbnd = 0.1_r_kind ! use fed == fed_lowbnd + real(r_kind),parameter:: fed_lowbnd2 = 0.1_r_kind ! use fed >= fed_lowbnd2 +! real(r_kind),parameter:: fed_highbnd = 18.0_r_kind ! 18 flashes/min from Sebok and Back (2021, unpublished) + real(r_kind),parameter:: fed_highbnd = 8.0_r_kind ! 8 flashes/min from Back (2023) for regional FV3 tests + +! +! For fed observations +! + integer(i_kind) nreal,nchanl + + integer(i_kind) ifn,i + + integer(i_kind) :: ilon,ilat + + logical :: fedobs, fedob + real(r_kind),allocatable,dimension(:,:):: cdata_out + real(r_kind) :: federr, thiserr + real(r_kind) :: hgt_fed(1) + data hgt_fed / 6500.0 / + + real(r_kind) :: i_maxloc,j_maxloc,k_maxloc + integer(i_kind) :: kint_maxloc + real(r_kind) :: fed_max + integer(i_kind) :: ndata2 + integer(i_kind) :: ppp + +! +! for read in bufr +! + real(r_kind) :: hdr(5),obs(1,3) + character(80):: hdrstr='SID XOB YOB DHR TYP' + character(80):: obsstr='FED' + + character(8) subset + character(8) station_id + real(r_double) :: rstation_id + equivalence(rstation_id,station_id) + integer(i_kind) :: lunin,idate + integer(i_kind) :: ireadmg,ireadsb + + integer(i_kind) :: maxlvl + integer(i_kind) :: numlvl,numfed,nmsgmax,maxobs + integer(i_kind) :: k,iret + integer(i_kind) :: nmsg,ntb + + real(r_kind),allocatable,dimension(:,:) :: fed3d_column ! 3D fed in column + real(r_kind),allocatable,dimension(:) :: utime ! time + + integer(i_kind) :: ikx + real(r_kind) :: timeo,t4dv + + character*128 :: myname='read_fed' + + real(r_kind) :: dlat, dlon ! rotated corrdinate + real(r_kind) :: dlat_earth, dlon_earth ! in unit of degree + real(r_kind) :: rlat00, rlon00 ! in unit of rad + + logical :: l_psot_fed + logical :: l_latlon_fedobs + logical :: outside + integer :: unit_table + +! for read netcdf + integer(i_kind) :: sec70,mins_an + integer(i_kind) :: varID, ncdfID, status + real(r_kind) :: timeb,twindm,rmins_an,rmins_ob + + + unit_table = 23 +!********************************************************************** +! +! END OF DECLARATIONS....start of program +! + write(6,*) "r_kind=",r_kind + l_psot_fed = .FALSE. + l_latlon_fedobs = .TRUE. + + fedob = obstype == 'fed' + if(fedob) then + nreal=25 + else + write(6,*) ' illegal obs type in read_fed : obstype=',obstype + call stop2(94) + end if + if(perturb_obs .and. fedob)nreal=nreal+1 + write(6,*)'read_fed: nreal=',nreal + + fedobs = .false. + ikx=0 + do i=1,nconvtype + if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then + fedobs=.true. + ikx=i + federr = oe_fed ! Obs error (flashes per minute) + thiserr = federr + exit ! Exit loop when finished with initial convinfo fields + else if (i == nconvtype ) then + write(6,*) 'read_fed: Obs Type for fed is not in CONVINFO !' + write(6,*) 'read_fed: PLEASE modify the CONVINFO file !' + write(6,*) 'read_fed: abort read_fed !' + return + endif + end do + write(6,'(1x,A,A30,I4,A15,F7.3,A7)') & + trim(myname),': fed in convinfo-->ikx=',ikx,' fed ob err:',thiserr," (fed)" + + nread=0 + ndata=0 + nchanl=0 + ifn = 15 + + if(fedobs) then + maxlvl= 1 ! fed only has one level + + if(trim(infile) .eq. "fedbufr") then ! prebufr or netcdf format + !! get message and subset counts + ! nmsgmax and maxobs are read in from BUFR data file, not pre-set. + call getcount_bufr(infile,nmsgmax,maxobs) + write(6,*)'read_fed: nmsgmax=',nmsgmax,' maxobs=',maxobs + +! read in fed obs in bufr code format + lunin = 10 + allocate(fed3d_column(maxlvl+2+2,maxobs)) + + open ( unit = lunin, file = trim(infile),form='unformatted',err=200) + call openbf ( lunin, 'IN', lunin ) + open(unit_table,file='prepobs_kr.bufrtable') !temporily dump the bufr table, which is already saved in file + call dxdump(lunin,unit_table) + call datelen ( 10 ) + + nmsg=0 + ntb = 0 + + ndata =0 + ppp = 0 + msg_report: do while (ireadmg(lunin,subset,idate) == 0) + nmsg=nmsg+1 + if (nmsg>nmsgmax) then + write(6,*)'read_fed: messages exceed maximum ',nmsgmax + call stop2(50) + endif + loop_report: do while (ireadsb(lunin) == 0) + ntb = ntb+1 + if (ntb>maxobs) then + write(6,*)'read_fed: reports exceed maximum ',maxobs + call stop2(50) + endif + + ! Extract type, date, and location information from BUFR file + call ufbint(lunin,hdr,5,1,iret,hdrstr) + if(hdr(3) .gt. r90 ) write(6,*) "Inside read_fed.f90, hdr(2)=",hdr(2),"hdr(3)=",hdr(3) + if ( l_latlon_fedobs ) then + if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_report + if(hdr(2)== r360)hdr(2)=hdr(2)-r360 + if(hdr(2) < zero)hdr(2)=hdr(2)+r360 + end if + +! check time window in subset + if (l4dvar.or.l4densvar) then + t4dv=hdr(4) + if (t4dvwinlen) then + write(6,*)'read_fed: time outside window ',& + t4dv,' skip this report' + cycle loop_report + endif + else + timeo=hdr(4) + if (abs(timeo)>ctwind(ikx) .or. abs(timeo) > twind) then + write(6,*)'read_fed: time outside window ',& + timeo,' skip this report' + cycle loop_report + endif + endif +! read in observations + call ufbint(lunin,obs,1,3,iret,obsstr) !Single level bufr data, Rong Kong + if(obs(1,1) .gt. 5 ) write(6,*) "Inside read_fed.f90, obs(1,1)=",obs(1,1) + numlvl=min(iret,maxlvl) + if (numlvl .ne. maxlvl) then + write(6,*)' read_fed: numlvl is not equalt to maxlvl:',numlvl,maxlvl + end if + if(hdr(3) .gt. 90) write(6,*) "hdr(3)=",hdr(3) + if ( l_latlon_fedobs ) then + if(hdr(2)>= r360)hdr(2)=hdr(2)-r360 + if(hdr(2) < zero)hdr(2)=hdr(2)+r360 + fed3d_column(1,ntb)=hdr(2) ! observation location, earth lon + fed3d_column(2,ntb)=hdr(3) ! observation location, earth lat +! write(6,*) "Inside read_fed.f90, fed3d_column(1,ntb)=",fed3d_column(1,ntb),"fed3d_column(2,ntb)=",fed3d_column(2,ntb) + else + fed3d_column(1,ntb)=hdr(2)*10.0_r_kind ! observation location, grid index i + fed3d_column(2,ntb)=hdr(3)*10.0_r_kind ! observation location, grid index j + end if + + if (l_psot_fed .and. .NOT. l_latlon_fedobs ) then + do k=1,numlvl + if (NINT(fed3d_column(1,ntb)) .eq. 175 .and. NINT(fed3d_column(2,ntb)) .eq. 105 .and. & + NINT(hgt_fed(k)) .ge. 100 ) then + write(6,*) 'read_fed: single point/column obs run on grid: 175 105' + write(6,*) 'read_fed: found the pseudo single(column) fed obs:',fed3d_column(1:2,ntb),hgt_fed(k) + else + obs(1,1) = -999.0 + end if + end do + end if + + fed3d_column(3,ntb)=obs(1,1) + fed3d_column(4,ntb)=obs(1,2) + fed3d_column(5,ntb)=obs(1,3) + if (obs(1,1) == fed_lowbnd .or. obs(1,1) >= fed_lowbnd2 ) then + if (obs(1,1) == 0.0) then + ppp = ppp + 1 + endif + ndata = ndata + 1 + endif + + enddo loop_report + enddo msg_report + + write(6,*)'read_fed: messages/reports = ',nmsg,'/',ntb + print*,'number of Z that is less than 0 is ppp = ', ppp + numfed=ntb + +! - Finished reading fed observations from BUFR format data file +! + call closbf(lunin) + close(lunin) + + else ! NETCDF format +!!!! Start reading fed observations from NETCDF format data file + ! CHECK IF DATA FILE EXISTS + + ! OPEN NETCDF FILE + status = nf90_open(TRIM(infile), NF90_NOWRITE, ncdfID) + print*, '*** OPENING GOES FED OBS NETCDF FILE: ', infile, status + + + !------------------------ + ! Get date information + !------------------------- + ! status = nf90_get_att( ncdfID, nf90_global, 'year', idate5s(1) ) + ! print*, 'year ',status + ! status = nf90_get_att( ncdfID, nf90_global, 'month', idate5s(2) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'day', idate5s(3) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'hour', idate5s(4) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'minute', idate5s(5) ) + ! read(idate5s(:) , *) idate5(:) + ! print*, idate5 + + !------------------------ + ! Get Dimension Info (1-D) + !------------------------- + status = nf90_inq_varid( ncdfID, 'numobs', varID ) + status = nf90_get_var( ncdfID, varID, maxobs ) + + !------------------------ + ! Allocate data arrays + !------------------------- + ALLOCATE( fed3d_column( 5, maxobs ) ) + ALLOCATE( utime( 1 ) ) ! seconds since from 2000-01-01 12:00 + + !------------------------ + ! Get useful data arrays + !------------------------- + ! LON + status = nf90_inq_varid( ncdfID, 'lon', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(1, :) ) + ! LAT + status = nf90_inq_varid( ncdfID, 'lat', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(2, :) ) + ! FED value + status = nf90_inq_varid( ncdfID, 'value', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(3, :) ) + ! TIME + status = nf90_inq_varid( ncdfID, 'time', varID ) + status = nf90_get_var( ncdfID, varID, utime ) + + ! CLOSE NETCDF FILE + status = nf90_close( ncdfID ) + + + !-Obtain analysis time in minutes since reference date + sec70 = 694267200.0 ! seconds since from 1978-01-01 00:00 to 2000-01-01 12:00 + ! because the official GOES prescribed epoch time for GLM data is 2000-01-01 12:00:00 + + call w3fs21(iadatemn,mins_an) !mins_an -integer number of mins snce 01/01/1978 + rmins_an=mins_an !convert to real number + + ! SINCE ALL OBS WILL HAVE THE SAME TIME, CHECK TIME HERE: + rmins_ob = ( utime(1) + sec70 )/60 !Convert to Minutes from seconds + twindm = twind*60. !Convert to Minutes from hours + timeb = rmins_ob-rmins_an + + if(abs(timeb) > abs(twindm)) then + print*, 'WARNING: ALL FED OBSERVATIONS OUTSIDE ASSIMILATION TIME WINDOW: ', timeb, twindm + ! goto 314 + endif + numfed = maxobs + do i=1,numfed + if (fed3d_column( 3, i ) >= fed_lowbnd2 .or. fed3d_column( 3, i ) == fed_lowbnd ) then + ndata = ndata + 1 + end if + end do + end if ! end if prebufr or netcdf format + + write(6,*)'read_fed: total no. of obs = ',ndata + nread=ndata + nodata=ndata +!!! - Finished reading fed observations from NETCDF format data file + + + + allocate(cdata_out(nreal,ndata)) +! +! + do i=1,numfed + do k=1,maxlvl + +! DCD 1 July 2021 + if (fed3d_column(k+2,i) .gt. fed_highbnd) fed3d_column(k+2,i) = fed_highbnd + + end do + end do + + write(6,*) ' ------- check max and min value of OBS: bufr fed -------' + write(6,*) ' level maxval(fed) minval(fed)' + do k=1,maxlvl + write(6,*) k,maxval(fed3d_column(k+2,:)),minval(fed3d_column(k+2,:)) + end do + + + i_maxloc=-1.0 + j_maxloc=-1.0 + k_maxloc=-1.0 + kint_maxloc=-1 + fed_max=-999.99 + ndata2=0 + do i=1,numfed + do k=1,maxlvl + if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd) then !Rong Kong + dlon_earth = fed3d_column(1,i) ! longitude (degrees) of observation + ! ilone=18 ! index of longitude (degrees) + dlat_earth = fed3d_column(2,i) ! latitude (degrees) of observation + ! ilate=19 ! index of latitude (degrees) + !-Check format of longitude and correct if necessary + if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 + if(dlon_earth=r360 .or. dlat_earth >90.0_r_kind) cycle + + !-Convert back to radians + rlon00 = dlon_earth*deg2rad + rlat00 = dlat_earth*deg2rad + call tll2xy(rlon00,rlat00,dlon,dlat,outside) + if (outside) cycle + + !If observation is outside the domain + ! then cycle, but don't increase + ! range right away. + ! Domain could be rectangular, so ob + ! may be out of + ! range at one end, but not the + ! other. + + ndata2=ndata2+1 + cdata_out( 1,ndata2) = thiserr ! obs error (flashes/min) - inflated/adjusted + + cdata_out( 2,ndata2) = dlon ! + + cdata_out( 3,ndata2) = dlat + + cdata_out( 4,ndata2) = hgt_fed(k) ! obs absolute height (m) above MSL + ! ipres=4 ! index of pressure + cdata_out( 5,ndata2) = fed3d_column(k+2,i) ! FED value + ! idbzob=5 ! index of dbz observation + cdata_out( 6,ndata2) = rstation_id ! station id (charstring equivalent to real double) + ! id=6 ! index of station id + + cdata_out( 7,ndata2) = 0.0_r_kind ! observation time in data array + ! itime=7 ! index of observation time in data array + cdata_out( 8,ndata2) = ikx ! ob type + ! ikxx=8 ! index of ob type + cdata_out( 9,ndata2) = thiserr*2.0_r_kind ! max error + ! iqmax=9 ! index of max error + cdata_out(10,ndata2) = 273.0_r_kind ! dry temperature + ! itemp=10 ! index of dry temperature + cdata_out(11,ndata2) = 1.0_r_kind ! quality mark + ! iqc=11 ! index of quality mark + cdata_out(12,ndata2) = thiserr ! original-original obs error ratio + ! ier2=12 ! index of original-original obs error ratio + cdata_out(13,ndata2) = icuse(ikx) ! index of use parameter + ! iuse=13 ! index of use parameter + cdata_out(14,ndata2) = icuse(ikx) ! dominant surface type + ! idomsfc=14 ! index of dominant surface type + cdata_out(15,ndata2) = 273.0_r_kind ! index of surface skin temperature + ! iskint=15 ! index of surface skin temperature + cdata_out(16,ndata2) = 0.5_r_kind ! 10 meter wind factor + ! iff10=16 ! index of 10 meter wind factor + cdata_out(17,ndata2) = 0.5_r_kind ! surface roughness + ! isfcr=17 ! index of surface roughness + + cdata_out(18,ndata2) = dlon_earth ! longitude (degrees) + + cdata_out(19,ndata2) = dlat_earth ! latitude (degrees) + + cdata_out(20,ndata2) = hgt_fed(k) ! station elevation (m) + ! istnelv=20 ! index of station elevation (m) + cdata_out(21,ndata2) = hgt_fed(k) ! observation height (m) + ! iobshgt=21 ! index of observation height (m) + cdata_out(22,ndata2) = hgt_fed(k) ! surface height + ! izz=22 ! index of surface height + cdata_out(23,ndata2) = fed3d_column(4,i) ! i index of obs grid for bufr resolution (i.e.,8km) + + cdata_out(24,ndata2) = fed3d_column(5,i) ! j index of obs grid for bufr resolution + + cdata_out(25,ndata2) = hgt_fed(k) ! data level category + ! icat =25 ! index of data level category + if(perturb_obs .and. fedob)then + cdata_out(26,ndata2) = 1.0_r_kind ! obs perturbation + ! iptrb=26 ! index of q perturbation + end if +! print*,'cdata_out(:,ndata2)=',cdata_out(:,ndata2) + if(fed3d_column(k+2,i) > fed_max)then + kint_maxloc=k + k_maxloc=real(k,r_kind) + j_maxloc=fed3d_column(2,i) + i_maxloc=fed3d_column(1,i) + fed_max =fed3d_column(k+2,i) + end if + endif + enddo + enddo + +!---all looping done now print diagnostic output + write(6,*)'READ_FED: Reached eof on FED file' + write(6,*)'READ_FED: # read in obs. number =',nread + write(6,*)'READ_FED: # read in obs. number for further processing =',ndata2 + ! write(6,*)'READ_FED: dlon_earth', cdata_out(18,10:15) + + ilon=2 ! array index for longitude + ilat=3 ! array index for latitude in obs information array + ndata=ndata2 + nodata=ndata2 + + !---Write observations to scratch file---! + +! if(ndata > 0 ) then + call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata + write(lunout) ((cdata_out(k,i),k=1,nreal),i=1,ndata) + ! print*,'cdata_out',cdata_out +! endif + + deallocate(cdata_out) + if (allocated(fed3d_column)) deallocate(fed3d_column) + + write(6,'(1x,A,F12.5,1x,A,3(1x,F8.3),1x,I4)') & + 'read_fed: max fed =',fed_max, '@ i j k =', & + i_maxloc,j_maxloc,k_maxloc,kint_maxloc + + end if +! close(lunout) ! ???? + return + +200 continue + write(6,*) 'read_fed, Warning : cannot find or open bufr fed data file: ', trim(infile) + +314 continue +print* ,'FINISHED WITH READ_FED' +end subroutine read_fed +! +! diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 86c7e4ce45..15476e2d04 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -192,6 +192,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) if ( .not. l_use_dbz_directDA) then if(trim(dtype) == 'dbz' )return end if + if(trim(dtype) == 'fed' )return ! Use routine as usual @@ -911,7 +912,8 @@ subroutine read_obs(ndata,mype) obstype == 'mitm' .or. obstype=='pmsl' .or. & obstype == 'howv' .or. obstype=='tcamt' .or. & obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. & - obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then + obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' .or. & + obstype=='fed') then ditype(i) = 'conv' else if (obstype == 'swcp' .or. obstype == 'lwcp') then ditype(i) = 'wcp' @@ -1302,6 +1304,10 @@ subroutine read_obs(ndata,mype) use_hgtl_full=.true. if(belong(i))use_hgtl_full_proc=.true. end if + if(obstype == 'fed')then + use_hgtl_full=.true. + if(belong(i))use_hgtl_full_proc=.true. + end if if(obstype == 'sst')then if(belong(i))use_sfc=.true. endif @@ -1639,6 +1645,12 @@ subroutine read_obs(ndata,mype) endif end if +! Process flash extent density + else if (obstype == 'fed' ) then + print *, "calling read_fed" + call read_fed(nread,npuse,nouse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i)) + string='READ_FED' + ! Process lagrangian data else if (obstype == 'lag') then call read_lag(nread,npuse,nouse,infile,lunout,obstype,& diff --git a/src/gsi/setupfed.f90 b/src/gsi/setupfed.f90 new file mode 100644 index 0000000000..cf6334e567 --- /dev/null +++ b/src/gsi/setupfed.f90 @@ -0,0 +1,1100 @@ +module fed_setup + implicit none + private + public:: setup + interface setup; module procedure setupfed; end interface + +contains +subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsave,init_pass) +!$$$ subprogram documentation block +! . . . . +! subprogram: setupfed compute rhs of flash extent density +! orig. prgmmr: +! Rong Kong CAPS/OU 2018-01-21 (modified based on setupdbz.f90) +! modified: +! Yaping Wang CIMMS/OU 2019-11-11 +! David Dowell (DCD) NOAA GSL 2021-07-01 +! - added a second option (tanh) for observation operator, based on the +! work of Sebok and Back (2021, unpublished) +! - capped maximum model FED +! +! + use mpeu_util, only: die,perr + 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 + use obsmod, only: oberror_tune + use m_obsNode, only: obsNode + use m_fedNode, only: fedNode + use m_fedNode, only: fedNode_appendto + use obsmod, only: luse_obsdiag, 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_obsLList, only: obsLList + use gsi_4dvar, only: nobs_bins,hr_obsbin + use oneobmod, only: oneobtest,maginnov,magoberr + use guess_grids, only: hrdifsig,nfldsig,ges_prsi + use guess_grids, only: ges_lnprsl, geop_hgtl + use gridmod, only: lat2, lon2 + use gridmod, only: nsig, get_ij,get_ijk,tll2xy + use constants, only: flattening,semi_major_axis,grav_ratio,zero,grav,wgtlim + use constants, only: half,one,two,grav_equator,eccentricity,somigliana + use constants, only: deg2rad,r60,tiny_r_kind,cg_term,huge_single + use constants, only: r10,r100,r1000 + use constants, only: grav,tpwcon + use qcmod, only: npres_print,ptopq,pbotq + use jfunc, only: jiter,last,miter + use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype + use convinfo, only: icsubtype + use converr, only: ptabl + use m_dtime, only: dtime_setup, dtime_check, dtime_show + use state_vectors, only: nsdim + + use gsi_bundlemod, only: GSI_BundleGetPointer + use gsi_metguess_mod, only: gsi_metguess_get, GSI_MetGuess_Bundle + + use netcdf + + + implicit none +! Declare passed variables + type(obsLList ),target,dimension(:),intent(in):: obsLL + type(obs_diags),target,dimension(:),intent(in):: odiagLL + + logical ,intent(in ) :: fed_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 + logical ,intent(in ) :: init_pass ! state of "setup" parameters + + +! Declare local parameters + integer(i_kind),parameter:: fed_obs_ob_shape = 2 ! 1 = linear (Allen et al.) + ! 2 = tanh (Sebok and Back) +! coefficients for tanh operator, from Sebok and Back (2021) +! real(r_kind),parameter:: a_coeff = 8.4_r_kind ! a (flashes/min) in tanh operator +! real(r_kind),parameter:: b_coeff = 12.248_r_kind ! b (flashes/min) in tanh operator +! real(r_kind),parameter:: c_coeff = 5.0e-10_r_kind ! c (radians/kg) in tanh operator +! real(r_kind),parameter:: d_coeff = 1.68e9_r_kind ! d (kg) in tanh operator +! real(r_kind),parameter:: fed_highbnd = 18.0_r_kind ! DCD: Sebok and Back (2021, unpublished) + +! coefficients for tanh operator, from work by A. Back with regional FV3 (2023) + real(r_kind),parameter:: a_coeff = -3.645_r_kind ! a (flashes/min) in tanh operator + real(r_kind),parameter:: b_coeff = 15.75_r_kind ! b (flashes/min) in tanh operator + real(r_kind),parameter:: c_coeff = 1.939e-10_r_kind ! c (radians/kg) in tanh operator + real(r_kind),parameter:: d_coeff = -1.215e9_r_kind ! d (kg) in tanh operator + real(r_kind),parameter:: fed_highbnd = 8.0_r_kind ! DCD: Back (2023, unpublished) for FV3 + + real(r_kind),parameter:: fed_height = 6500.0_r_kind ! assumed height (m) of FED observations + real(r_kind),parameter:: r0_001 = 0.001_r_kind + real(r_kind),parameter:: r8 = 8.0_r_kind + real(r_kind),parameter:: ten = 10.0_r_kind + real(r_kind),parameter:: D608=0.608_r_kind + character(len=*),parameter:: myname='setupfed' + +! Declare external calls for code analysis + external:: tintrp2a1 + external:: tintrp2a11 + external:: tintrp2a1116 + external:: tintrp31 + external:: grdcrd1 + external:: stop2 + +! Declare local variables + real(r_kind) rlow,rhgh,rsig + real(r_kind) dz + real(r_kind) jqg_num,jqg + real(r_kind) dlnp,pobl,zob + real(r_kind) sin2,termg,termr,termrg + real(r_kind) psges,zsges + real(r_kind),dimension(nsig):: zges,hges + real(r_kind) prsltmp(nsig) + real(r_kind) sfcchk + real(r_kind) residual,obserrlm,obserror,ratio,scale,val2 + real(r_kind) ress,ressw + real(r_kind) val,valqc,rwgt + real(r_kind) cg_w,wgross,wnotgross,wgt,arg,exp_arg,term,rat_err2 + real(r_double) rstation_id + real(r_kind) dlat,dlon,dtime,dpres,ddiff,error,slat,dlat8km,dlon8km + real(r_kind) ratio_errors + real(r_kind) presw + 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(lat2,lon2,nfldsig)::rp + real(r_single),allocatable,dimension(:,:)::rdiagbuf + real(r_kind),allocatable,dimension(:,:,: ) :: ges_z + real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_qg,ges_qg_mask + + real(r_kind) :: presq + real(r_kind) :: T1D,RHO + real(r_kind) :: glmcoeff = 2.088_r_kind*10.0**(-8.0) ! Allen et al. (2016,MWR) + real(r_kind) :: CM = 0.5_r_kind ! tuning factor in eq. 14 of Kong et al. 2020 + + integer(i_kind) i,nchar,nreal,k,j,k1,ii,jj + integer(i_kind) mm1,k2 + integer(i_kind) jsig,ikxx,nn,ibin,ioff,ioff0 + integer(i_kind) ier,ilat,ilon,ifedob,ikx,itime,iuse + integer(i_kind) id,ilone,ilate + integer(i_kind) ier2 + + integer(i_kind) nlat_ll,nlon_ll,nsig_ll,nfld_ll + + integer(i_kind) ipres,iqmax,iqc,icat,itemp + integer(i_kind) istnelv,iobshgt,izz,iprvd,isprvd,iptrb + integer(i_kind) idomsfc,iskint,isfcr,iff10 + + character(8) station_id + character(8),allocatable,dimension(:):: cdiagbuf + character(80):: string + character(128):: diag_file + logical :: diagexist + logical,dimension(nobs):: luse,muse + integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID + logical proceed + + equivalence(rstation_id,station_id) + integer(i_kind) numequal,numnotequal + + type(fedNode),pointer:: my_head + type(obs_diag),pointer:: my_diag + type(obs_diags),pointer:: my_diagLL + + real(r_kind),dimension(nsig+1):: prsitmp + + +!------------------------------------------------! + + integer(i_kind) :: itmp,jtmp + + integer(i_kind), parameter :: ntimesfed=1 + integer(i_kind),parameter :: nxfed=99, nyfed=99, nzfed=1, nfldfed=3 + real(r_kind),dimension(nobs) :: FEDMdiag,FEDMdiagTL + integer(i_kind) :: npt + real(r_kind) :: dlat_earth,dlon_earth + +! YPW added the next lines + logical :: l_set_oerr_ratio_fed=.False. + logical :: l_gpht2gmht = .True. + real(r_kind),dimension(nobs) :: dlatobs,dlonobs + integer(i_kind):: ngx,ngy,igx,jgy + real(r_kind):: dx_m, dy_m + + type(obsLList),pointer,dimension(:):: fedhead + fedhead => obsLL(:) + +!============================================================================================ +! Read and reformat observations in work arrays. + read(lunin)data,luse,ioid + + write(6,*)myname,'(pe=',mype,') nele nobs =',nele,nobs, & + 'luse_obsdiag=',luse_obsdiag,'lat2,lon2=',lat2,lon2 + + 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 + ifedob=5 ! index of fed observation + id=6 ! index of station id + itime=7 ! index of observation time in data array + ikxx=8 ! index of ob type + iqmax=9 ! index of max error + itemp=10 ! index of dry temperature + 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) + iobshgt=21 ! index of observation height (m) + izz=22 ! index of surface height + iprvd=23 ! index of observation provider + isprvd=24 ! index of observation subprovider + icat =25 ! index of data level category + iptrb=26 ! index of fed perturbation + do i=1,nobs + muse(i)=nint(data(iuse,i)) <= jiter + end do + + numequal=0 + numnotequal=0 + +! +! If requested, save select data for output to diagnostic file + if(fed_diagsave)then + ii=0 + nchar=1_i_kind + ioff0=26_i_kind ! 21 + 5 (22->Zr; 23->Zs; 24->Zg; 25->tsenges;26->RHO;) + nreal=ioff0 + if (lobsdiagsave) nreal=nreal+4*miter+1 + allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs)) + rdiagbuf=zero + if(netcdf_diag) call init_netcdf_diag_ + end if + mm1=mype+1 + scale=one + rsig=nsig + + + + !============================================================================================ +! +! Check to see if required guess fields are available +! vars. list: ps, z, q +! vars. list: qr, qs, qg + !============================================================================================ + + call check_vars_(proceed) + if(.not.proceed) then + write(6,*) myname,': some or all necessary variables are not available for fed obs operator. Quit!' + return ! not all vars available, simply return + end if + +! If require guess vars available, extract from bundle ... + call init_vars_ +! qscalar=zero + + !============================================================================================ + ! 1) Calculate the graupel-mass and graupel-volume based flash extent density + ! (FED) on model space, added by R. Kong, 07/05/2018 + !============================================================================================ + ges_qg_mask=ges_qg + where(ges_qg>0.0005) !Count the volume where qg > 0.5/kg + ges_qg_mask=1.0 + elsewhere + ges_qg_mask=0.0 + endwhere + + ! Operator start here + ! set ngx and ngy =2, so the integrated domain is 15kmx15km + ngx = 2 + ngy = 2 + dx_m = 3000. + dy_m = 3000. + print*,'FED Operator start here!,ngx=',ngx,'ngy=',ngy + rp=zero + + print*, 'mype = ', mype + print*, 'nfldsig = ', nfldsig + print*, 'nsig = ', nsig + print*, 'lon2 = ', lon2 + print*, 'lat2 = ', lat2 + +! compute graupel mass, in kg per 15 km x 15 km column + do jj=1,nfldsig + do k=1,nsig + do i=1,lon2 + do j=1,lat2 !How to handle MPI???? + do igx=1,2*ngx+1 !horizontal integration of qg around point to calculate FED + do jgy=1,2*ngy+1 + itmp = i-ngx+igx-1 + jtmp = j-ngy+jgy-1 + itmp = min(max(1,itmp),lon2) + jtmp = min(max(1,jtmp),lat2) + rp(j,i,jj)=rp(j,i,jj) + ges_qg(jtmp,itmp,k,jj)* & + dx_m*dy_m*(ges_prsi(jtmp,itmp,k,jj)-ges_prsi(jtmp,itmp,k+1,jj))*& + tpwcon * r10 + end do !igx + end do !jgy + end do !j + end do !i + end do !k + end do !jj + +! compute FED, in flashes/min + do jj=1,nfldsig + do i=1,lon2 + do j=1,lat2 + if (fed_obs_ob_shape .eq. 1) then + rp(j,i,jj) = CM * glmcoeff * rp(j,i,jj) + else if (fed_obs_ob_shape .eq. 2) then + rp(j,i,jj) = a_coeff + b_coeff & + * tanh(c_coeff * (rp(j,i,jj) - d_coeff)) + else + write(6,*) ' unknown fed_obs_ob_shape: ', fed_obs_ob_shape + write(6,*) ' aborting setupfed' + call stop2(999) + end if + if (rp(j,i,jj) .gt. fed_highbnd) rp(j,i,jj) = fed_highbnd + end do !j + end do !i + end do !jj + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + write(6,*) 'fed_obs_ob_shape=',fed_obs_ob_shape + if (fed_obs_ob_shape .eq. 2) then + write(6,*) 'a_coeff=',a_coeff + write(6,*) 'b_coeff=',b_coeff + write(6,*) 'c_coeff=',c_coeff + write(6,*) 'd_coeff=',d_coeff + end if + write(6,*) 'fed_highbnd=',fed_highbnd + write(6,*) 'maxval(ges_qg)=',maxval(ges_qg),'pe=',mype + + + !============================================================================================ + + nlat_ll=size(ges_qg,1) + nlon_ll=size(ges_qg,2) + nsig_ll=size(ges_qg,3) + nfld_ll=size(ges_qg,4) + + +! - Observation times are checked in read routine - comment out for now + +! call dtime_setup() + +!print*,"maxval(data(ifedob,:)),mmaxval(data(ilat,:))=",minval(data(ifedob,:)),maxval(data(ifedob,:)),maxval(data(ilat,:)) +!write(6,*) "OKOKOKOKOK, nobs=", nobs + do i=1,nobs + dtime=data(itime,i) + dlat=data(ilat,i) + dlon=data(ilon,i) + + dlon8km=data(iprvd,i) !iprvd=23 + dlat8km=data(isprvd,i) !isprvd=24 + + dpres=data(ipres,i) ! from rdararef_mosaic2: this height abv MSL + ikx = nint(data(ikxx,i)) + error=data(ier2,i) + slat=data(ilate,i)*deg2rad ! needed when converting geophgt to + dlon_earth = data(ilone,i) !the lontitude and latitude on the obs pts. + dlat_earth = data(ilate,i) + ! geometric hgh (hges --> zges below) + + if (nobs_bins>1) then + ibin = NINT( dtime/hr_obsbin ) + 1 + else + ibin = 1 + end if + + 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) + end if + +! Interpolate terrain height(model elevation) to obs location. + call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) +! print*,'i,after tintrp2all',i,mype,dlat,zsges +! 1. dpres (MRMS obs height is height above MSL) is adjusted by zsges, so it +! is changed to height relative to model elevation (terrain). +! because in GSI, geop_hgtl is the height relative to terrain (ges_z) +! (subroutine guess_grids) + dpres=dpres-zsges + if (dpres rsig)ratio_errors = zero + +!----------------------------------------------------------------------------! +! ! +! Implementation of forward operator for flash extend densit ----------------! +! ! +!----------------------------------------------------------------------------! + + !============================================================================================ + ! 3) H(x), interpolate the FED from model space on the local domain to obs space (FEDMdiag) + !============================================================================================ + + npt = 0 + FEDMdiag(i) = 0. + call tintrp2a11(rp,FEDMdiag(i),dlat,dlon,dtime,hrdifsig,mype,nfldsig) + dlonobs(i) = dlon_earth + dlatobs(i) = dlat_earth + + ! also Jacobian used for TLM and ADM + !FEDMdiagTL, used for gsi-3dvar,will be implemented in future...... + FEDMdiagTL(i) = 0. + jqg_num = FEDMdiagTL(i) !=dFED/Dqg + jqg = jqg_num + + + !end select + + if(FEDMdiag(i)==data(ifedob,i)) then + numequal=numequal+1 + else + numnotequal=numnotequal+1 + end if + +!!!!!!!!!!!!!!!!!END H(x)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Compute innovations + !--------------Calculate departure from observation----------------! + + ddiff = data(ifedob,i) - FEDMdiag(i) + +! If requested, setup for single obs test. +! Note: do not use this way to run single obs test for fed in the current version. (g.zhao) + if (oneobtest) then + ddiff=maginnov +! if (trim(adjustl(oneob_type))=='fed') then +! data(ifedob,i) = maginnov +! ddiff = data(ifedob,i) - FEDMdiag(i) +! end if + error=one/(magoberr) + ratio_errors=one + end if + + +! Gross error checks + obserror = one/max(ratio_errors*error,tiny_r_kind) + obserrlm = max(cermin(ikx),min(cermax(ikx),obserror)) + + residual = abs(ddiff) != y-H(xb) + ratio = residual/obserrlm != y-H(xb)/sqrt(R) + + if (l_set_oerr_ratio_fed) then + if (ratio > cgross(ikx) .or. ratio_errors < tiny_r_kind) then + if (luse(i)) awork(4) = awork(4)+one + error = zero + ratio_errors = zero + end if + else + ratio_errors = one + end if + + if (ratio_errors*error <=tiny_r_kind) muse(i)=.false. +! if (nobskeep>0 .and. luse_obsdiag) muse(i)=obsdiags(i_fed_ob_type,ibin)%tail%muse(nobskeep) + if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i)) + + val = error*ddiff !=y-H(xb)/sqrt(R) + +! Compute penalty terms (linear & nonlinear qc). + if(luse(i))then + exp_arg = -half*val**2 + rat_err2 = ratio_errors**2 + val2=val*val !(o-g)**2/R, would be saved in awork + if (cvar_pg(ikx) > tiny_r_kind .and. error > tiny_r_kind) then + arg = exp(exp_arg) + wnotgross= one-cvar_pg(ikx) + cg_w=cvar_b(ikx) + wgross = cg_term*cvar_pg(ikx)/(cg_w*wnotgross) + term = log((arg+wgross)/(one+wgross)) + wgt = one-wgross/(arg+wgross) + rwgt = wgt/wgtlim + else + term = exp_arg + wgt = wgtlim + rwgt = wgt/wgtlim + end if + valqc = -two*rat_err2*term + +! print*,'Compute penalty terms' +! Accumulate statistics for obs belonging to this task + if (muse(i)) then + if(rwgt < one) awork(21) = awork(21)+one + jsig = dpres + jsig=max(1,min(jsig,nsig)) + awork(6*nsig+jsig+100)=awork(6*nsig+jsig+100)+val2*rat_err2 + awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+one + awork(3*nsig+jsig+100)=awork(3*nsig+jsig+100)+valqc + end if +! Loop over pressure level groupings and obs to accumulate +! statistics as a function of observation type. + 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 + do k = 1,npres_print +! if(presw >=ptop(k) .and. presw<=pbot(k))then + if(presq >=ptopq(k) .and. presq<=pbotq(k))then + 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 + end if + + if(luse_obsdiag)then + call obsdiagNode_set(my_diag, luse=luse(i), wgtjo=(error*ratio_errors)**2, & + jiter=jiter, muse=muse(i), nldepart=ddiff) + end if + +! 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) ! YPW added + call fedNode_appendto(my_head,fedhead(ibin)) + + my_head%idv=is + my_head%iob=i + my_head%elat= data(ilate,i) + my_head%elon= data(ilone,i) + + my_head%dlev= dpres + call get_ijk(mm1,dlat,dlon,dpres,my_head%ij,my_head%wij) + my_head%res = ddiff ! Observation - ges + my_head%err2 = error**2 + my_head%raterr2 = ratio_errors**2 + my_head%time = dtime + my_head%b = cvar_b(ikx) + my_head%pg = cvar_pg(ikx) + my_head%luse = luse(i) + + if(oberror_tune) then + ! my_head%fedpertb=data(iptrb,i)/error/ratio_errors + my_head%kx=ikx + if(presq > ptabl(2))then + my_head%k1=1 + else if( presq <= ptabl(33)) then + my_head%k1=33 + else + k_loop: do k=2,32 + if(presq > ptabl(k+1) .and. presq <= ptabl(k)) then + my_head%k1=k + exit k_loop + end if + end do k_loop + end if + end if +!------------------------------------------------- + + 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 + end if + + my_head => null() + end if + +! Save select output for diagnostic file + if(.not.luse(i))write(6,*)' luse, mype',luse(i),mype + if(fed_diagsave .and. luse(i))then + ii=ii+1 + 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 + end if + 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 do + +! Release memory of local guess arrays + call final_vars_ + +! Write information to diagnostic file + if(fed_diagsave .and. netcdf_diag) call nc_diag_write + if(fed_diagsave .and. binary_diag .and. ii>0)then + + write(string,600) jiter +600 format('fed_',i2.2) + diag_file=trim(dirname) // trim(string) + if(init_pass) then + open(66,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + else + inquire(file=trim(diag_file),exist=diagexist) + if (diagexist) then + open(66,file=trim(diag_file),form='unformatted',status='old',position='append') + else + open(66,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + end if + end if + if(init_pass .and. mype == 0) then + write(66) ianldate + write(6,*)'SETUPFED: write time record to file ',& + trim(diag_file), ' ',ianldate + end if + +! call dtime_show(myname,'diagsave:fed',i_fed_ob_type) + write(66)'fed',nchar,nreal,ii,mype,ioff0 + write(66)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + ! write(6,*)'fed,nchar,nreal,ii,mype',nchar,nreal,ii,mype + deallocate(cdiagbuf,rdiagbuf) + close(66) + 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::q' , ivar, istatus ) + proceed=proceed.and.ivar>0 +! call gsi_metguess_get ('var::tv' , ivar, istatus ) +! proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::qs', ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::qg', ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::qr', ivar, istatus ) + proceed=proceed.and.ivar>0 + end subroutine check_vars_ + + + subroutine init_vars_ + +! use radaremul_cst, only: mphyopt + + 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) + end if + 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 + end do + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + end if +! 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) + end if + 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 + end do + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + end if +! get q ... + varname='q' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_q))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + end if + allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_q(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_q(:,:,:,ifld)=rank3 + end do + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + end if +! 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) +! end if +! 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 +! ges_tv(:,:,:,ifld)=rank3 +! end do +! else +! write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus +! call stop2(999) +! end if +! get qr ... +! get qg ... + varname='qg' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_qg))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + end if + allocate(ges_qg(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + if(.not. allocated(ges_qg_mask))then + allocate(ges_qg_mask(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + end if + + ges_qg(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_qg(:,:,:,ifld)=rank3 + end do + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + end if + else + write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& + nfldsig,size(gsi_metguess_bundle) + call stop2(999) + end if + 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_fed_',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 + end if + end if + + call nc_diag_init(diag_conv_file, append=append_diag) + + if (.not. append_diag) then ! don't write headers on append - the module will break? + call nc_diag_header("date_time",ianldate ) + call nc_diag_header("Number_of_state_vars", nsdim ) + end if + 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) = presq ! observation pressure (hPa) + rdiagbuf(7,ii) = fed_height ! observation height (meters) + rdiagbuf(8,ii) = (dtime*r60)-time_offset ! obs time (sec 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 + end if + + rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight + rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (dBZ)**-1 + rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (dBZ)**-1 + rdiagbuf(16,ii) = errinv_final ! final inverse observation error (dBZ)**-1 + rdiagbuf(17,ii) = data(ifedob,i) ! radar reflectivity observation (dBZ) + rdiagbuf(18,ii) = ddiff ! obs-ges (dBZ) + rdiagbuf(19,ii) = data(ifedob,i)-FEDMdiag(i) ! obs-ges w/o bias correction (dBZ) (future slot) + rdiagbuf(20,ii) = dlat8km ! j-index on 8km bufr obs grid + rdiagbuf(21,ii) = dlon8km ! i-index on 8km bufr obs grid + +! print*,'data(ilat,i)=',data(ilat,i),'data(ilon,i)=',data(ilon,i) + + rdiagbuf(22,ii) = FEDMdiag(i) ! dBZ from rain water + + rdiagbuf(23,ii) = T1D ! temperature (sensible, K) + rdiagbuf(24,ii) = RHO ! air density (kg/m**3) + + if (lobsdiagsave) then + write(6,*)'wrong here, stop in setupfed.f90 ' + stop + ioff=nreal + do jj=1,miter + ioff=ioff+1 + if (odiag%muse(jj)) then + rdiagbuf(ioff,ii) = one + else + rdiagbuf(ioff,ii) = -one + end if + end do + do jj=1,miter+1 + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%nldepart(jj) + end do + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%tldepart(jj) + end do + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%obssen(jj) + end do + end if + + end subroutine contents_binary_diag_ + subroutine contents_netcdf_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag +! Observation class + character(7),parameter :: obsclass = ' fed' + 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(presq) ) + call nc_diag_metadata("Height", sngl(fed_height) ) + call nc_diag_metadata("Time", sngl(dtime*r60-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) ) + end if + + 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(data(ifedob,i)) ) + call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ddiff) ) + call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(data(ifedob,i)-FEDMdiag(i)) ) + if (lobsdiagsave) then + do jj=1,miter + if (odiag%muse(jj)) then + obsdiag_iuse(jj) = one + else + obsdiag_iuse(jj) = -one + end if + end do + + 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 ) + end if + + end subroutine contents_netcdf_diag_ + + subroutine final_vars_ + if(allocated(ges_z )) deallocate(ges_z ) + if(allocated(ges_q )) deallocate(ges_q ) +! if(allocated(ges_tv)) deallocate(ges_tv) + if(allocated(ges_ps)) deallocate(ges_ps) + if(allocated(ges_qg)) deallocate(ges_qg) + end subroutine final_vars_ + + subroutine init_qcld(t_cld, qxmin_cld, icat_cld, t_dpnd) + use kinds, only: r_kind,r_single,r_double,i_kind + implicit none + real(r_kind), intent(in ) :: t_cld + real(r_kind), intent(inout) :: qxmin_cld + integer, intent(in ) :: icat_cld + logical, intent(in ) :: t_dpnd +! +! local variables + real :: tr_ll, qrmin_ll, tr_hl, qrmin_hl + real :: ts_ll, qsmin_ll, ts_hl, qsmin_hl + real :: tg_ll, qgmin_ll, tg_hl, qgmin_hl + real :: qr_min, qs_min, qg_min +!------------------------------------------------------ + + qr_min = 5.0E-6_r_kind + qs_min = 5.0E-6_r_kind + qg_min = 5.0E-6_r_kind + tr_ll = 275.65; qrmin_ll = 5.0E-6_r_kind; + tr_hl = 270.65; qrmin_hl = 1.0E-8_r_kind; + ts_ll = 275.65; qsmin_ll = 1.0E-8_r_kind; + ts_hl = 270.65; qsmin_hl = 5.0E-6_r_kind; + tg_ll = 275.65; qgmin_ll = 1.0E-6_r_kind; + tg_hl = 270.65; qgmin_hl = 5.0E-6_r_kind; + + select case (icat_cld) + case (1) + if ( t_dpnd ) then + if (t_cld <= tr_hl) then + qxmin_cld = qrmin_hl + else if (t_cld >= tr_ll) then + qxmin_cld = qrmin_ll + else + qxmin_cld = (qrmin_hl + qrmin_ll) * 0.5 + end if + else + qxmin_cld = qr_min + end if + case default + write(6,*) 'wrong cloud hydrometer category ID',icat_cld + end select + + return + + end subroutine init_qcld + +end subroutine setupfed +end module fed_setup diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 3efcb69859..8075956431 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -168,6 +168,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) i_gust,i_vis,i_pblh,i_wspd10m,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 use m_rhs, only: i_light use m_gpsStats, only: gpsStats_genstats ! was genstats_gps() @@ -625,7 +626,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) 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_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_dbz, & + i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz, & size(awork1,2),bwork1,awork1,ndata) ! Compute and print statistics for "lightning" data diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90 index a01675d8d0..0da8606f24 100644 --- a/src/gsi/statsconv.f90 +++ b/src/gsi/statsconv.f90 @@ -2,7 +2,7 @@ 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_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,& - i_swcp,i_lwcp,i_dbz,i_ref,bwork,awork,ndata) + i_swcp,i_lwcp,i_fed,i_dbz,i_ref,bwork,awork,ndata) !$$$ subprogram documentation block ! . . . . ! subprogram: statconv prints statistics for conventional data @@ -74,6 +74,7 @@ subroutine statsconv(mype,& ! i_vwnd10m- index in awork array holding vwnd10m info ! i_swcp - index in awork array holding swcp info ! i_lwcp - index in awork array holding lwcp info +! i_fed - index in awork array holding fed info ! i_dbz - index in awork array holding dbz info ! i_ref - size of second dimension of awork array ! bwork - array containing information for statistics @@ -96,12 +97,12 @@ subroutine statsconv(mype,& iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_td2m,& iout_mxtm,iout_mitm,iout_pmsl,iout_howv,iout_tcamt,iout_lcbas,iout_cldch,& iout_uwnd10m,iout_vwnd10m,& - iout_dbz,iout_swcp,iout_lwcp,& + 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_pmsl,mype_howv,mype_tcamt,mype_lcbas,mype_cldch,mype_uwnd10m,mype_vwnd10m,& - mype_dbz,mype_swcp,mype_lwcp + mype_fed,mype_dbz,mype_swcp,mype_lwcp use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq use jfunc, only: first,jiter use gridmod, only: nsig @@ -112,7 +113,7 @@ subroutine statsconv(mype,& 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_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_dbz,i_ref + 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 integer(i_kind),dimension(ndat,3) ,intent(in ) :: ndata @@ -136,6 +137,7 @@ subroutine statsconv(mype,& real(r_kind) dwqcplty,tqcplty,qctt,qctrw,rwqcplty,qctdw,qqcplty,qctgps real(r_kind) gpsqcplty,tpw3,pw3,qctq real(r_kind) tswcp3,tlwcp3,qctdbz,dbzqcplty + real(r_kind) fedmplty,tfed,qctfed,fedqcplty real(r_kind),dimension(1):: pbotall,ptopall logical,dimension(nconvtype):: pflag @@ -1325,6 +1327,68 @@ subroutine statsconv(mype,& end if end if +! Summary report for flash extent density + if(mype==mype_fed) then + nread=0 + nkeep=0 + do i=1,ndat + if(dtype(i)== 'fed')then + nread=nread+ndata(i,2) + nkeep=nkeep+ndata(i,3) + end if + end do + if(nread > 0)then + if(first)then + open(iout_fed) + else + open(iout_fed,position='append') + end if + + fedmplty=zero; fedqcplty=zero ; ntot=0 + tfed=zero ; qctfed=zero + if(nkeep > 0)then + mesage='current vfit of flash extent density, ranges in flashes per minute$' + do j=1,nconvtype + pflag(j)=trim(ioctype(j)) == 'fed' + end do + call dtast(bwork,npres_print,pbot,ptop,mesage,jiter,iout_fed,pflag) + + numgross=nint(awork(4,i_fed)) + numfailqc=nint(awork(21,i_fed)) + do k=1,nsig + num(k)=nint(awork(k+5*nsig+100,i_fed)) + rat=zero + rat3=zero + if(num(k) > 0) then + rat=awork(6*nsig+k+100,i_fed)/float(num(k)) + rat3=awork(3*nsig+k+100,i_fed)/float(num(k)) + end if + ntot=ntot+num(k) + fedmplty=fedmplty+awork(6*nsig+k+100,i_fed) + fedqcplty=fedqcplty+awork(3*nsig+k+100,i_fed) + write(iout_fed,240) 'r',num(k),k,awork(6*nsig+k+100,i_fed), & + awork(3*nsig+k+100,i_fed),rat,rat3 + end do + if(ntot > 0) then + tfed=fedmplty/float(ntot) + qctfed=fedqcplty/float(ntot) + end if + write(iout_fed,925) 'fed',numgross,numfailqc + numlow = nint(awork(2,i_fed)) + numhgh = nint(awork(3,i_fed)) + nhitopo = nint(awork(5,i_fed)) + ntoodif = nint(awork(6,i_fed)) + write(iout_fed,900) 'fed',numhgh,numlow + write(iout_fed,905) 'fed',nhitopo,ntoodif + end if + write(iout_fed,950) 'fed',jiter,nread,nkeep,ntot + write(iout_fed,951) 'fed',fedmplty,fedqcplty,tfed,qctfed + + close(iout_fed) + end if + end if + + if(mype==mype_tcp) then nread=0 nkeep=0