Skip to content

Add Dynamic lightning scheme #2276

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
weiwangncar wants to merge 6 commits into
base: develop
Choose a base branch
from
+894 −6
Open

Add Dynamic lightning scheme #2276

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c75c7e3
Add the Dynamic Lightning Scheme from Barry Lynn, Weather It Is, LTD …
weiwangncar Jan 16, 2026
c7aa6c7
Add modified files to the dynamic lightning scheme
weiwangncar Jan 16, 2026
769e915
make new variables optional so that wrfplus can compile
weiwangncar Jan 17, 2026
4d3774a
Add new namelists to README file
weiwangncar Jan 25, 2026
f5cfd96
Update solve_em.F
weiwangncar Jan 29, 2026
f6b33ec
remove unused var (l_obs), commented out print, and a minor change
weiwangncar Jan 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 32 additions & 0 deletions Registry/registry.dyn_light
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
rconfig integer dyn_lightning_option namelist,physics max_domains 0 rh "dynlight" "Dynamic Lightning Scheme option, 1: on"


#dynamic lightning
rconfig integer num_light_periods namelist,physics 1 1 irh "num_light_periods" "" ""
rconfig real coul_pos namelist,physics max_domains 0.25E-4 rh "Constant for Positive Charging of Clouds" "Coulombs" ""
rconfig real coul_neg namelist,physics max_domains 0.25E-4 rh "Constant for Negative Charging of Clouds" "Coulombs" ""
rconfig real coul_neu namelist,physics max_domains 0.25E-4 rh "Constant for Charging of Intracloud Lightning" "Coulombs" ""
rconfig real j_pos namelist,physics max_domains 5.E9 rh "Threshold for Producing Positive Event" "J" ""
rconfig real j_neg namelist,physics max_domains 1.E9 rh "Threshold for Producing Negative Event" "J" ""
rconfig real j_neu namelist,physics max_domains 1.E9 rh "Threshold for Producing IC Event" "" "J"
# end dynamic lightning

# state real - ikjftb scalar 1 - - -
state real light_ne ikjftb scalar 1 - \
i0rhusdf=(bdy_interp:dt) "LIGHTNING_NE" "LIGHTNING NEG TOTAL ENERGY " "J"
state real light_pe ikjftb scalar 1 - \
i0rhusdf=(bdy_interp:dt) "LIGHTNING_PE" "LIGHTNING POS TOTAL ENERGY" "J"
state real light_neu ikjftb scalar 1 - \
i0rhusdf=(bdy_interp:dt) "LIGHTNING_NEU" "LIGHTNING NEU TOTAL ENERGY" "J"
#
state real LPOS ij misc 1 - irh05du "LPOS" "Positive Cloud to Ground Lightning Density" "# time-l"
state real LNEG ij misc 1 - irh05du "LNEG" "Negative Cloud to Ground Lightning Density" "# time-1"
state real LNEU ij misc 1 - irh05du "LNEU" "Intra-Cloud Lightning Density" "# time-1"
state real LPI2D ij misc 1 - rh05 "LPI2D" "Lightning Potential Index (2D)" "J/kg"
state real LPI3D ikj misc 1 - rh05 "LPI3D" "Lightning Potential Index (3D)" "J/kg"
state real l_obs i{lp}j misc 1 - irh05 "L_OBS" "Lightning OBS" " "

#package dynlight_output dyn_lightning_option==1 - scalar:light_pe,light_ne,light_neu;state:lneu,lneg,lpos,l_obs,lpi2d,lpi3d
package dynlight_output dyn_lightning_option==1 - scalar:light_pe,light_ne,light_neu;state:lneu,lneg,lpos,lpi2d,lpi3d


