Add lightning data assimilation

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