SUBROUTINE convert_lghtn2ref_nmmb(mype,nlon,nlat,nsig,ref_mos_3d,lightning,h_bk)
!
!$$$ subprogram documentation block
! . . . .
! subprogram: convert_lghtn2ref convert lightning stroke rate to radar reflectivity
!
! PRGMMR: Ming Hu ORG: GSD/AMB DATE: 2012-10-16
!
! ABSTRACT:
! This subroutine converts lightning flash density to radar reflectivity based
! on Jing's statistic analysis
! PROGRAM HISTORY LOG:
! 2015-10-06 S.Liu Add NCO document block
! 2015-10-06 s.liu -add new algorithm from Jing Her to retrieve REF from lghtn for NMMB
! 2015-10-26 s.liu -reduce estimated reflectivity, appears the current
! algorithm overestimated ref (5dBz)
!
!
!
! input argument list:
! mype - processor ID
! nlon - no. of lons on subdomain (buffer points on ends)
! nlat - no. of lats on subdomain (buffer points on ends)
! nsig - no. of levels
! ref_mos_3d - 3D reflectivity in analysis grid
! lightning - 2D lightning flash rate in analysis grid
! h_bk - 3D height
!
! output argument list:
! ref_mos_3d - 3D reflectivity in analysis grid
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,i_kind,r_single
implicit none
INTEGER(i_kind),intent(in) :: mype
INTEGER(i_kind),intent(in) :: nlon,nlat,nsig
real(r_single), intent(in) :: h_bk(nlon,nlat,nsig) ! height
real(r_kind), intent(inout) :: lightning(nlon,nlat)
real(r_kind), intent(inout):: ref_mos_3d(nlon,nlat,nsig) ! reflectivity in grid
!
! local
!
real(r_kind) :: dbz_lightning(nlon,nlat)
real(r_kind) :: table_lghtn2ref_winter(30) ! table content the map from lightning strakes
! to maximum reflectivity
DATA table_lghtn2ref_winter/ &
32.81,33.98,34.93,36.26,36.72,37.07,37.93,38.79,39.65,40.10, &
40.42,41.42,41.90,42.04,42.19,42.45,42.90,43.20,43.50,43.80, &
44.10,44.66,44.84,45.56,45.64,45.80,45.95,46.11,46.32,46.50/
real(r_kind) :: table_lghtn2ref_summer(30) ! table content the map from lightning strakes
! to maximum reflectivity
DATA table_lghtn2ref_summer/ &
30.13,31.61,32.78,33.86,34.68,35.34,36.13,36.15,37.02,37.04, &
37.74,38.00,38.56,38.85,39.10,39.37,39.78,39.98,40.64,41.33, &
41.50,41.65,41.85,42.08,42.77,43.03,43.26,43.53,43.74,43.73/
integer(i_kind) :: maxlvl
parameter (maxlvl=31)
real(r_kind) :: newlvlAll(maxlvl) ! vertical levels of reflectivity statistic profile
DATA newlvlAll/0.2, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, &
3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, &
9, 10, 11, 12, 13, 14, 15, 16/
real(r_kind) :: refprofile_winter(maxlvl,4) ! statistic reflectivity profile used to
! retrieve vertical ref based on lightning
! max reflectivity 30-35 dbz
DATA refprofile_winter(:,1) / &
0.966,0.958,0.977,0.989,0.998,1.000,0.997,0.992,0.981,0.962, &
0.933,0.898,0.826,0.752,0.687,0.626,0.578,0.547,0.522,0.526, &
0.519,0.501,0.482,0.464,0.437,0.430,0.454,0.539,0.662,0.742, &
0.793/
! max reflectivity 35-40 dbz
DATA refprofile_winter(:,2) / &
0.947,0.953,0.980,0.994,1.000,0.996,0.987,0.974,0.956,0.928, &
0.891,0.848,0.761,0.679,0.613,0.559,0.522,0.491,0.473,0.462, &
0.451,0.433,0.415,0.403,0.382,0.380,0.406,0.482,0.603,0.707, &
0.723/
! max reflectivity 40-45 dbz
DATA refprofile_winter(:,3) / &
0.937,0.955,0.986,1.000,0.997,0.995,0.988,0.978,0.957,0.920, &
0.871,0.824,0.735,0.654,0.584,0.518,0.465,0.442,0.435,0.412, &
0.398,0.385,0.376,0.360,0.340,0.350,0.377,0.446,0.551,0.625, &
0.656/
! max reflectivity 45-50 dbz
DATA refprofile_winter(:,4) / &
0.900,0.949,0.982,0.995,1.000,0.998,0.983,0.954,0.914,0.874, &
0.834,0.793,0.721,0.664,0.612,0.565,0.530,0.496,0.460,0.431, &
0.402,0.383,0.370,0.354,0.335,0.321,0.347,0.342,0.441,0.510, &
0.548/
real(r_kind) :: refprofile_summer(maxlvl,4) ! statistic reflectivity profile used to
! retrieve vertical ref based on lightning
! max reflectivity 30-35 dbz
DATA refprofile_summer(:,1) / &
0.870,0.885,0.914,0.931,0.943,0.954,0.967,0.975,0.982,0.989, &
0.995,1.000,0.998,0.973,0.918,0.850,0.791,0.735,0.690,0.657, &
0.625,0.596,0.569,0.544,0.510,0.479,0.461,0.460,0.477,0.522, &
0.570/
! max reflectivity 35-40 dbz
DATA refprofile_summer(:,2) / &
0.871,0.895,0.924,0.948,0.961,0.971,0.978,0.983,0.988,0.992, &
0.997,1.000,0.995,0.966,0.913,0.848,0.781,0.719,0.660,0.611, &
0.576,0.542,0.523,0.513,0.481,0.448,0.416,0.402,0.417,0.448, &
0.491/
! max reflectivity 40-45 dbz
DATA refprofile_summer(:,3) / &
0.875,0.895,0.914,0.936,0.942,0.951,0.964,0.979,0.990,0.998, &
1.000,0.992,0.961,0.905,0.834,0.772,0.722,0.666,0.618,0.579, &
0.545,0.518,0.509,0.483,0.419,0.398,0.392,0.403,0.423,0.480, &
0.440/
! max reflectivity 45-50 dbz
DATA refprofile_summer(:,4) / &
0.926,0.920,0.948,0.975,0.988,0.989,0.995,0.997,1.000,1.000, &
0.997,0.991,0.970,0.939,0.887,0.833,0.788,0.741,0.694,0.655, &
0.611,0.571,0.551,0.537,0.507,0.470,0.432,0.410,0.420,0.405, &
0.410/
INTEGER(i_kind) :: season ! 1= summer, 2=winter
INTEGER(i_kind) :: num_lightning
INTEGER(i_kind) :: i,j, k2, k, mref
REAL(r_kind) :: heightGSI,upref,downref,wght
INTEGER(i_kind) :: ilvl,numref
REAL(r_kind) :: lowest,highest,tempref, tempprofile(maxlvl)
!
! map lightning strokes to maximum reflectiivty
!
season=1
dbz_lightning = -9999.0_r_kind
DO j=2,nlat-1
DO i=2,nlon-1
if(lightning(i,j) > 1.0_r_kind ) then
num_lightning = max(1,min(30,int(lightning(i,j))))
if(season== 2 ) then
dbz_lightning(i,j) = 7.62*log10(lightning(i,j))+30.0
else if(season== 1 ) then
dbz_lightning(i,j) = 7.62*log10(lightning(i,j))+30.0
endif
endif
ENDDO
ENDDO
lightning = -999.0
DO j=2,nlat-1
DO i=2,nlon-1
lightning(i,j) = dbz_lightning(i,j)
ENDDO
ENDDO
!
! vertical reflectivity distribution
!
DO k=1,maxlvl
newlvlAll(k)=newlvlAll(k)*1000.0_r_kind
ENDDO
! ref_mos_3d=-9999.0
DO j=2,nlat-1
DO i=2,nlon-1
if( dbz_lightning(i,j) > 30 ) then
mref = min(4,(int((dbz_lightning(i,j) - 30.0_r_kind)/5.0_r_kind) + 1 ))
if(season== 2 ) then
DO k=1,maxlvl
tempprofile(k)=refprofile_winter(k,mref)*dbz_lightning(i,j)
enddo
lowest=newlvlAll(2)
highest=7000.0_r_kind
else if(season== 1 ) then
DO k=1,maxlvl
tempprofile(k)=refprofile_summer(k,mref)*dbz_lightning(i,j)
enddo
lowest=newlvlAll(3)
highest=12000.0_r_kind
endif
DO k2=1,nsig
heightGSI=h_bk(i,j,k2)
if(heightGSI >= lowest .and. heightGSI < highest) then ! lower 12km ?
do k=1,maxlvl-1
if( heightGSI >=newlvlAll(k) .and. heightGSI < newlvlAll(k+1) ) ilvl=k
enddo
upref=tempprofile(ilvl+1)
downref=tempprofile(ilvl)
wght=(heightGSI-newlvlAll(ilvl))/(newlvlAll(ilvl+1)-newlvlAll(ilvl))
tempref=(1-wght)*downref + wght*upref
ref_mos_3d(i,j,k2) = max(ref_mos_3d(i,j,k2),tempref)
endif
ENDDO
endif
ENDDO
ENDDO
END SUBROUTINE convert_lghtn2ref_nmmb