1 change: 1 addition & 0 deletions Registry/registry.em_shared_collection
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,5 @@ include registry.trad_fields
include registry.solar_fields
include registry.diags
include registry.iau
include registry.dyn_light
include registry.CMAQ
10 changes: 8 additions & 2 deletions dyn_em/solve_em.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3983,8 +3983,14 @@ END SUBROUTINE CMAQ_DRIVER
& ,pert_thom_qc=config_flags%pert_thom_qc &
& ,pert_thom_qi=config_flags%pert_thom_qi &
& ,pert_thom_qs=config_flags%pert_thom_qs &
& ,pert_thom_ni=config_flags%pert_thom_ni &
& ,cloudnc=grid%cloudnc &
& ,pert_thom_ni=config_flags%pert_thom_ni &
& ,cloudnc=grid%cloudnc &
! Dynamic Lightning Scheme
& ,light_pe_curr=scalar(ims,kms,jms,P_light_pe), F_light_pe=F_light_pe &
& ,light_ne_curr=scalar(ims,kms,jms,P_light_ne), F_light_ne=F_light_ne &
& ,light_neu_curr=scalar(ims,kms,jms,P_light_neu), F_light_neu=F_light_neu &
& ,lpos=grid%lpos,lneg=grid%lneg,lneu=grid%lneu &
& ,lpi2d=grid%lpi2d,lpi3d=grid%lpi3d &
)


Expand Down
3 changes: 3 additions & 0 deletions phys/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,9 @@ MODULES = \
module_ltng_cpmpr92z.o \
module_ltng_crmpr92.o \
module_ltng_iccg.o \
module_calc_lpi_new.o \
module_ltng_pe.o \
module_ltng_strokes.o \
module_fdda_psufddagd.o \
module_fdda_spnudging.o \
module_fddagd_driver.o \
Expand Down
215 changes: 215 additions & 0 deletions phys/module_calc_lpi_new.F
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
MODULE module_calc_lpi_new
CONTAINS
SUBROUTINE calc_lpi_new(qc, qr, qi, qs, qg,qh &
,w,dz8w,pi_phy,th_phy,p_phy,rho_phy &
,lpi2d,lpi3d &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
)

IMPLICIT NONE
!-------------------------------------------------------------------
!
!
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN) :: &
qc, &
qi, &
qr, &
qs, &
qg, &
qh
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT),OPTIONAL :: &
lpi3d
! ,&
! light_ne, &

REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: w
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN) :: th_phy
!

!
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: &
rho_phy, &
dz8w, &
pi_phy, &
p_phy

! REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme):: &
! & prec_ice,wqg_m15
REAL, DIMENSION(ims:ime,jms:jme):: &
& lpi2d
! REAL, DIMENSION(ims:ime,jms:jme):: &
! & LPI,LPOS,LNEG,LNEU,plpi,nlpi,neulpi,prec_ice,wqg_m15




REAL, DIMENSION(kms:kme):: tempk,rh,qtot,qitot
REAL, DIMENSION(kms:kme):: p1d,rho1d,qti1d
REAL, DIMENSION(kms:kme):: temp,qc1d,ql1d,qi1d,qs1d,qg1d,lpi1d,del_z
REAL, DIMENSION(0:kme):: w1d,height
REAL, DIMENSION(kms:kme):: e1d,height_t,w1d_t
REAL z_full,qrs,teten,RELHUM,LOC,Td_850,Td_700,PC_DWPT
INTEGER level
REAL :: t_base,t_top
real top, bot
real num,den,ave_z
real num_s,den_s
real num_i,den_i
real q_isg,del_z_tot

INTEGER I_COLLAPSE
LOGICAL LOOK_T
INTEGER I_START,I_END,J_START,J_END


INTEGER :: i,j,k
!-------------------------------------------------------------------
! ! print*,'ims,ime,jms,jme,kms,kme = ',ims,ime,jms,jme,kms,kme
! print*,'its,ite,jts,jte,kts,kte = ',its,ite,jts,jte,kts,kte
! parameter(t_pos=0.12,t_neg=0.06,t_neu=0.03)
! note: t_pos is modified based on the height of the top of the cloud
! I believe that the value of 60000 m/s might be a typical speed
! See: https\://en.wikipedia.org/wiki/Lightning
! print*,'dt = ',dt
! print*,'dx = ',dx
! print*,'coul_pos = ',coul_pos
! print*,'coul_neg = ',coul_neg
! print*,'coul_neu = ',coul_neu
! return (1)
! if (MAXVAL(light_pe).ne.0)print*,'maxval light_pe'
! if (MAXVAL(light_pe).ne.0)write(MAXVAL(light_pe))
! if (MAXVAL(light_ne).ne.0)print*,'maxval light_ne'
! if (MAXVAL(light_ne).ne.0)write(MAXVAL(light_ne))
!-------------------------------------------------------------------
t_base = 0
t_top = -20.

