Skip to content

Commit

Permalink
Gsi fed (#590)
Browse files Browse the repository at this point in the history
<!-- PLEASE READ -->
<!--
Before opening a PR, please note these guidelines:

- Each PR should only address ONE topic and have an associated issue
- No hardcoded or paths to personal directories should be present
- No temporary or backup files should be committed
- Any code that was disabled by being commented out should be removed
-->

**Description**

Initialization of the operational RRFSv1 will include assimilation of
flash-extent density (FED) observations from the GOES Geostationary
Lightning Mapper (GLM). The current PR is the first of at least 3 that
will be needed to introduce the capability of FED assimilation into the
code and regional workflow. The new capabilities that are added to GSI
are:

* reading NetCDF FED observations
* applying an observation operator that maps the model state to FED.

Much of the code was originally developed by Rong Kong at OU-CAPS (Kong
et al. 2020, Wang et al. 2021, Kong et al. 2022;
https://doi.org/10.1175/MWR-D-19-0192.1,
https://doi.org/10.1175/MWR-D-20-0406.1,
https://doi.org/10.1175/MWR-D-21-0326.1). Recently, the observation
operator has been modified by Amanda Back and Ashley Sebok based on
tests with regional, convection-allowing FV3 forecasts. The new
observation operator includes a cap of 8 flashes / minute for both the
observed and simulated FED.

The observation operator is specific to the 3-km regional FV3
application in RRFS. Development of a more general observation operator
is left to future work.

Fixes #588 

<!-- Please include relevant motivation and context. -->
<!-- Please include a summary of the change and which issue is fixed.
-->
<!-- List any dependencies that are required for this change. -->

<!-- Please provide reference to the issue this pull request is
addressing. -->
<!-- For e.g. Fixes #IssueNumber -->

**Type of change**

Please delete options that are not relevant.

- [ ] Bug fix (non-breaking change which fixes an issue)
- [X] New feature (non-breaking change which adds functionality)
- [ ] Breaking change (fix or feature that would cause existing
functionality to not work as expected)
- [ ] This change requires a documentation update

**How Has This Been Tested?**

Initial tests were with NOAA-EMC GSI-EnKF code obtained in April 2023
and modified to include the assimilation of FED observations. A
prototype of RRFSv1 was cycled hourly for 2.5 days, and the EnKF
assimilation included FED data assimilation.

For the current PR, only the GSI observer with FED (and radar
reflectivity) observations was tested. It produces identical results to
those obtained in April 2023.

<!-- Please describe the tests that you ran to verify your changes and
on the platforms these tests were conducted. -->
<!-- Provide instructions so we can reproduce. -->
<!-- Please also list any relevant details for your test configuration
-->
  
**Checklist**

- [ ] My code follows the style guidelines of this project
- [X] I have performed a self-review of my own code
- [ ] I have commented my code, particularly in hard-to-understand areas
- [ ] New and existing tests pass with my changes
- [ ] Any dependent changes have been merged and published

**DUE DATE for this PR is 8/24/2023.** If this PR is not merged into
`develop` by this date, the PR will be closed and returned to the
developer.

---------

Co-authored-by: Ming Hu <Ming.Hu@noaa.gov>
  • Loading branch information
daviddowellNOAA and hu5970 authored Sep 8, 2023
1 parent be4a3d9 commit d7ac706
Show file tree
Hide file tree
Showing 14 changed files with 2,169 additions and 13 deletions.
174 changes: 174 additions & 0 deletions src/gsi/gsi_fedOper.F90
Original file line number Diff line number Diff line change
@@ -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
4 changes: 4 additions & 0 deletions src/gsi/gsi_files.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -532,6 +535,7 @@ setupco.f90
setupdbz.f90
setupdbz_lib.f90
setupdw.f90
setupfed.f90
setupgust.f90
setuphowv.f90
setuplag.f90
Expand Down
7 changes: 7 additions & 0 deletions src/gsi/gsi_obOperTypeManager.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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_
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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.
Expand Down
7 changes: 5 additions & 2 deletions src/gsi/gsimod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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,&
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/gsi/intjo.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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, &
Expand Down
Loading

0 comments on commit d7ac706

Please sign in to comment.