phys/module_ltng_lpi.F

   1 MODULE module_ltng_lpi
   2 !Yair, Y., B. Lynn, C. Price, V. Kotroni, K. Lagouvardos, E. Morin,
   3 !A. Magnai, and M. del Carmen Llasat (2010), Predicting the potential for
   4 !lightning activity in Mediterranean storms based on the Weather
   5 !Research and Forecasting (WRF) model dynamic and microphysical
   6 !fields, J. Geophys. Res., 115, D04205, doi:10.1029/2008JD010868.
   7 ! However, we don't check for collapsing cell (so as not to require use of halo).
   8 ! This means that lpi is also calculated in cells that are no longer (on average) growing
   9 ! For a "complete" lightning forecast scheme, please see:
  10 !https://doi.org/10.1175/WAF-D-11-00144.1
  11 !(Predicting Cloud-to-Ground and Intracloud Lightning in Weather Forecast Models)
  12
  13 CONTAINS
  14 !===================================================================
  15 !
  16   SUBROUTINE calclpi(qv,qc, qr, qi, qs, qg                         &
  17                  ,w,z,dz8w,pi_phy,th_phy,p_phy,rho_phy             &
  18                  ,lpi                                              &
  19                  ,ids,ide, jds,jde, kds,kde                        &
  20                  ,ims,ime, jms,jme, kms,kme                        &
  21                  ,its,ite, jts,jte, kts,kte                        &
  22                                                                    )
  23 !-------------------------------------------------------------------
  24   IMPLICIT NONE
  25 !-------------------------------------------------------------------
  26 !
  27 !
  28   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
  29                                       ims,ime, jms,jme, kms,kme , &
  30                                       its,ite, jts,jte, kts,kte
  31   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
  32         INTENT(IN) ::                                          &
  33                                                               qv, &
  34                                                               qc, &
  35                                                               qi, &
  36                                                               qr, &
  37                                                               qs, &
  38                                                               qg
  39
  40       REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
  41          INTENT(IN ) ::  w, z
  42       REAL, INTENT(IN),     DIMENSION(ims:ime, kms:kme, jms:jme)::      &
  43      &                      dz8w,pi_phy,p_phy,rho_phy
  44       REAL, INTENT(IN),  DIMENSION(ims:ime, kms:kme, jms:jme)::      &
  45      &                      th_phy
  46       REAL, INTENT(INOUT),  DIMENSION(ims:ime,jms:jme)::      &
  47      &                      LPI
  48
  49
  50
  51
  52       REAL, DIMENSION(kms:kme)::    tempk,rh
  53       REAL, DIMENSION(kms:kme):: qv1d,p1d,rho1d,qti1d
  54       REAL, DIMENSION(kms:kme):: temp,qc1d,ql1d,qi1d,qs1d,qg1d,lpi1d
  55       REAL, DIMENSION(0:kme):: w1d,height
  56       REAL, DIMENSION(kms:kme):: e1d,height_t,w1d_t
  57       REAL z_full,qrs,teten,RELHUM,LOC,Td_850,Td_700,PC_DWPT
  58       INTEGER level
  59       REAL :: dt_lpi,t_base,t_top
  60       INTEGER I_COLLAPSE
  61       LOGICAL LOOK_T
  62       INTEGER I_START,I_END,J_START,J_END
  63
  64
  65   INTEGER ::               i,j,k
  66 !-------------------------------------------------------------------
  67       DO j = jts,jte
  68       DO i = its,ite
  69         z_full=0.
  70         height(0)=z_full
  71         w1d(0)=w(i,1,j)
  72       DO k = kts,kte-1
  73           if (k.lt.kte-1)then
  74            w1d(k)=w(i,k+1,j)
  75           else
  76            w1d(k)=0.
  77           end if
  78           temp(k) = th_phy(i,k,j)*pi_phy(i,k,j)-273.16
  79           tempk(k) = th_phy(i,k,j)*pi_phy(i,k,j)
  80           qv1d(k)=qv(i,k,j)
  81           p1d(k)=p_phy(i,k,j)
  82           rho1d(k)=rho_phy(i,k,j)
  83           z_full=z_full+dz8w(i,k,j)
  84           height(k)=z_full
  85           qc1d(k)=qc(i,k,j)
  86           ql1d(k)=qc(i,k,j)+qr(i,k,j)
  87           qi1d(k)=qi(i,k,j)
  88           qti1d(k)=qi(i,k,j)+qs(i,k,j)+qg(i,k,j)
  89           qs1d(k)=qs(i,k,j)
  90 !         qg1d(k)=qg(i,k,j)+qh(i,k,j)
  91 ! Hail doesn't usually charge
  92           qg1d(k)=qg(i,k,j)
  93 ! For conservative advection multiply by rho1d and divide by it below
  94       ENDDO
  95       do k = kts,kte-1
  96        height_t(k)=0.5*(height(k-1)+height(k))
  97        w1d_t(k)=0.5*(w1d(k-1)+w1d(k))
  98       end do
  99       t_base=-0
 100       t_top=-20
 101       call calc_lpi(ql1d,qi1d,qs1d,qg1d,w1d,temp,height,lpi(i,j),t_base,t_top,kme,kte)
 102       END DO
 103       END DO
 104       return
 105       end subroutine calclpi
 106
 107       subroutine &
 108      &  calc_lpi(ql3d,qi3d,qs3d,qg3d,w3d,t3d,height,lpi,t_base,t_top,nk,nke)
 109       implicit none
 110       integer nk,nke
 111       real t_base,t_top
 112       real ql3d(nk)
 113       real qg3d(nk)
 114       real qi3d(nk)
 115       real qs3d(nk)
 116       real w3d(0:nk)
 117       real t3d(nk)
 118       real height(0:nk)
 119       real lpi
 120       real del_z(nk)
 121       real w_ave(nk)
 122       integer ic,jc,icnt,i,j,k,i_collapse
 123       real i_dist,j_dist,del_z_tot
 124       real top, bot
 125       real num,den,ave_z
 126       real num_s,den_s
 127       real num_i,den_i
 128       real q_isg
 129       real small_num
 130
 131       small_num=1.e-10
 132       icnt=0
 133       do k=1,nke
 134         top=height(k)
 135         bot=height(k-1)
 136         del_z(k)=top-bot
 137         w_ave(k)=0.5*(w3d(k)+w3d(k-1))
 138       end do
 139 !
 140 !     Check for collapsing cell
 141 ! Here, we don't check, since it requires a halo.
 142       ave_z=0
 143       del_z_tot=0
 144       lpi=0
 145
 146       do k=1,nke-1
 147        if (t3d(k).le.t_base.and.t3d(k).gt.t_top)then ! set temp range
 148
 149         den_i = qi3d(k)+qg3d(k)
 150         den_s = qs3d(k)+qg3d(k)
 151         if (qs3d(k).le.small_num.or.qg3d(k).le.small_num)then !checks for zeroes
 152          den_s=10000.
 153          num_s = 0.
 154         else
 155          num_s = sqrt(qs3d(k)*qg3d(k))
 156         end if
 157         if (qi3d(k).le.small_num.or.qg3d(k).le.small_num)then  ! checks for zeroes
 158          den_i=10000.
 159          num_i = 0.
 160         else
 161          num_i = sqrt(qi3d(k)*qg3d(k))
 162         end if
 163         q_isg = qg3d(k)*(num_i/den_i+num_s/den_s)  ! ice "fract"-content
 164
 165         if (ql3d(k).le.small_num.or.q_isg.eq.0)then
 166           num=0
 167           den=10000.
 168         else
 169          num = sqrt(ql3d(k)*q_isg)
 170          den = ql3d(k)+q_isg
 171         end if
 172         del_z_tot=del_z_tot+del_z(k)
 173         if (num.gt.0)then
 174          ave_z=ave_z+del_z(k)*(2.*num/den)*w_ave(k)**2 ! lightning potential index J/unit-mass
 175         end if
 176        end if
 177       end do
 178 !
 179       if (del_z_tot.eq.0)del_z_tot=100000
 180       lpi=ave_z/del_z_tot
 181 !
 182       return
 183       end subroutine calc_lpi
 184   END MODULE module_ltng_lpi