diff --git a/Registry/registry.dyn_light b/Registry/registry.dyn_light new file mode 100644 index 0000000000..e5dd202f64 --- /dev/null +++ b/Registry/registry.dyn_light @@ -0,0 +1,31 @@ +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 + + diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 0761df036f..88f1ba77ed 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -183,8 +183,7 @@ SUBROUTINE solve_em ( grid , config_flags & LOGICAL :: leapfrog INTEGER :: l,kte,kk LOGICAL :: f_flux ! flag for computing averaged fluxes in cu_gd - REAL :: curr_secs, curr_secs2, curr_mins2 - REAL(8) :: curr_secs_r8, curr_secs2_r8 + REAL :: curr_secs, curr_secs2 INTEGER :: num_sound_steps INTEGER :: idex, jdex REAL :: max_msft @@ -199,7 +198,6 @@ SUBROUTINE solve_em ( grid , config_flags & TYPE(WRFU_TimeInterval) :: tmpTimeInterval, tmpTimeInterval2 REAL :: real_time - REAL(8) :: real_time_r8 LOGICAL :: adapt_step_flag LOGICAL :: fill_w_flag @@ -333,9 +331,6 @@ END SUBROUTINE CMAQ_DRIVER tmpTimeInterval2 = domain_get_current_time ( grid ) - domain_get_start_time ( grid ) curr_secs = real_time(tmpTimeInterval) curr_secs2 = real_time(tmpTimeInterval2) - curr_secs_r8 = real_time_r8(tmpTimeInterval) - curr_secs2_r8 = real_time_r8(tmpTimeInterval2) - curr_mins2 = REAL( curr_secs2_r8 / 60.0d0 ) old_dt = grid%dt ! store old time step for flux averaging code at end of RK loop @@ -816,7 +811,7 @@ END SUBROUTINE CMAQ_DRIVER , ph_tendf, mu_tendf & , tke_tend & , config_flags%use_adaptive_time_step & - , curr_secs, curr_mins2 & + , curr_secs & , psim , psih , gz1oz0 & , chklowq & , cu_act_flag , hol , th_phy & @@ -3729,7 +3724,6 @@ END SUBROUTINE CMAQ_DRIVER & ,P8W=p8w ,P=p_phy ,PI_PHY=pi_phy & & ,RHO=grid%rho ,SPEC_ZONE=grid%spec_zone & & ,SR=grid%sr ,TH=th_phy & - & ,ssat=grid%ssat, ssati=grid%ssati & & ,refl_10cm=grid%refl_10cm & ! hm, 9/22/09 for refl & ,vmi3d=grid%vmi3d & ! for P3 & ,di3d=grid%di3d & ! for P3 @@ -3763,7 +3757,7 @@ END SUBROUTINE CMAQ_DRIVER & ,DGNUM4D=grid%dgnum4d,DGNUMWET4D=grid%dgnumwet4d & !====================== #endif - & ,XLAND=grid%xland,SNOWH=grid%SNOW,XICE=grid%XICE & + & ,XLAND=grid%xland,SNOWH=grid%SNOW & !PMA & ,SPECIFIED=specified_bdy, CHANNEL_SWITCH=channel_bdy & & ,F_RAIN_PHY=grid%f_rain_phy & & ,F_RIMEF_PHY=grid%f_rimef_phy & @@ -3983,9 +3977,16 @@ 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 & +! 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 & + & ,l_obs=grid%l_obs& + & ,lpi2d=grid%lpi2d,lpi3d=grid%lpi3d & + & ,num_light_periods=config_flags%num_light_periods & + & ,curr_secs=curr_secs ) BENCH_END(micro_driver_tim) diff --git a/phys/Makefile b/phys/Makefile index 5eda61c111..aedb7dcff0 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -43,9 +43,9 @@ MODULES = \ module_bl_myjpbl.o \ module_bl_qnsepbl.o \ module_bl_acm.o \ - module_bl_mynnedmf_common.o \ - module_bl_mynnedmf.o \ - module_bl_mynnedmf_driver.o \ + module_bl_mynn_common.o \ + module_bl_mynn.o \ + module_bl_mynn_wrapper.o \ module_bl_fogdes.o \ module_bl_gwdo.o \ module_bl_gwdo_gsl.o \ @@ -97,7 +97,6 @@ MODULES = \ module_mp_etanew.o \ module_mp_fer_hires.o \ module_mp_thompson.o \ - module_mp_rcon.o \ module_fire_emis.o \ module_mp_SBM_polar_radar.o \ module_mp_full_sbm.o \ @@ -114,7 +113,6 @@ MODULES = \ module_mp_wdm5.o \ module_mp_wdm6.o \ module_mp_wdm7.o \ - module_mp_udm.o \ module_mp_ntu.o \ module_mp_cammgmp_driver.o \ module_ra_sw.o \ @@ -197,6 +195,9 @@ MODULES = \ module_radiation_driver.o \ module_surface_driver.o \ module_lightning_driver.o \ + module_calc_lpi_new.o \ + module_ltng_pe.o \ + module_ltng_strokes.o \ module_ltng_cpmpr92z.o \ module_ltng_crmpr92.o \ module_ltng_iccg.o \ @@ -251,7 +252,6 @@ LIBTARGET = physics TARGETDIR = ./ $(LIBTARGET) : - (cd .. && ./tools/manage_externals/checkout_externals --externals ./arch/Externals.cfg) $(MAKE) $(J) non_nmm ; \ $(AR) $(ARFLAGS) ../main/$(LIBWRFLIB) $(MODULES) $(OBJS) \ $(FIRE_MODULES) $(DIAGNOSTIC_MODULES_EM) $(PHYSMMM_MODULES) @@ -273,17 +273,7 @@ submodules : else \ echo No action required for NoahMP submodule ; \ fi - @if [ \( ! -f module_bl_mynnedmf.F \) -o \( ! -f module_bl_mynedmf_common.F \) -o \ - \( ! -f module_bl_mynnedmf_driver.F \) ] ; then \ - echo Pulling in MYNN-EDMF submodule ; \ - ( cd .. ; git submodule update --init --recursive ) ; \ - ln -sf MYNN-EDMF/module_bl_mynnedmf.F90 module_bl_mynnedmf.F ; \ - ln -sf MYNN-EDMF/WRF/module_bl_mynnedmf_common.F90 module_bl_mynnedmf_common.F ; \ - ln -sf MYNN-EDMF/WRF/module_bl_mynnedmf_driver.F90 module_bl_mynnedmf_driver.F ; \ - else \ - echo No action required for MYNN-EDMF submodule ; \ - fi - + clean: @ echo 'use the clean script' diff --git a/phys/module_calc_lpi_new.F b/phys/module_calc_lpi_new.F new file mode 100644 index 0000000000..97e7630516 --- /dev/null +++ b/phys/module_calc_lpi_new.F @@ -0,0 +1,214 @@ +MODULE module_calc_lpi_new +CONTAINS + SUBROUTINE calc_lpi_new(qc, qr, qi, qs, qg & + ,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 + 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) + 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) +! 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 diff --git a/phys/module_ltng_pe.F b/phys/module_ltng_pe.F new file mode 100644 index 0000000000..ee7fce1ae9 --- /dev/null +++ b/phys/module_ltng_pe.F @@ -0,0 +1,329 @@ +MODULE module_ltng_pe +CONTAINS +!https://journals.ametsoc.org/doi/10.1175/WAF-D-11-00144.1 +! coul_pos, etc, are charging constants, in units of Coulombs.. +! coul_pos, etc, now set in Registry.EM_COMMON; can be added to namelist.input +! The values set work well for the NSSL scheme, and probably for Thompson +! Aerosol Aware. They are suggested values and can be modified by the user +! based on a set of case study comparisons. +! For the WSM6, it is suggested to double the values. + + SUBROUTINE calc_dyn_pe(dx,dy,qc, qr, qi, qs, qg & + ,w,dz8w,pi_phy,th_phy,p_phy,rho_phy & + ,light_pe & + ,light_ne & + ,light_neu & + ,dt & + ,coul_pos,coul_neg,coul_neu & + ,ids,ide, jds,jde, kds,kde & + ,ims,ime, jms,jme, kms,kme & + ,its,ite, jts,jte, kts,kte & + ) + + IMPLICIT NONE +!------------------------------------------------------------------- +! +! + REAL, INTENT(IN) :: coul_pos,coul_neg,coul_neu + 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 + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & + INTENT(INOUT),OPTIONAL :: & + light_pe,& + light_ne, & + light_neu + + REAL,INTENT(IN):: dx,dy + 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):: & + & plpi,nlpi,neulpi +! 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 + 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 :: dt,dxdy,t_base,t_top + INTEGER I_COLLAPSE + LOGICAL LOOK_T + real t_pos,t_neg,t_neu + 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 + t_pos=0.12 + t_neg=0.06 + t_neu=0.03 + dxdy=dx*dy +! 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)) + 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 + if (k.lt.kte)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) + qs1d(k)=qs(i,k,j) +! qg1d(k)=qg(i,k,j)+qh(i,k,j) + qg1d(k)=qg(i,k,j) +! qtot(k)=qc(i,k,j)+qr(i,k,j)+qi(i,k,j)+qs(i,k,j)+qg(i,k,j)+qh(i,k,j) + qtot(k)=qc(i,k,j)+qr(i,k,j)+qi(i,k,j)+qs(i,k,j)+qg(i,k,j) + qitot(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j) +! qitot(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j)+qh(i,k,j) +! For conservative advection multiply by rho1d and divide by it below + qitot(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j) + ENDDO + do k = kts,kte + height_t(k)=0.5*(height(k-1)+height(k)) + w1d_t(k)=0.5*(w1d(k-1)+w1d(k)) + end do + t_base=-10 + t_top=-30 + do k = kts,kte + e1d(k)=light_pe(i,k,j)*rho1d(k) + end do + call calc_lpi(ql1d,qi1d,qs1d,qg1d,w1d,temp,height,plpi(i,j),lpi1d,t_base,t_top,kme,kte) + look_t=.true. + do k=kte,kts,-1 + if (look_t)then + if (temp(k).ge.-27.and.qtot(k).gt.0.00001)then + t_pos=height_t(k)/60000. + look_t=.false. + end if + end if + end do + + call power_index(temp,qitot,rho1d,plpi(i,j),lpi1d,e1d,dt,dxdy,t_pos,t_base,t_top,coul_pos,kme,kte) + do k = kts,kte + light_pe(i,k,j)=e1d(k) + end do + do k = kts,kte + e1d(k)=light_ne(i,k,j)*rho1d(k) +! e1d(k)=0. +! IF (ISNAN(e1d(k))) STOP 'BARRY WANTS TO STOP HERE 18' + end do + t_base=-0 + t_top=-20 + call calc_lpi(ql1d,qi1d,qs1d,qg1d,w1d,temp,height,nlpi(i,j),lpi1d,t_base,t_top,kme,kte) +! OLD 2 + look_t=.true. + do k=kte,kts,-1 + if (look_t)then + if (temp(k).ge.-2.and.qtot(k).gt.0.0001)then + t_neg=height_t(k)/60000. + look_t=.false. + end if + end if + end do + call power_index(temp,qitot,rho1d,nlpi(i,j),lpi1d,e1d,dt,dxdy,t_neg,t_base,t_top,coul_neg,kme,kte) + do k = kts,kte + light_ne(i,k,j)=e1d(k) + end do +! OLD 3 + do k = kts,kte + e1d(k)=light_neu(i,k,j)*rho1d(k) + end do + t_base=0. + t_top=-30. + call calc_lpi(ql1d,qi1d,qs1d,qg1d,w1d,temp,height,neulpi(i,j),lpi1d,t_base,t_top,kme,kte) +! if(t_pos.ne.0.12)then +! print*,'t_pos = ',t_pos +! end if +! print*,'t_neg = ',t_neg + t_neu=0.5*(t_pos-t_neg) +! print*,'t_neu = ',t_neu + if (abs(t_pos-t_neg).le.0.01)t_neu=0.01 + call power_index(temp,qitot,rho1d,neulpi(i,j),lpi1d,e1d,dt,dxdy,t_neu,t_base,t_top,coul_neu,kme,kte) + do k = kts,kte + light_neu(i,k,j)=e1d(k) + end do + + + END DO + END DO +! print*,'after call' +! if (MAXVAL(light_neu).ne.0)print*,'maxval light_neu' +! if (MAXVAL(light_neu).ne.0)write(MAXVAL(light_neu)) +! 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)) + return + end subroutine calc_dyn_pe + subroutine & + & calc_lpi(ql3d,qi3d,qs3d,qg3d,w3d,t3d,height,lpi,lpi1d,t_base,t_top,nk,nke) + implicit none + integer nk,nke + real t_base,t_top + real lpi1d(nk) + real ql3d(nk) + real qg3d(nk) + real qi3d(nk) + real qs3d(nk) + real w3d(0:nk) + real t3d(nk) + real height(0:nk) + real lpi,plpi + real del_z(nk) + real w_ave(nk) + integer ic,jc,icnt,i,j,k,i_collapse + real i_dist,j_dist,del_z_tot + real top, bot + real num,den,ave_z + real num_s,den_s + real num_i,den_i + real q_isg + icnt=0 + do k=1,nke + lpi1d(k)=0 + top=height(k) + bot=height(k-1) + del_z(k)=top-bot + w_ave(k)=0.5*(w3d(k)+w3d(k-1)) + end do + lpi1d(nk)=0 +! +! Check for collapsing cell + ave_z=0 + del_z_tot=0 + lpi=0 + do k=1,nke + if (t3d(k).le.t_base.and.t3d(k).gt.t_top)then ! set temp range + + den_i = qi3d(k)+qg3d(k) + den_s = qs3d(k)+qg3d(k) + if (qs3d(k).eq.0.or.qg3d(k).eq.0.)then !checks for zeroes + den_s=10000. + num_s = 0. + else + num_s = sqrt(qs3d(k)*qg3d(k)) + end if + if (qi3d(k).eq.0.or.qg3d(k).eq.0.)then ! checks for zeroes + den_i=10000. + num_i = 0. + else + num_i = sqrt(qi3d(k)*qg3d(k)) + end if + q_isg = qg3d(k)*(num_i/den_i+num_s/den_s) ! ice "fract"-content + + if (ql3d(k).eq.0.or.q_isg.eq.0)then + num=0 + den=10000. + else + num = sqrt(ql3d(k)*q_isg) + den = ql3d(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)*w_ave(k)**2 ! power index J/unit-mass + lpi1d(k)=del_z(k)*(2.*num/den)*w_ave(k)**2 ! power index J/unit-mass + end if + end if + end do +! + if (del_z_tot.eq.0)del_z_tot=100000 + lpi=ave_z/del_z_tot + +! + return + end subroutine calc_lpi + subroutine & + & power_index(t3d,qitot,rho1d,lpi,lpi1d,le1d,dt,dxdy,t_pos,t_base,t_top,coul_value,nk,nke) + implicit none + integer nk,nke + real coul_value + real t_pos,t_top,t_base + real lpi1d(nk),lpi,dxdy + real le1d(nk) + real rho1d(nk) + real ql3d(nk) + real qitot(nk) + real qg3d(nk) + real qi3d(nk) + real qs3d(nk) + real w3d(nk) + real t3d(nk) + real del_z(nk) + real w_ave(nk) + real power + integer ic,jc,icnt,i,j,k,i_collapse + real i_dist,j_dist,del_z_tot + real top, bot + real num,den,ave_z + real num_s,den_s + real num_i,den_i + real q_isg + real dt + lpi=0 + le1d(nk)=0 + do k=1,nke + if (t3d(k).le.t_base.and.t3d(k).gt.t_top)then ! set temp range + power=coul_value/t_pos*qitot(k)*rho1d(k)*lpi1d(k) + le1d(k)=le1d(k)+power*dt*dxdy + lpi = lpi+le1d(k) + end if + end do + return + end subroutine power_index + END MODULE module_ltng_pe diff --git a/phys/module_ltng_strokes.F b/phys/module_ltng_strokes.F new file mode 100644 index 0000000000..e0c8b288fe --- /dev/null +++ b/phys/module_ltng_strokes.F @@ -0,0 +1,212 @@ +MODULE module_ltng_strokes +!https://journals.ametsoc.org/doi/10.1175/WAF-D-11-00144.1 + +! j_pos,j_neg,j_neu:Energy per flash; threshold values set in Registry.EM_COMMON +! Values can be modified through namelist.input +! Sample values used: parameter (j_pos=5.E9,j_neg=1.E9,j_neu=1.E9) +! References for this and accompanying program module_calc_lpi.F +! Note, the calculation of lpi is different than in the published paper by Yair +! et al., JGR: +! https://scholar.google.co.il/citations?user=vzBSZmEAAAAJ&hl=en#d=gs_md_cita-d&u=%2Fcitations%3Fview_op%3Dview_citation%26hl%3Den%26user%3DvzBSZmEAAAAJ%26citation_for_view%3DvzBSZmEAAAAJ%3AULOm3_A8WrAC%26tzom%3D-120 +! Lynn, B. H., Y. Yair, C. Price, G. Kelman, A. Clark, +! 2012: Predicting Cloud-to-Ground and Intracloud Lightning in Weather +!Forecast Models. Wea. Forecast. 27, 1470-1486. +! Lynn, B. H., G. Kelman, and G. Ellrod, 2015: An Evaluation of the +! Efficacy of Using Observed Lightning to Improve Convective Lightning Forecasts. +! Wea. Forecasting, 30, 405–423. +! Lynn, B.H, 2016: The Usefulness and Economic Value of Total Lightning +! Forecasts made with a "Dynamic" Lightning Scheme coupled with Lightning Data +! Assimilation. Wea. Forecasting. 32, 645 – 663. +! Contact: barry.lynn@weather-it-is.com or barry.h.lynn@gmail.com + +CONTAINS +!=================================================================== +! + SUBROUTINE calcstroke(dx,dy & + ,w,dz8w,pi_phy,th_phy,p_phy,rho_phy & + ,light_pe & + ,light_ne & + ,light_neu & + ,lpos,lneg,lneu& + ,j_pos,j_neg,j_neu & + ,ids,ide, jds,jde, kds,kde & + ,ims,ime, jms,jme, kms,kme & + ,its,ite, jts,jte, kts,kte & + ) +!=================================================================== +! + IMPLICIT NONE +!------------------------------------------------------------------- +! +! + REAL, INTENT(IN):: j_pos,j_neg,j_neu + real j_sa ! scale aware + 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(INOUT),OPTIONAL :: & + light_pe,& + light_ne, & + light_neu + REAL,INTENT(IN):: dx,dy + 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):: & + & LPOS,LNEG,LNEU +! REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme):: & +! & LPI,LPOS,LNEG,LNEU,plpi,nlpi,neulpi + + + + REAL, DIMENSION(kms:kme):: tempk,rh,qtot,qitot + REAL, DIMENSION(kms:kme):: qv1d,p1d,rho1d,qti1d + REAL, DIMENSION(kms:kme):: temp,qc1d,ql1d,qi1d,qs1d,qg1d,lpi1d + REAL, DIMENSION(0:kme):: w1d,height + REAL, DIMENSION(kms:kme):: e1d,height_t,w1d_t + REAL z_full + INTEGER level + INTEGER I_COLLAPSE + LOGICAL LOOK_T + real dxkm,dykm + real check_light_b,check_light_a + real check_light_bb,check_light_aa + real max,sum + + INTEGER :: i,j,k +!------------------------------------------------------------------- + +! NOTE: dxkm is a "scale-aware" term. +! The volume change from grid to grid depends on the square of the area. +! The vertical velocity, though, will increase, for example, by deltax +! the ratio of w to volume is deltax/deltax-sq or 1./delta_x +! 4 km is the reference value for deltax (based on paper published in 2012). +! print*,'j_pos = ',j_pos +! print*,'j_neg = ',j_neg +! print*,'j_neu = ',j_neu +! print*,'dx = ',dx + dxkm = dx*1.E-3 + dykm = dy*1.E-3 + DO j = jts,jte + DO i = its,ite + z_full=0. + height(0)=z_full + do k=kts,kte + z_full=z_full+dz8w(i,k,j) + height(k)=z_full + temp(k) = th_phy(i,k,j)*pi_phy(i,k,j)-273.16 + rho1d(k)=rho_phy(i,k,j) + end do +! note: light_neu, etc, are multiplied by rho in calc_lpi.F, for converting from +! /kg to /m^3 + e1d(kme) = 0 + do k=kts,kte + e1d(k)=light_neu(i,k,j) + end do +! j_neu=1.E9 + j_sa=j_neu*(dxkm/4.)*(dykm/4.) +! j_sa=j_neu*(dxkm/4.) +! check_light_b = maxval(e1d) +! check_light_bb = sum(e1d) + call flash(height,temp,e1d,j_sa,lneu(i,j),kme,kte) +! check_light_a = maxval(e1d) +! check_light_aa = sum(e1d) +! if (check_light_a.lt.check_light_b)then +! print*,'max i,k,j = ',i,k,j,check_light_a,check_light_b +! end if +! if (check_light_aa.lt.check_light_bb)then +! print*,'sum,i,k,j = ',i,k,j,check_light_aa,check_light_bb +! end if + do k=kts,kte + light_neu(i,k,j) = e1d(k)/rho1d(k) + end do + light_neu(i,kme,j) = 0 + e1d(kme) = 0 + do k = kts,kte + e1d(k)=light_ne(i,k,j) + end do + j_sa=j_neg*(dxkm/4.)*(dykm/4.) +! j_sa=j_neg*(dxkm/4.) + call flash(height,temp,e1d,j_sa,lneg(i,j),kme,kte) + do k = kts,kte + light_ne(i,k,j) = e1d(k)/rho1d(k) + end do + light_ne(i,kme,j) = 0 +! j_pos=5.E9 + j_sa=j_pos*(dxkm/4.)*(dykm/4.) +! j_sa=j_pos*(dxkm/4.) + e1d(kme) = 0 + do k = kts,kte + e1d(k)=light_pe(i,k,j) + end do + call flash(height,temp,e1d,j_sa,lpos(i,j),kme,kte) + do k = kts,kte + light_pe(i,k,j)=e1d(k)/rho1d(k) + end do + light_pe(i,kme,j) = 0 + END DO + END DO + + + return + end subroutine calcstroke + subroutine & + & flash(height,t3d,pe1d,j_pos,lpos,nk,nke) + implicit none + integer nk,nke + real lpos,j_pos + real pe1d(nk) + integer i,j,k + integer num_ran + parameter (num_ran=2) + real ran_num(num_ran) + real i_seed + real t3d(nk) + real height(0:nk) + real del_z(nk) + real top, bot + real n,ave_z,del_z_tot + real residual,e_p,e_n,e_neu,loss,e_min + del_z_tot=0 + do k=1,nke-1 + top=0.5*(height(k+1)+height(k)) + bot=0.5*(height(k)+height(k-1)) + del_z(k)=top-bot + end do +! plpi=0 + e_p=0 + do k=1,nke + e_p=pe1d(k)+e_p + end do +! FLASHES AND LOSSES +! POSITIVE + pe1d(1)=0 + pe1d(nk)=0 + pe1d(nke)=0 + do k=2,nke-1 + del_z_tot=del_z_tot+del_z(k) + end do + if (e_p.ge.j_pos)then + lpos=lpos+int(e_p/j_pos) + loss = int(e_p/j_pos) + residual=e_p-j_pos*loss + do k=2,nke-1 + e_min=min(pe1d(k),e_p) + pe1d(k)=pe1d(k)-e_min+residual*del_z(k)/del_z_tot + end do + end if + return + end subroutine flash +END MODULE module_ltng_strokes diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 94004191cb..d6a76c4207 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -16,7 +16,7 @@ SUBROUTINE microphysics_driver( & ,chem_opt, progn & ,cldfra, cldfra_old, exch_h, nsource & ,qlsink, precr, preci, precs, precg & - ,xland,snowh,xice,itimestep & + ,xland,snowh,itimestep & ,f_ice_phy,f_rain_phy,f_rimef_phy & ,lowlyr,sr, id & ,ids,ide, jds,jde, kds,kde & @@ -111,7 +111,6 @@ SUBROUTINE microphysics_driver( & #endif ,qnwfa2d, qnifa2d, qnbca2d & ! for water/ice-friendly/black carbon aerosols ,qnocbb2d, qnbcbb2d & ! for biomass burning aerosols - ,ssat,ssati & ,refl_10cm & ! HM, 9/22/09, add for refl ,vmi3d & ! for P3 ,di3d & ! for P3 @@ -162,16 +161,27 @@ SUBROUTINE microphysics_driver( & ,perts_qsnow, perts_ni & ,pert_thom_qv,pert_thom_qc,pert_thom_qi & ,pert_thom_qs,pert_thom_ni & - ,cloudnc & +!Dynamic Lightning + ,curr_secs & + ,light_pe_curr & + ,light_ne_curr & + ,light_neu_curr & + ,lpi2d & + ,lpi3d & + ,f_light_pe & + ,f_light_ne & + ,f_light_neu & + ,lpos,lneg,lneu & + ,l_obs,num_light_periods & ) ! Framework USE module_state_description, ONLY : & KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME & - ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, RCON_MP_SCHEME, THOMPSONGH, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & + ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, THOMPSONGH, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & ,GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, NSSL_2MOM, MADWRF_MP & ,FER_MP_HIRES_ADVECT & - ,WSM7SCHEME, WDM7SCHEME, UDMSCHEME & + ,WSM7SCHEME, WDM7SCHEME & ,NUWRF4ICESCHEME & ,MILBRANDT2MOM , CAMMGMPSCHEME,FULL_KHAIN_LYNN, P3_1CATEGORY, P3_1CATEGORY_NC, P3_2CATEGORY, P3_1CAT_3MOM & ,MORR_TM_AERO, JENSEN_ISHMAEL, SPRINKLER, NTU !,MILBRANDT3MOM @@ -207,6 +217,9 @@ SUBROUTINE microphysics_driver( & #endif ! *** add new modules of schemes here + USE module_calc_lpi_new + USE module_ltng_pe + USE module_ltng_strokes USE module_mp_kessler #if ( WRFPLUS == 1 ) @@ -222,7 +235,6 @@ SUBROUTINE microphysics_driver( & USE module_mp_wsm7 USE module_mp_etanew USE module_mp_fer_hires - USE module_mp_rcon USE module_mp_thompson USE module_mp_full_sbm #if ( BUILD_SBM_FAST == 1 ) @@ -240,7 +252,6 @@ SUBROUTINE microphysics_driver( & USE module_mp_wdm5 USE module_mp_wdm6 USE module_mp_wdm7 - USE module_mp_udm USE module_mp_milbrandt2mom # if (EM_CORE == 1) USE module_mp_cammgmp_driver, ONLY: CAMMGMP ! CAM5's microphysics driver @@ -599,8 +610,16 @@ SUBROUTINE microphysics_driver( & ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: XLAND +!Lightning + REAL , INTENT(IN ) :: curr_secs + REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT),OPTIONAL :: LPOS,LNEG,LNEU + INTEGER, INTENT(IN ):: num_light_periods + REAL, DIMENSION( ims:ime , 1:num_light_periods , jms:jme ), & + INTENT(IN) :: l_obs + REAL, DIMENSION( ims:ime , jms:jme ),INTENT(INOUT) :: lpi2d + REAL, DIMENSION( ims:ime ,kms:kme, jms:jme ),INTENT(INOUT) :: lpi3d + REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN), OPTIONAL :: SNOWH - REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: XICE REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: SR @@ -611,7 +630,7 @@ SUBROUTINE microphysics_driver( & ! ! Optional ! - REAL, OPTIONAL, DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(OUT) :: refl_10cm,ssat,ssati + REAL, OPTIONAL, DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(OUT) :: refl_10cm REAL, OPTIONAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: & ! for ntu3m qdcn_curr,qtcn_curr,qccn_curr,qrcn_curr,qnin_curr, & ! for ntu3m fi_curr,fs_curr,vi_curr,vs_curr,vg_curr,ai_curr, & ! for ntu3m @@ -653,7 +672,8 @@ SUBROUTINE microphysics_driver( & ,kext_ft_qic,kext_ft_qip,kext_ft_qid & ,kext_ft_qs,kext_ft_qg & ,qnwfa_curr,qnifa_curr,qnbca_curr & ! Added by G. Thompson - ,qvolg_curr,qvolh_curr, qrimef_curr + ,qvolg_curr,qvolh_curr, qrimef_curr & + ,light_pe_curr,light_neu_curr,light_ne_curr REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & @@ -689,7 +709,6 @@ SUBROUTINE microphysics_driver( & ,GRAUPELNCV & ,HAILNC & ,HAILNCV & - ,CLOUDNC & ,hail_maxk1, hail_max2d #if ( WRF_CHEM == 1) @@ -730,7 +749,8 @@ SUBROUTINE microphysics_driver( & ,f_qvoli,f_qaoli & ! for Jensen ISHMAEL ,f_qvoli2,f_qaoli2 & ! for Jensen ISHMAEL ,f_qi3,f_qni3,f_qvoli3,f_qaoli3 & ! for Jensen ISHMAEL - ,f_qnwfa, f_qnifa, f_qnbca ! Added by G. Thompson + ,f_qnwfa, f_qnifa, f_qnbca & ! Added by G. Thompson + ,f_light_pe,f_light_ne,f_light_neu LOGICAL, OPTIONAL, INTENT(IN) :: diagflag @@ -1031,90 +1051,6 @@ SUBROUTINE microphysics_driver( & CALL wrf_error_fatal ( 'arguments not present for calling mkessler' ) ENDIF #endif -! -! - CASE (RCON_MP_SCHEME) - - CALL wrf_debug ( 100 , 'microphysics_driver: calling RCON' ) - IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. & - PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. & - PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. & - PRESENT( QNR_CURR) .AND. PRESENT ( QNI_CURR) .AND. & - PRESENT( QNC_CURR) .AND. PRESENT ( QNWFA_CURR) .AND. & - PRESENT( QNIFA_CURR).AND.PRESENT ( QNWFA2D) .AND. & - PRESENT( QNIFA2D) .AND. & - PRESENT( SNOWNC) .AND. PRESENT ( SNOWNCV) .AND. & - PRESENT( GRAUPELNC).AND. PRESENT ( GRAUPELNCV) .AND. & - PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN -#if ( WRF_CHEM == 1 ) - qv_b4mp(its:ite,kts:kte,jts:jte) = qv_curr(its:ite,kts:kte,jts:jte) - qc_b4mp(its:ite,kts:kte,jts:jte) = qc_curr(its:ite,kts:kte,jts:jte) - qi_b4mp(its:ite,kts:kte,jts:jte) = qi_curr(its:ite,kts:kte,jts:jte) - qs_b4mp(its:ite,kts:kte,jts:jte) = qs_curr(its:ite,kts:kte,jts:jte) -#endif - CALL mp_rcon_driver( & - QV=qv_curr, & - QC=qc_curr, & - QR=qr_curr, & - QI=qi_curr, & - QS=qs_curr, & - QG=qg_curr, & - NI=qni_curr, & - NR=qnr_curr, & - NC=qnc_curr, & - NWFA=qnwfa_curr, & - NIFA=qnifa_curr, & - NBCA=qnbca_curr, & - NWFA2D=qnwfa2d, & - NIFA2D=qnifa2d, & - NBCA2D=qnbca2d, & - aer_init_opt=config_flags%aer_init_opt, & - wif_input_opt=config_flags%wif_input_opt, & - TH=th, & - PII=pi_phy, & - P=p, & - W=w, & - DZ=dz8w, & - DT_IN=dt, & - ITIMESTEP=itimestep, & - RAINNC=RAINNC, & - RAINNCV=RAINNCV, & - CLOUDNC=CLOUDNC, & - SNOWNC=SNOWNC, & - SNOWNCV=SNOWNCV, & - GRAUPELNC=GRAUPELNC, & - GRAUPELNCV=GRAUPELNCV, & - SR=SR, & -#if ( WRF_CHEM == 1 ) - WETSCAV_ON=config_flags%wetscav_onoff == 1, & - RAINPROD=rainprod, & - EVAPPROD=evapprod, & -#endif - REFL_10CM=refl_10cm, & - diagflag=diagflag, & - ke_diag = ke_diag, & - do_radar_ref=do_radar_ref, & - re_cloud=re_cloud, & - re_ice=re_ice, & - re_snow=re_snow, & - has_reqc=has_reqc, & - has_reqi=has_reqi, & - has_reqs=has_reqs, & - IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, & - IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, & - ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte) - - - IF (config_flags%aer_fire_emit_opt.gt.0) then - CALL wrf_debug ( 200 , ' call fire_emis_simple_plumerise' ) - CALL fire_emis_simple_plumerise (config_flags%wif_fire_inj, config_flags%aer_fire_emit_opt & - ,z_at_mass, pblh, qnwfa_curr, qnbca_curr & - ,qnocbb2d, qnbcbb2d, dt, ids, ide, jds, jde, kds, kde & - ,ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte ) - ENDIF - ELSE - CALL wrf_error_fatal ( 'arguments not present for calling rcon' ) - ENDIF ! CASE (THOMPSONAERO) if (pert_thom .and. multi_perturb == 1) then @@ -2078,9 +2014,6 @@ SUBROUTINE microphysics_driver( & GRPLNCV = GRAUPELNCV, & SR=SR, & dbz = refl_10cm, & - ssat3d = ssat, & - ssati = ssati, & - nssl_ssat_output = config_flags%nssl_ssat_output, & #if ( WRF_CHEM == 1 ) WETSCAV_ON = config_flags%wetscav_onoff == 1, & EVAPPROD=evapprod,RAINPROD=rainprod, & @@ -2659,58 +2592,6 @@ SUBROUTINE microphysics_driver( & CALL wrf_error_fatal ( 'arguments not present for calling wdm7') ENDIF - CASE (UDMSCHEME) - CALL wrf_debug ( 100 , 'microphysics_driver: calling udm' ) - IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. & - PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. & - PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. & - PRESENT( QH_CURR ) .AND. & - PRESENT( QNN_CURR ) .AND. PRESENT ( QNC_CURR ) .AND. & - PRESENT( QNR_CURR ).AND. & - PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN - CALL udm( & - TH=th & - ,Q=qv_curr & - ,QC=qc_curr & - ,QR=qr_curr & - ,QI=qi_curr & - ,QS=qs_curr & - ,QG=qg_curr & - ,QH=qh_curr & - ,NN=qnn_curr & - ,NC=qnc_curr & - ,NR=qnr_curr & - ,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w & - ,DELT=dt,G=g,CPD=cp,CPV=cpv,CCN0=ccn_conc & ! RAS - ,RD=r_d,RV=r_v,T0C=svpt0 & - ,EP1=ep_1, EP2=ep_2, QMIN=epsilon & - ,XLS=xls, XLV0=xlv, XLF0=xlf & - ,DEN0=rhoair0, DENR=rhowater & - ,CLIQ=cliq,CICE=cice,PSAT=psat & - ,xland=xland, xice=xice & ! land mask, 1: land, 2: water - ,RAIN=rainnc ,RAINNCV=rainncv & - ,SNOW=snownc ,SNOWNCV=snowncv & - ,SR=sr & - ,REFL_10CM=refl_10cm & ! added for radar reflectivity - ,diagflag=diagflag & ! added for radar reflectivity - ,do_radar_ref=do_radar_ref & ! added for radar reflectivity - ,GRAUPEL=graupelnc ,GRAUPELNCV=graupelncv & - ,HAIL=hailnc ,HAILNCV=hailncv & - ,ITIMESTEP=itimestep & - ,has_reqc=has_reqc & ! for radiation + - ,has_reqi=has_reqi & - ,has_reqs=has_reqs & - ,re_cloud=re_cloud & - ,re_ice=re_ice & - ,re_snow=re_snow & ! for radiation - - ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & - ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & - ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & - ) - ELSE - CALL wrf_error_fatal ( 'arguments not present for calling udm') - ENDIF - CASE (ETAMPNEW) !-- Operational 4-km High-Resolution Window (HRW) version CALL wrf_debug ( 100 , 'microphysics_driver: calling etampnew') @@ -2901,7 +2782,68 @@ SUBROUTINE microphysics_driver( & CALL wrf_error_fatal ( wrf_err_message ) END SELECT micro_select +#if (EM_CORE==1) + IF ( PRESENT ( QC_CURR ) .AND. & + PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. & + PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ))Then +! print*,'config_flags%dyn_lightning_option = ',config_flags%dyn_lightning_option + if (config_flags%dyn_lightning_option.eq.1)then + CALL calc_lpi_new( & + W=w, & + PI_PHY=pi_phy, RHO_PHY=rho, & + TH_PHY=TH,P_PHY=p, & + DZ8w=dz8w, & + QC=qc_curr, & + QR=qr_curr, & + QI=qi_curr, & + QS=qs_curr, & + QG=qg_curr, & + lpi2d = lpi2d, & + lpi3d = lpi3d & + ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & + ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & + ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte) + CALL calc_dyn_pe(dx=dx,dy=dy, & + W=w, & + PI_PHY=pi_phy, RHO_PHY=rho, & + TH_PHY=TH,P_PHY=p, & + DZ8w=dz8w, & + QC=qc_curr, & + QR=qr_curr, & + QI=qi_curr, & + QS=qs_curr, & + QG=qg_curr, & + light_pe=light_pe_curr, & + light_ne=light_ne_curr, & + light_neu=light_neu_curr, & + dt=dt & + ,coul_pos=config_flags%coul_pos & + ,coul_neg=config_flags%coul_neg & + ,coul_neu=config_flags%coul_neu & + ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & + ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & + ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte) + CALL calcstroke(dx=dx,dy=dy,& + W=w, & + PI_PHY=pi_phy, RHO_PHY=rho, & + TH_PHY=TH,P_PHY=p, & + DZ8w=dz8w, & + light_pe=light_pe_curr, & + light_ne=light_ne_curr, & + light_neu=light_neu_curr, & + lpos=lpos,lneg=lneg,lneu=lneu & + ,j_pos=config_flags%j_pos & + ,j_neg=config_flags%j_neg & + ,j_neu=config_flags%j_neu & + ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & + ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & + ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte) + END IF + ELSE + CALL wrf_error_fatal ( 'arguments not present for calling dynamic lightning scheme' ) + ENDIF +#endif ENDDO !$OMP END PARALLEL DO