Add lightning data assimilation (#1962)

TYPE: new feature

KEYWORDS: lightning data, vertical velocity, pseudo divergence, pseudo water vapor

SOURCE: Zhixiong Chen (zchen@fjnu.edu.cn, Fujian Normal University, China)

DESCRIPTION OF CHANGES:
A new lightning data assimilation scheme has been implemented in the WRFDA. In this lightning DA scheme, the lightning flash rate is first converted to the maximum vertical velocity using an empirical relationship and then expanded to profiles of vertical velocities. The profiles of vertical velocity can be directly assimilated with the control variable of W. Another way is to assimilate pseudo divergence converted from vertical velocity to adjust the Kinematic states. A scheme to assimilate pseudo water vapor retrievals from lightning flash rate observations is also added. In the humidity assimilation scheme, the mixing ratio profiles of pseudo water vapor are retrieved from the locations where a rapid increase in lightning flash rate is detected. 

LIST OF MODIFIED FILES: 
M       Registry/registry.var
M       frame/module_driver_constants.F
M       var/build/da.make
M       var/build/depend.txt
M       var/da/da_control/da_control.f90
M       var/da/da_define_structures/da_allocate_observations.inc
M       var/da/da_define_structures/da_allocate_y.inc
A       var/da/da_define_structures/da_allocate_y_lightning.inc
M       var/da/da_define_structures/da_deallocate_observations.inc
M       var/da/da_define_structures/da_deallocate_y.inc
M       var/da/da_define_structures/da_define_structures.f90
M       var/da/da_define_structures/da_zero_y.inc
A       var/da/da_lightning/da_ao_stats_lightning.inc
A       var/da/da_lightning/da_calculate_grady_lightning.inc
A       var/da/da_lightning/da_check_max_iv_lightning.inc
A       var/da/da_lightning/da_div_profile.inc
A       var/da/da_lightning/da_div_profile_adj.inc
A       var/da/da_lightning/da_div_profile_tl.inc
A       var/da/da_lightning/da_get_innov_vector_lightning.inc
A       var/da/da_lightning/da_jo_and_grady_lightning.inc
A       var/da/da_lightning/da_lightning.f90
A       var/da/da_lightning/da_oi_stats_lightning.inc
A       var/da/da_lightning/da_print_stats_lightning.inc
A       var/da/da_lightning/da_residual_lightning.inc
A       var/da/da_lightning/da_transform_xtoy_lightning.inc
A       var/da/da_lightning/da_transform_xtoy_lightning_adj.inc
M       var/da/da_main/da_wrfvar_top.f90
M       var/da/da_minimisation/da_calculate_grady.inc
M       var/da/da_minimisation/da_calculate_residual.inc
M       var/da/da_minimisation/da_get_innov_vector.inc
M       var/da/da_minimisation/da_get_var_diagnostics.inc
M       var/da/da_minimisation/da_jo_and_grady.inc
M       var/da/da_minimisation/da_minimisation.f90
M       var/da/da_minimisation/da_write_diagnostics.inc
M       var/da/da_obs/da_fill_obs_structures.inc
A       var/da/da_obs/da_fill_obs_structures_lightning.inc
M       var/da/da_obs/da_obs.f90
M       var/da/da_obs/da_obs_sensitivity.inc
M       var/da/da_obs/da_transform_xtoy.inc
M       var/da/da_obs/da_transform_xtoy_adj.inc
M       var/da/da_obs_io/da_final_write_obs.inc
M       var/da/da_obs_io/da_obs_io.f90
M       var/da/da_obs_io/da_read_iv_for_multi_inc.inc
A       var/da/da_obs_io/da_read_obs_lightning.inc
M       var/da/da_obs_io/da_read_omb_tmp.inc
A       var/da/da_obs_io/da_scan_obs_lightning.inc
M       var/da/da_obs_io/da_search_obs.inc
M       var/da/da_obs_io/da_write_iv_for_multi_inc.inc
M       var/da/da_obs_io/da_write_obs.inc
M       var/da/da_setup_structures/da_setup_obs_structures.inc
M       var/da/da_setup_structures/da_setup_obs_structures_ascii.inc
A       var/da/da_setup_structures/da_setup_obs_structures_lightning.inc
M       var/da/da_setup_structures/da_setup_structures.f90
M       var/da/da_statistics/da_analysis_stats.inc
M       var/da/da_statistics/da_statistics.f90
M       var/da/da_test/da_check_xtoy_adjoint.inc
A       var/da/da_test/da_check_xtoy_adjoint_lightning.inc
M       var/da/da_test/da_get_y_lhs_value.inc
M       var/da/da_test/da_test.f90

TESTS CONDUCTED: 
1. Tested it with 3DVar on Cheyenne using the intel compiler.

RELEASE NOTE: 
A lightning data assimilation scheme has been added to assimilate pseudo vertical velocity, divergence fields, or water vapor retrievals from lightning flash rate data.
Chen, Z., Sun, J., Qie, X., Zhang, Y., Ying, Z., Xiao, X., & Cao, D. (2020). A method to update model kinematic states by assimilating satellite-observed total lightning data to improve convective analysis and forecasting. Journal of Geophysical Research: Atmospheres, 125, e2020JD033330. https://doi.org/10.1029/2020JD033330

Registry/registry.var

@@ -177,6 +177,12 @@ rconfig   real      rfmin                   namelist,wrfvar4  1  0.0      - "min
 rconfig   integer   rf_noice                namelist,wrfvar4  1  0        - "disable ice phace in H"   ""  ""
 rconfig   real      radar_rf_rscl           namelist,wrfvar4  1  1.0      - "weight of rf"             ""  ""
 rconfig   real      radar_rv_rscl           namelist,wrfvar4  1  1.0      - "weight of rv"             ""  ""
+rconfig   logical   use_lightningobs        namelist,wrfvar4  1  .false.  - "use_lightningobs"         ""  ""
+rconfig   logical   use_lightning_w         namelist,wrfvar4  1  .false.  - "use_lightning_w"          ""  ""
+rconfig   logical   use_lightning_qv        namelist,wrfvar4  1  .false.  - "use_lightning_qv"         ""  ""
+rconfig   logical   use_lightning_div       namelist,wrfvar4  1  .false.  - "use_lightning_div"        ""  ""
+rconfig   real      min_flashrate           namelist,wrfvar4  1   2.0     - "min_flashrate"            ""  ""
+rconfig   real      lightning_min_rh        namelist,wrfvar4  1   85.     - "lightning_min_rh"         ""  ""
 rconfig   logical   use_rainobs             namelist,wrfvar4  1  .false.  - "use_rainobs"              ""  ""   
 rconfig   logical   use_hirs2obs            namelist,wrfvar4  1  .false.  - "use_hirs2obs"             ""  ""
 rconfig   logical   use_hirs3obs            namelist,wrfvar4  1  .false.  - "use_hirs3obs"             ""  ""
@@ -230,6 +236,9 @@ rconfig   real      max_error_buv           namelist,wrfvar5  1   500.0   - "max
 rconfig   real      max_error_bt            namelist,wrfvar5  1   500.0   - "max_error_bt"             ""  ""
 rconfig   real      max_error_bq            namelist,wrfvar5  1   500.0   - "max_error_bq"             ""  ""
 rconfig   real      max_error_slp           namelist,wrfvar5  1   500.0   - "max_error_slp"            ""  ""
+rconfig   real      max_error_lda_w         namelist,wrfvar5  1   5.0     - "max_error_lda_w"          ""  ""
+rconfig   real      max_error_lda_div       namelist,wrfvar5  1   5.0     - "max_error_lda_div"        ""  ""
+rconfig   real      max_error_lda_qv        namelist,wrfvar5  1   5.0     - "max_error_lda_qv"         ""  ""
 rconfig   logical   check_buddy             namelist,wrfvar5  1  .false.  - "check_buddy"              ""  ""
 rconfig   logical   put_rand_seed           namelist,wrfvar5  1  .false.  - "put_rand_seed"            ""  ""
 rconfig   logical   omb_set_rand            namelist,wrfvar5  1  .false.  - "omb_set_rand"             ""  ""

frame/module_driver_constants.F

@@ -85,9 +85,9 @@ MODULE module_driver_constants
    !  The maximum number of obs indexes (for conventional DA obs) 
 
 #if (WRF_CHEM == 1)
-   INTEGER , PARAMETER :: num_ob_indexes  =   30
+   INTEGER , PARAMETER :: num_ob_indexes  =   31
 #else
-   INTEGER , PARAMETER :: num_ob_indexes  =   29
+   INTEGER , PARAMETER :: num_ob_indexes  =   30
 #endif
 
 

var/build/da.make

@@ -29,6 +29,7 @@ WRFVAR_OBJS = \
    da_pilot.o \
    da_radar.o \
    da_rain.o \
+   da_lightning.o \
    da_gpspw.o \
    da_gpsref.o \
    da_gpseph.o \

var/da/da_control/da_control.f90

@@ -240,6 +240,7 @@ module da_control
    real, parameter    :: typical_rv_rms = 1.0    ! m/s
    real, parameter    :: typical_rf_rms = 1.0    ! dBZ
    real, parameter    :: typical_rain_rms = 1.0   ! mm  
+   real, parameter    :: typical_div_rms = 0.001 
 
    ! The following typical mean squared values depend on control variable. They   
    ! are calculated in da_setup_background_errors and used in the VvToVp adjoint 
@@ -487,7 +488,7 @@ module da_control
 
    ! rtm_init setup parameter
 
-   integer, parameter            :: maxsensor = 30
+   integer, parameter            :: maxsensor = 31
 
    integer, parameter :: npres_print = 12
 
@@ -525,8 +526,9 @@ module da_control
    integer, parameter :: tamdar_sfc = 27
    integer, parameter :: rain      = 28
    integer, parameter :: gpseph    = 29
+   integer, parameter :: lightning = 30
 #if (WRF_CHEM == 1)
-   integer, parameter :: chemic_surf = 30
+   integer, parameter :: chemic_surf = 31
 #endif
 
    character(len=14), parameter :: obs_names(num_ob_indexes) = (/ &
@@ -558,7 +560,8 @@ module da_control
       "tamdar        ", &
       "tamdar_sfc    ", &
       "rain          ", &
-      "gpseph        "  &
+      "gpseph        ", &
+      "lightning     "  & 
 #if (WRF_CHEM == 1)
      ,"chemic_surf   " &
 #endif

var/da/da_define_structures/da_allocate_observations.inc

@@ -36,6 +36,7 @@ subroutine da_allocate_observations (iv)
    if (iv%info(profiler)%nlocal  > 0) allocate(iv%profiler (1:iv%info(profiler)%nlocal))
    if (iv%info(buoy)%nlocal      > 0) allocate(iv%buoy     (1:iv%info(buoy)%nlocal))
    if (iv%info(radar)%nlocal     > 0) allocate(iv%radar    (1:iv%info(radar)%nlocal))
+   if (iv%info(lightning)%nlocal > 0) allocate(iv%lightning(1:iv%info(lightning)%nlocal))
    if (iv%info(bogus)%nlocal     > 0) allocate(iv%bogus    (1:iv%info(bogus)%nlocal))
    if (iv%info(airsr)%nlocal     > 0) allocate(iv%airsr    (1:iv%info(airsr)%nlocal))
 

var/da/da_define_structures/da_allocate_y.inc

@@ -211,6 +211,19 @@ subroutine da_allocate_y (iv, y)
       end do
    end if
 
+   if (y % nlocal(lightning) > 0) then
+      allocate (y % lightning(1:y % nlocal(lightning)))
+      do n = 1, y % nlocal(lightning)
+         nlevels = iv%info(lightning)%levels(n)
+         allocate (y % lightning(n) % w(1:nlevels))
+         allocate (y % lightning(n) % div(1:nlevels))
+         allocate (y % lightning(n) % qv(1:nlevels))
+         y % lightning(n) % w(1:nlevels) = 0.0
+         y % lightning(n) % div(1:nlevels) = 0.0
+         y % lightning(n) % qv(1:nlevels) = 0.0 
+      end do
+   end if
+
    if (y % nlocal(airep) > 0) then
       allocate (y % airep(1:y % nlocal(airep)))
       do n = 1, y % nlocal(airep)

var/da/da_define_structures/da_allocate_y_lightning.inc

@@ -0,0 +1,44 @@
+subroutine da_allocate_y_lightning (iv, y)
+
+   !---------------------------------------------------------------------------
+   ! Purpose: Allocate arrays used in y and residual obs structures.
+   ! Authors: Z Chen (zchen@fjnu.edu.cn), Jenny Sun (NCAR), X Qie (IAP)   
+   !---------------------------------------------------------------------------
+
+   implicit none
+
+   type (iv_type), intent(in)    :: iv      ! Ob type input.
+   type (y_type),  intent(inout) :: y       ! Residual type structure.
+
+   integer                       :: n       ! Loop counter.
+   integer                       :: nlevels ! Number of levels.
+
+   !---------------------------------------------------------------------------
+   !  [1.0] Copy number of observations:
+   !---------------------------------------------------------------------------
+
+   if (trace_use) call da_trace_entry("da_allocate_y_lightning")
+
+   y % nlocal(lightning) = iv%info(lightning)%nlocal
+   y % ntotal(lightning) = iv%info(lightning)%ntotal
+
+  !---------------------------------------------------------------------------
+  ! [2.0] Allocate:
+  !---------------------------------------------------------------------------
+
+   if (y % nlocal(lightning) > 0) then
+      allocate (y % lightning(1:y % nlocal(lightning)))
+      do n = 1, y % nlocal(lightning)
+         nlevels = iv%info(lightning)%levels(n)
+         allocate (y % lightning(n) % w(1:nlevels))
+         allocate (y % lightning(n) % div(1:nlevels))
+         allocate (y % lightning(n) % qv(1:nlevels))
+         y % lightning(n) % w(1:nlevels) = 0.0
+         y % lightning(n) % div(1:nlevels) = 0.0
+         y % lightning(n) % qv(1:nlevels) = 0.0
+	   end do
+   end if
+
+   if (trace_use) call da_trace_exit("da_allocate_y_lightning")
+
+end subroutine da_allocate_y_lightning

var/da/da_define_structures/da_deallocate_observations.inc

@@ -157,6 +157,15 @@ subroutine da_deallocate_observations (iv)
       deallocate (iv%radar)
    end if
 
+   if (iv%info(lightning)%nlocal > 0) then
+      do n = 1, iv%info(lightning)%nlocal
+         deallocate (iv%lightning(n) % w)
+         deallocate (iv%lightning(n) % div)
+         deallocate (iv%lightning(n) % qv)
+      end do
+      deallocate (iv%lightning)
+   end if
+
    if (iv%info(rain)%nlocal > 0) deallocate (iv%rain)
    if (iv%info(gpspw)%nlocal > 0) deallocate (iv%gpspw)
 

var/da/da_define_structures/da_deallocate_y.inc

@@ -81,18 +81,26 @@ subroutine da_deallocate_y(y)
       deallocate (y % bogus)
    end if
 
-    if (y % nlocal(radar) > 0) then
-       do n = 1, y % nlocal(radar)
-          deallocate (y % radar(n)%rv)
-          deallocate (y % radar(n)%rf)
-          if (associated(y%radar(n)%rqv)) deallocate(y%radar(n)%rqv)
-          if (associated(y%radar(n)%rgr)) deallocate(y%radar(n)%rgr)
-          if (associated(y%radar(n)%rsn)) deallocate(y%radar(n)%rsn)
-          if (associated(y%radar(n)%rrn)) deallocate(y%radar(n)%rrn)
-       end do
-       deallocate (y % radar)
-    end if
+   if (y % nlocal(radar) > 0) then
+      do n = 1, y % nlocal(radar)
+         deallocate (y % radar(n)%rv)
+         deallocate (y % radar(n)%rf)
+         if (associated(y%radar(n)%rqv)) deallocate(y%radar(n)%rqv)
+         if (associated(y%radar(n)%rgr)) deallocate(y%radar(n)%rgr)
+         if (associated(y%radar(n)%rsn)) deallocate(y%radar(n)%rsn)
+         if (associated(y%radar(n)%rrn)) deallocate(y%radar(n)%rrn)
+      end do
+      deallocate (y % radar)
+   end if
 
+   if (y % nlocal(lightning) > 0) then
+      do n = 1, y % nlocal(lightning)
+         deallocate (y % lightning(n)%w)
+         deallocate (y % lightning(n)%div)
+         deallocate (y % lightning(n)%qv)
+      end do
+      deallocate (y % lightning)
+   end if
 
    if (y % nlocal(airep) > 0) then
       do n = 1, y % nlocal(airep)

var/da/da_define_structures/da_define_structures.f90

@@ -21,7 +21,7 @@ module da_define_structures
       put_rand_seed, seed_array1, seed_array2, missing_r, &
       sound, synop, pilot, satem, geoamv, polaramv, airep, gpspw, gpsref, gpseph, &
       metar, ships, ssmi_rv, ssmi_tb, ssmt1, ssmt2, qscat, profiler, buoy, bogus, &
-      mtgirs, tamdar, tamdar_sfc, pseudo, radar, radiance, airsr, sonde_sfc, rain, &
+      mtgirs, tamdar, tamdar_sfc, pseudo, radar, lightning, radiance, airsr, sonde_sfc, rain, &
 #if (WRF_CHEM == 1)
       chemic_surf, chem_cv_options, &
 #endif
@@ -318,6 +318,42 @@ module da_define_structures
       type (rain_each_type)                   :: each(1)
    end type rain_single_level_type
 
+   type lightning_stn_type
+      character (len = 5)    :: platform      ! Data type
+      character (len = 12)   :: name          ! Station name
+      character (len = 19)   :: date_char     ! CCYY-MM-DD_HH:MM:SS date
+      integer                :: numobs        ! number of Obs
+      integer                :: levels        ! number of levels
+      real                   :: lat           ! Latitude in degree
+      real                   :: lon           ! Longitude in degree
+      real                   :: elv           ! Elevation in 
+   end type lightning_stn_type
+
+   type lightning_type
+      type (stn_loc_type)             :: stn_loc  
+      real                  , pointer :: height   (:) ! Height in m
+      integer               , pointer :: height_qc(:) ! Height QC	  
+      type (field_type)     , pointer :: w(:)         ! Retrieved vertical velocity from flash rate
+      type (field_type)     , pointer :: div(:)       ! Retrieved convergence fileds from vertical velocity
+      type (field_type)     , pointer :: qv(:)        ! Retrieved vapor mixing ratio from flash rate         
+   end type lightning_type
+
+   type lightning_each_level_type
+      real                   :: height         ! Height in m
+      integer                :: height_qc      ! Height QC
+      real                   :: zk             ! MM5 k-coordinates
+      type (field_type)      :: w
+      type (field_type)      :: div
+      type (field_type)      :: qv
+   end type lightning_each_level_type
+
+   type lightning_multi_level_type
+      type (lightning_stn_type)               :: stn
+      type (info_type)                        :: info
+      type (model_loc_type)                   :: loc
+      type (lightning_each_level_type)        :: each(max_ob_levels)
+   end type lightning_multi_level_type
+
 #if (WRF_CHEM == 1)
 
    type chemic_surf_type
@@ -702,6 +738,7 @@ module da_define_structures
       real    :: bogus_ef_u, bogus_ef_v, bogus_ef_t, bogus_ef_p, bogus_ef_q, bogus_ef_slp
       real    :: airsr_ef_t,  airsr_ef_q
       real    :: rain_ef_r
+      real    :: lightning_ef_w, lightning_ef_div, lightning_ef_qv
 #if (WRF_CHEM == 1)
       real    :: chemic_surf_ef
 #endif
@@ -737,6 +774,7 @@ module da_define_structures
       type (tamdar_type)   , pointer :: tamdar(:)
       type (synop_type)    , pointer :: tamdar_sfc(:)
       type (rain_type)     , pointer :: rain(:)
+      type (lightning_type), pointer :: lightning(:)
 #if (WRF_CHEM == 1)
       type (chemic_surf_type), pointer :: chemic_surf(:)
 #endif
@@ -782,6 +820,8 @@ module da_define_structures
       type (bad_info_type)       :: slp
       type (bad_info_type)       :: rad
       type (bad_info_type)       :: rain
+      type (bad_info_type)       :: w
+      type (bad_info_type)       :: div
 #if (WRF_CHEM == 1)
       type (bad_info_type)       :: chemic_surf
 #endif
@@ -928,6 +968,12 @@ module da_define_structures
       real, pointer :: rqv(:) => null()
    end type residual_radar_type
 
+   type residual_lightning_type
+      real, pointer :: w(:)
+      real, pointer :: div(:)
+      real, pointer :: qv(:)
+   end type residual_lightning_type
+
    type residual_instid_type
       integer                          :: num_rad
       integer                          :: nchan
@@ -985,6 +1031,7 @@ module da_define_structures
       type (residual_radar_type),    pointer :: radar(:)
       type (residual_instid_type),   pointer :: instid(:)
       type (residual_rain_type),     pointer :: rain(:)
+      type (residual_lightning_type),pointer :: lightning(:)
 #if (WRF_CHEM == 1)
       type (residual_chem_surf_type),pointer :: chemic_surf(:)
 #endif
@@ -1039,6 +1086,7 @@ module da_define_structures
       real                :: bogus_u, bogus_v, bogus_t, bogus_q, bogus_slp
       real                :: airsr_t, airsr_q
       real                :: rain_r
+      real                :: lightning_w, lightning_div, lightning_qv
 #if (WRF_CHEM == 1)
       real                :: chemic_surf
 #endif
@@ -1216,6 +1264,7 @@ module da_define_structures
 #endif
 #include "da_allocate_y.inc"
 #include "da_allocate_y_radar.inc"
+#include "da_allocate_y_lightning.inc"
 #include "da_allocate_y_rain.inc"
 #if (WRF_CHEM == 1)
 #include "da_allocate_y_chem_sfc.inc"

var/da/da_define_structures/da_zero_y.inc

@@ -284,6 +284,17 @@ subroutine da_zero_y( iv, y, value )
       end do
    end if
 
+   ! Initialize lightning:
+   if ( y % nlocal(lightning) > 0 ) then
+      do n = 1, y % nlocal(lightning)
+         nlevels = iv % info(lightning) % levels(n)
+
+         y % lightning(n) % w(1:nlevels) = value
+         y % lightning(n) % div(1:nlevels) = value
+         y % lightning(n) % qv(1:nlevels) = value
+      end do
+   end if
+
    ! Initialize rain:
    if ( y % nlocal(rain) > 0 ) then
           y % rain(1:y % nlocal(rain)) % rain = value