DO j = jts,jte
DO i = its,ite
z_full=0.
height(0)=z_full
w1d(0)=w(i,1,j)
do k = kts,kte-1
if (k.lt.kte-1)then
w1d(k)=w(i,k+1,j)
else
w1d(k)=0.
end if
temp(k) = th_phy(i,k,j)*pi_phy(i,k,j)-273.16
tempk(k) = th_phy(i,k,j)*pi_phy(i,k,j)
p1d(k)=p_phy(i,k,j)
rho1d(k)=rho_phy(i,k,j)
z_full=z_full+dz8w(i,k,j)
height(k)=z_full
qc1d(k)=qc(i,k,j)
ql1d(k)=qc(i,k,j)+qr(i,k,j)
qi1d(k)=qi(i,k,j)
! qti1d(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j)+qh(i,k,j)
qti1d(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j)+qh(i,k,j)
qs1d(k)=qs(i,k,j)
! qg1d(k)=qg(i,k,j)+qh(i,k,j)
! Hail doesn't usually charge
qg1d(k)=qg(i,k,j) + qh(i,k,j)
! For conservative advection multiply by rho1d and divide by it below
enddo
do k = kts,kte-1
w1d_t(k)=0.5*(w1d(k-1)+w1d(k))
end do
do k=kts,kte
lpi3d(i,k,j) = 0
top=height(k+1)
bot=height(k)
del_z(k)=top-bot
! if (max_w(i,j).gt.0.5)print*,"max_w > 0.5"
end do

ave_z = 0
del_z_tot = 0
lpi2d(i,j) = 0

do k=kts,kte
if (temp(k).le.t_base.and.temp(k).gt.t_top)then ! set t1d range

den_i = qi1d(k)+qg1d(k)
den_s = qs1d(k)+qg1d(k)
if (qs1d(k).eq.0.or.qg1d(k).eq.0.)then !checks for zeroes
den_s=10000.
num_s = 0.
else
num_s = sqrt(qs1d(k)*qg1d(k))
end if
if (qi1d(k).eq.0.or.qg1d(k).eq.0.)then ! checks for zeroes
den_i=10000.
num_i = 0.
else
num_i = sqrt(qi1d(k)*qg1d(k))
end if
q_isg = qg1d(k)*(num_i/den_i+num_s/den_s) ! ice "fract"-content

if (ql1d(k).eq.0.or.q_isg.eq.0)then
num=0
den=10000.
else
num = sqrt(ql1d(k)*q_isg)
den = ql1d(k)+q_isg
end if
del_z_tot=del_z_tot+del_z(k)
if (num.gt.0)then
ave_z=ave_z+del_z(k)*(2.*num/den)*w1d_t(k)**2 ! power index J/unit-mass
end if
end if
if (lpi2d(i,j).lt.0)lpi2d(i,j)=0
end do
!
if (del_z_tot.eq.0)del_z_tot=100000
lpi2d(i,j)=ave_z/del_z_tot
ave_z = 0.
do k=kts,kte

den_i = qi1d(k)+qg1d(k)
den_s = qs1d(k)+qg1d(k)
if (qs1d(k).eq.0.or.qg1d(k).eq.0.)then !checks for zeroes
den_s=10000.
num_s = 0.
else
num_s = sqrt(qs1d(k)*qg1d(k))
end if
if (qi1d(k).eq.0.or.qg1d(k).eq.0.)then ! checks for zeroes
den_i=10000.
num_i = 0.
else
num_i = sqrt(qi1d(k)*qg1d(k))
end if
q_isg = qg1d(k)*(num_i/den_i+num_s/den_s) ! ice "fract"-content

if (ql1d(k).eq.0.or.q_isg.eq.0)then
num=0
den=10000.
else
num = sqrt(ql1d(k)*q_isg)
den = ql1d(k)+q_isg
end if
del_z_tot=del_z_tot+del_z(k)
if (num.gt.0)then
ave_z=(2.*num/den)*w1d_t(k)**2 ! power index J/unit-mass
lpi3d(i,k,j) = ave_z
end if
end do
end do
end do
! check cell for active cell
! go to 10!
! Check within five grid elements that a majority of max_w > 0.5 m/s
return
end subroutine calc_lpi_new
END MODULE module_calc_lpi_new
Loading