var/da/da_obs_io/da_read_obs_lightning.inc

   1 subroutine da_read_obs_lightning (iv, filename, grid)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Read the lightning observation file
   5    ! Authors: Z Chen (zchen@fjnu.edu.cn), Jenny Sun (NCAR), X Qie (IAP)
   6    !----------------------------------------------------------------------------------------!
   7
   8    implicit none
   9
  10    type (iv_type),    intent(inout) :: iv
  11    character(len=*),  intent(in)    :: filename
  12    type(domain),     intent(in)     :: grid     ! first guess state.
  13
  14    integer                       :: i, n, iost
  15    integer                       :: iunit
  16
  17    integer                          :: i1, j1, k  ! Index dimension.
  18    real                             :: dx, dxm  ! Interpolation weights
  19    real                             :: dy, dym  ! Interpolation weights
  20    real                             :: zlcl
  21
  22    type (lightning_multi_level_type) :: platform
  23
  24    character (len = 120)         :: char_total_lightning, char_total_levels
  25    character (len = 160)         :: info_string
  26    integer                       :: total_lightning, nlevels, lightning_qc, rh_indicator
  27    real                          :: flashrate, wmax, lightning_error
  28    real, allocatable,dimension(:):: height, coff
  29    logical                       :: outside, outside_all
  30    integer                       :: nlocal
  31
  32    if (trace_use) call da_trace_entry("da_read_obs_lightning")
  33
  34    nlocal = 0
  35
  36    ! 1. open file
  37    call da_get_unit(iunit)
  38    open(unit   = iunit,     &
  39         FILE   = trim(filename), &
  40         FORM   = 'FORMATTED',  &
  41         ACCESS = 'SEQUENTIAL', &
  42         iostat =  iost,     &
  43         STATUS = 'OLD')
  44
  45    if (iost /= 0) then
  46       write(unit=message(1),fmt='(A,I5,A)') &
  47          "Error",iost," opening lightning obs file "//trim(filename)
  48       call da_warning(__FILE__,__LINE__,message(1:1))
  49       call da_free_unit(iunit)
  50       if (trace_use) call da_trace_exit("da_read_obs_lightning")
  51       return
  52    end if
  53
  54    ! 2. read basic info
  55
  56    ! 2.1 read the number of total lightning observation and vertical layers
  57    read (unit=iunit, fmt = '(A)', iostat = iost) char_total_lightning
  58    read (unit=iunit, fmt = '(A)', iostat = iost) char_total_levels
  59    read (unit=char_total_levels(9:15),fmt='(I7)', iostat = iost) nlevels
  60
  61    ! skip one line
  62    read (unit=iunit, fmt = '(A)', iostat = iost)
  63
  64    ! 2.2 read height and coefficient
  65    allocate(height(nlevels))
  66    allocate(coff(nlevels))
  67    do i = 1, nlevels
  68       read (unit = iunit, iostat = iost, fmt = '(2F12.3)') height(i), coff(i)
  69    end do
  70
  71    ! 2.3 read header info
  72    head_info: do
  73       read (unit=iunit, fmt = '(A)', iostat = iost) info_string
  74       if (iost /= 0) then
  75          write(unit=message(1),fmt='(A,I3,A,I3)') &
  76             "Error",iost,"reading lightning obs header on unit",iunit
  77          call da_warning(__FILE__,__LINE__,message(1:1))
  78       if (trace_use) call da_trace_exit("da_scan_obs_lightning")
  79       return
  80       end if
  81       if (info_string(1:6) == 'data  ') exit
  82    end do head_info
  83
  84    ! 2.4 read total lightning data info
  85    read (unit=char_total_lightning (8:14),fmt='(I7)', iostat = iost) total_lightning
  86
  87    ! 2.5 skip one line
  88    read (unit=iunit, fmt = '(A)', iostat = iost)
  89
  90    ! 3. read lightning data
  91    reports:   do n = 1, total_lightning
  92       ! 3.1 read station general info
  93       read (unit = iunit, iostat = iost, &
  94                    fmt = '(A12,1X,A19,1X,I6,2(F12.3,2X),F8.1,1X,A5)') &
  95                    platform%info%platform,  &
  96                    platform%info%date_char, &
  97                    platform%info%levels,    &
  98                    platform%info%lat,       &
  99                    platform%info%lon,       &
 100                    platform%info%elv,       &
 101                    platform%info%id
 102
 103      call da_llxy (platform%info, platform%loc, outside, outside_all)
 104
 105 ! Height information is from xb, get horizontal interpolation weights:
 106       i1   = platform%loc%i
 107       j1   = platform%loc%j
 108       dx   = platform%loc%dx
 109       dy   = platform%loc%dy
 110       dxm  = platform%loc%dxm
 111       dym  = platform%loc%dym
 112
 113       ! 3.2 read lightning flash rate and its qc and error info
 114       read (unit = iunit, fmt = '(F12.3,I4,F12.3,I4)') flashrate, lightning_qc, lightning_error, rh_indicator
 115
 116       !turn lighting flash rate into the maximum wmax
 117       if(flashrate .ge. min_flashrate .and. flashrate .le. min_flashrate+10.0)then
 118          wmax = 5!14.6
 119       end if
 120       if(flashrate .gt.min_flashrate+10.0 .and. flashrate .le. min_flashrate+20.0)then
 121          wmax = 8!17.07
 122       end if
 123       if(flashrate .gt.min_flashrate+20.0 .and. flashrate.le.min_flashrate+30.0)then
 124          wmax = 12!18.67
 125       end if
 126       if(flashrate .gt.min_flashrate+30.0)then
 127          wmax = 15!24.4 !m/s
 128       end if
 129
 130       zlcl = 125.0*(grid%xb%t(i1,j1,1)-grid%xb%td(i1,j1,1)) + grid%xb%terr(i1,j1)
 131       zlcl = amax1(grid%xb%terr(i1,j1)+1000.,zlcl)
 132       zlcl = amin1(3000.     ,zlcl)
 133
 134       do i = 1, nlevels !vertical layers
 135          platform%each(i) = lightning_each_level_type(missing_r, missing, -1.0,       &
 136                field_type(missing_r, missing, missing_r, missing, missing_r), & ! w
 137                field_type(missing_r, missing, missing_r, missing, missing_r), & ! div
 138                field_type(missing_r, missing, missing_r, missing, missing_r))   ! qv
 139
 140         platform%each(i)%height = grid%xb%h(i1,j1,k)    !height(i)
 141
 142         if(flashrate .ge. min_flashrate .and. i .gt. 1) then
 143            ! vertical velocity
 144            platform%each(i)%w%inv = wmax*coff(i)
 145            platform%each(i)%w%qc = 0
 146            platform%each(i)%w%error = amax1(1.0, 0.20*abs(platform%each(i)%w%inv))
 147            ! divergence
 148            platform%each(i)%div%inv = -wmax*(coff(i)-coff(i-1))/(height(i)-height(i-1))
 149            platform%each(i)%div%qc = 0
 150            platform%each(i)%div%error = amax1(0.0001, 0.20*abs(platform%each(i)%div%inv))
 151         else
 152            platform%each(i)%w%qc = missing_data
 153            platform%each(i)%w%error = missing_r
 154            platform%each(i)%div%qc = missing_data
 155            platform%each(i)%div%error = missing_r
 156         end if
 157
 158
 159         if(flashrate .ge. 10.0 .and. rh_indicator .gt. -1 .and. height(i) .ge. zlcl .and. height(i) .le.15000)then
 160            platform%each(i)%qv%inv = 0.01*(2*rh_indicator+lightning_min_rh)*grid%xb%qs(i1,j1,i)
 161            platform%each(i)%qv%qc = 0
 162            platform%each(i)%qv%error = amax1(0.001,0.20*platform%each(i)%qv%inv)
 163         else
 164            platform%each(i)%qv%qc = missing_data
 165            platform%each(i)%qv%error = missing_r
 166         end if
 167       end do  !vertical layers
 168
 169       if(outside)then
 170         cycle
 171       end if
 172       nlocal = nlocal+1
 173
 174       iv%info(lightning)%levels(nlocal)    = nlevels
 175       iv%info(lightning)%name(nlocal)      = platform%info%name
 176       iv%info(lightning)%platform(nlocal)  = platform%info%platform
 177       iv%info(lightning)%id(nlocal)        = platform%info%id
 178       iv%info(lightning)%date_char(nlocal) = platform%info%date_char
 179       iv%info(lightning)%lat(:,nlocal)     = platform%info%lat
 180       iv%info(lightning)%lon(:,nlocal)     = platform%info%lon
 181       iv%info(lightning)%elv(nlocal)       = platform%info%elv
 182       iv%info(lightning)%pstar(nlocal)     = platform%info%pstar
 183
 184       iv%info(lightning)%slp(nlocal)       = platform%loc%slp
 185       iv%info(lightning)%pw(nlocal)        = platform%loc%pw
 186       iv%info(lightning)%x(:,nlocal)       = platform%loc%x
 187       iv%info(lightning)%y(:,nlocal)       = platform%loc%y
 188       iv%info(lightning)%i(:,nlocal)       = platform%loc%i
 189       iv%info(lightning)%j(:,nlocal)       = platform%loc%j
 190       iv%info(lightning)%dx(:,nlocal)      = platform%loc%dx
 191       iv%info(lightning)%dxm(:,nlocal)     = platform%loc%dxm
 192       iv%info(lightning)%dy(:,nlocal)      = platform%loc%dy
 193       iv%info(lightning)%dym(:,nlocal)     = platform%loc%dym
 194       iv%info(lightning)%proc_domain(:,nlocal) = platform%loc%proc_domain
 195       iv%info(lightning)%obs_global_index(nlocal) = nlocal
 196
 197       allocate(iv%lightning(nlocal)%height(1:nlevels))
 198       allocate(iv%lightning(nlocal)%height_qc(1:nlevels))
 199       allocate(iv%lightning(nlocal)%w(1:nlevels))
 200       allocate(iv%lightning(nlocal)%qv(1:nlevels))
 201       allocate(iv%lightning(nlocal)%div(1:nlevels))
 202
 203       do i = 1, nlevels
 204         iv%lightning(nlocal)%height(i)   = platform%each(i)%height
 205         iv%lightning(nlocal)%height_qc(i)= platform%each(i)%height_qc
 206         iv%lightning(nlocal)%w(i)        = platform%each(i)%w
 207         iv%lightning(nlocal)%qv(i)       = platform%each(i)%qv
 208         iv%lightning(nlocal)%div(i)      = platform%each(i)%div
 209       end do
 210
 211    end do reports
 212    deallocate(height)
 213    deallocate(coff)
 214    close(iunit)
 215    call da_free_unit(iunit)
 216
 217    if (trace_use) call da_trace_exit("da_read_obs_lightning")
 218
 219
 220 end subroutine da_read_obs_lightning