SCALE-RM
scale_atmos_phy_lt_sato2019.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_phy_lt_sato2019
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use mpi
18  use scale_precision
19  use scale_io
20  use scale_prof
21  !-----------------------------------------------------------------------------
22  implicit none
23  private
24  !-----------------------------------------------------------------------------
25  !
26  !++ Public procedure
27  !
28  public :: atmos_phy_lt_sato2019_setup
29  public :: atmos_phy_lt_sato2019_adjustment
30  public :: atmos_phy_lt_sato2019_select_dqcrg_from_lut
31 ! public :: ATMOS_PHY_LT_sato2019_mkinit
32 ! public :: ATMOS_PHY_LT_electric_field
33 ! public :: ATMOS_PHY_LT_neutralization_MG2001
34 ! public :: ATMOS_PHY_LT_neutralization_F2013
35 ! public :: ATMOS_PHY_LT_judge_absE
36  !-----------------------------------------------------------------------------
37  !
38  !++ Public parameters & variables
39  !
40  !-----------------------------------------------------------------------------
41  !
42  !++ Private procedure
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private parameters & variables
47  !
48  !--- for Bi-CGSTAB to solve the poisson equation
49  integer, private :: ITMAX = 3000
50  real(RP), private :: epsilon = 0.1_rp ** (rp*2)
51  character(len=64),private :: NUTR_TYPE='F2013'
52  character(len=64),private,save :: NUTR_qhyd='Each_POLARITY2'
53  integer, private :: MAX_NUTR = 100
54  real(RP), private :: Eint = 150.0e+3_rp ! [V/m]
55  real(RP), private :: delEint = 10.0e+3_rp ! [V/m]
56  real(RP), private :: Estop = 15.0e+3_rp ! [V/m]
57  real(RP), private :: qrho_chan = 0.5_rp ! [nC/m3]
58  real(RP), private :: qrho_neut = 0.5_rp ! [nC/m3]
59  real(RP), private :: fp = 0.3_rp ! [-]
60  real(RP), private :: zcg = 0.0_rp ! lowest grid point of lightning path for MG2001 [m]
61  integer, private :: NUTR_ITMAX = 1000
62  integer, private :: KIJMAXG
63  character(len=H_LONG) :: ATMOS_PHY_LT_LUT_FILENAME !--- LUT file name
64  character(len=H_LONG) :: fname_lut_lt="LUT_TK1978_v.txt" !--- LUT file name
65  integer, save :: fid_lut_lt
66  real(RP), private :: R_neut = 2000.0_rp ! [m]
67  real(RP), private :: rho0 = 1.225_rp ! [kg/m3]
68  real(RP), save :: flg_eint_hgt
69  logical, private :: LT_DO_Lightning = .true.
70 
71  !--- for history output
72  real(RP), allocatable, private :: d_QCRG_TOT(:,:,:)
73  real(RP), allocatable, private :: LT_PATH_TOT(:,:,:,:)
74  real(RP), allocatable, private :: fls_int_p_tot(:,:,:)
75 
76  !---
77  integer, parameter, private :: nxlut_lt = 200, nylut_lt = 200
78  real(RP), private :: dq_chrg( nxlut_lt,nylut_lt ) !--- charge separation [fC]
79  real(RP), private :: grid_lut_t( nxlut_lt )
80  real(RP), private :: grid_lut_l( nylut_lt )
81  real(RP), private :: tcrglimit
82  logical, private :: Hgt_dependency_Eint = .false.
83 
84  !--- Index for local variable
85  integer, private, parameter :: I_lt_x = 1
86  integer, private, parameter :: I_lt_y = 2
87  integer, private, parameter :: I_lt_z = 3
88  integer, private, parameter :: I_lt_abs = 4
89 
90  !--- For history output
91  integer, private, parameter :: w_nmax = 10
92  integer, private, parameter :: I_Ex = 1
93  integer, private, parameter :: I_Ey = 2
94  integer, private, parameter :: I_Ez = 3
95  integer, private, parameter :: I_Eabs = 4
96  integer, private, parameter :: I_Epot = 5
97  integer, private, parameter :: I_Qneut = 6
98  integer, private, parameter :: I_LTpath = 7
99  integer, private, parameter :: I_PosFLASH = 8
100  integer, private, parameter :: I_NegFLASH = 9
101  integer, private, parameter :: I_FlashPoint = 10
102  integer, private :: HIST_id(w_nmax)
103  character(len=H_SHORT), private :: w_name(w_nmax)
104  character(len=H_MID), private :: w_longname(w_nmax)
105  character(len=H_SHORT), private :: w_unit(w_nmax)
106  data w_name / 'Ex', &
107  'Ey', &
108  'Ez', &
109  'Eabs', &
110  'Epot', &
111  'Qneut', &
112  'LTpath', &
113  'PosFLASH', &
114  'NegFLASH', &
115  'FlashPoint' /
116  data w_longname / &
117  'X component of Electrical Field', &
118  'Y component of Electrical Field', &
119  'Z component of Electrical Field', &
120  'Absolute value of Electrical Field', &
121  'Electric Potential', &
122  'Cumulative Neutralizated charge', &
123  'Cumulative Number of flash path', &
124  'Cumulative Number of Positive flash', &
125  'Cumulative Number of Negative flash', &
126  'Cumulative Number of Flash point' /
127  data w_unit / &
128  'kV/m', &
129  'kV/m', &
130  'kV/m', &
131  'kV/m', &
132  'V', &
133  'nC/m3', &
134  'num', &
135  'num', &
136  'num', &
137  'num' /
138  !-----------------------------------------------------------------------------
139 contains
140  !-----------------------------------------------------------------------------
142  subroutine atmos_phy_lt_sato2019_setup( KA, KS, KE, &
143  IA, IS, IE, &
144  JA, JS, JE, &
145  IMAXG, &
146  JMAXG, &
147  KMAX, &
148  MP_TYPE, &
149  CDX, CDY )
150  use scale_prc, only: &
151  prc_abort, &
152  prc_ismaster, &
153  prc_myrank
154  use scale_comm_cartesc, only: &
155  comm_bcast
156  use scale_atmos_grid_cartesc, only: &
157  cz => atmos_grid_cartesc_cz
158  use scale_const, only: &
159  t00 => const_tem00
160  use scale_file_history, only: &
161  file_history_reg
162  implicit none
163  integer, intent(in) :: ka, ks, ke
164  integer, intent(in) :: ia, is, ie
165  integer, intent(in) :: ja, js, je
166  integer, intent(in) :: imaxg
167  integer, intent(in) :: jmaxg
168  integer, intent(in) :: kmax
169 
170  character(len=*), intent(in) :: mp_type
171  real(rp), intent(in) :: cdx(ia)
172  real(rp), intent(in) :: cdy(ja)
173 
174  integer :: n, myu, ip
175  integer :: ierr
176 
177  namelist / param_atmos_phy_lt_sato2019 / &
178  nutr_type, &
179  atmos_phy_lt_lut_filename, &
180  itmax, &
181  epsilon, &
182  eint, &
183  deleint, &
184  estop, &
185  qrho_chan, &
186  qrho_neut, &
187  fp, &
188  nutr_itmax, &
189  max_nutr, &
190  nutr_qhyd, &
191  zcg, &
192  r_neut, &
193  hgt_dependency_eint, &
194  lt_do_lightning
195 
196  !--- read namelist
197  rewind(io_fid_conf)
198  read(io_fid_conf,nml=param_atmos_phy_lt_sato2019,iostat=ierr)
199 
200  if( ierr < 0 ) then !--- missing
201  log_info("ATMOS_PHY_LT_sato2019_setup",*) '*** Not found namelist. Default used.'
202  elseif( ierr > 0 ) then !--- fatal error
203  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_LT_SATO2019. Check!'
204  call prc_abort
205  endif
206  log_nml(param_atmos_phy_lt_sato2019)
207 
208  !--- If zcg is lower than lowest grid point, zcg set lowest grid point
209  if( zcg <= cz(ks) ) then
210  zcg = cz(ks)
211  endif
212 
213  if( r_neut <= minval( cdx,1 ) .or. r_neut <= minval( cdy,1 ) ) then
214  r_neut = 2.d0 * min( minval( cdx,1 ), minval( cdy,1 ) )
215  endif
216 
217  if ( prc_ismaster ) then
218  fname_lut_lt = atmos_phy_lt_lut_filename
219  fid_lut_lt = io_get_available_fid()
220  !--- open parameter of cloud microphysics
221  open ( fid_lut_lt, file = fname_lut_lt, form = 'formatted', status = 'old', iostat=ierr )
222 
223  if ( ierr == 0 ) then
224  log_info("ATMOS_PHY_LT_sato2019_setup",'(2A)') 'Read LUT of TK78 table from ', trim(fname_lut_lt)
225  read( fid_lut_lt,* )
226  do n = 1, nylut_lt
227  do myu = 1, nxlut_lt
228  read( fid_lut_lt,* ) grid_lut_t( myu ), grid_lut_l( n ), dq_chrg( myu,n )
229  enddo
230  enddo
231  !--- LUT file does not exist
232  else
233  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx LUT for LT is requied when ATMOS_PHY_LT_TYPE = SATO2019, stop!'
234  call prc_abort
235  endif
236 
237  if( nutr_type /= 'MG2001' .and. nutr_type /= 'F2013' ) then
238  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx NUTR_TYPE should be MG2001 or F2013 stop!'
239  call prc_abort
240  endif
241  endif
242 
243  call comm_bcast( dq_chrg, nxlut_lt,nylut_lt )
244  call comm_bcast( grid_lut_t,nxlut_lt )
245  call comm_bcast( grid_lut_l,nylut_lt )
246 
247 ! KIJMAXG = (IEG-ISG+1)*(JEG-JSG+1)*(KE-KS+1)
248  kijmaxg = imaxg*jmaxg*kmax
249 
250  ! for history output
251  allocate( d_qcrg_tot(ka,ia,ja) )
252  allocate( lt_path_tot(ka,ia,ja,3) )
253  allocate( fls_int_p_tot(ka,ia,ja) )
254  d_qcrg_tot(:,:,:) = 0.0_rp
255  lt_path_tot(:,:,:,:) = 0.0_rp
256  fls_int_p_tot(:,:,:) = 0.0_rp
257 
258  tcrglimit = -60.0_rp+t00
259 
260  if( nutr_type == 'F2013' ) then
261  log_info("ATMOS_PHY_LT_sato2019_setup",'(A,F15.7,A)') 'Radius of neutralization is ', r_neut, "[m]"
262  endif
263 
264  flg_eint_hgt = 0.0_rp
265  if( hgt_dependency_eint .and. nutr_type == 'MG2001') then
266  flg_eint_hgt = 1.0_rp
267  endif
268 
269  if( nutr_qhyd /= 'TOTAL' .and. nutr_qhyd /= 'Each_POLARITY' .and. nutr_qhyd /= 'Each_POLARITY2' ) then
270  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx NUTR_qhyd should be TOTAL, or Each_POLARITY, or Each_POLARITY2, stop!'
271  call prc_abort
272  endif
273 
274 ! if( ( NUTR_qhyd == 'Each_POLARITY' .or. NUTR_qhyd /= 'Each_POLARITY2' ) .and. MP_TYPE == 'SUZUKI10' ) then
275 ! LOG_ERROR("ATMOS_PHY_LT_sato2019_setup",*) 'NUTR_qhyd = Each_POLARITY, Each_POLARITY2 is not supported for MP_TYPE SUZUKI10'
276 ! call PRC_abort
277 ! endif
278 
279  do ip = 1, w_nmax
280  call file_history_reg( w_name(ip), w_longname(ip), w_unit(ip), & ! [IN]
281  hist_id(ip) ) ! [OUT]
282  end do
283 
284  return
285  end subroutine atmos_phy_lt_sato2019_setup
286 
287  !-----------------------------------------------------------------------------
289  subroutine atmos_phy_lt_sato2019_adjustment( &
290  KA, KS, KE, &
291  IA, IS, IE, &
292  JA, JS, JE, &
293  KIJMAX, &
294  IMAX, &
295  JMAX, &
296  QA_LT, &
297  DENS, &
298  RHOT, &
299  QHYD, &
300  Sarea, &
301  dt_LT, &
302  QTRC, &
303  Epot )
304  use scale_const, only: &
305  small => const_eps
306  use scale_prc, only: &
307  prc_ismaster, &
308  prc_abort
309  use scale_comm_cartesc, only: &
310  comm_wait, &
311  comm_vars8
312  use scale_file_history, only: &
313  file_history_query, &
314  file_history_put
315  implicit none
316 
317  integer, intent(in) :: ka, ks, ke
318  integer, intent(in) :: ia, is, ie
319  integer, intent(in) :: ja, js, je
320  integer, intent(in) :: kijmax
321  integer, intent(in) :: imax
322  integer, intent(in) :: jmax
323 ! character(len=H_SHORT)
324  integer, intent(in) :: qa_lt
325  real(rp), intent(in) :: dens(ka,ia,ja)
326  real(rp), intent(in) :: rhot(ka,ia,ja)
327  real(rp), intent(in) :: qhyd(ka,ia,ja)
328  real(rp), intent(in) :: sarea(ka,ia,ja,qa_lt) !--- Surface area and that of each catergory [m2]
329  real(dp), intent(in) :: dt_lt
330  real(rp), intent(inout) :: qtrc(ka,ia,ja,qa_lt)
331  real(rp), intent(inout) :: epot(ka,ia,ja) !--- Electrical potential at previous time step [V]
332 
333  real(rp) :: qhyd_mass(ka,ia,ja) !--- Mass of total hydrometeor [kg/m3]
334  real(rp) :: rhoq0 !--- Tracer after lightning component
335  real(rp) :: qcrg(ka,ia,ja) !--- Total charge density [nC/m3]
336  real(rp) :: efield(ka,ia,ja,i_lt_abs) !--- Electrical field (1-3 ->x,y,z, 4->abs. )
337  real(rp) :: num_end(ka,ia,ja,3) !--- Number of each flash type (1->negative, 2->ground, 3->positive)
338  real(rp) :: d_qcrg(ka,ia,ja) !--- Change of charge by charge neutralization [fC/m3]
339  real(rp) :: lt_path(ka,ia,ja) !--- Lightning path (0-> no path, 1-> path)
340  real(rp) :: fls_int_p(ka,ia,ja)
341  real(rp) :: total_sarea(2) !--- Sum of surface area and that of each catergory [m2]
342  real(rp) :: neg_crg, pos_crg
343  real(rp) :: frac, r_totalsarea(2)
344  real(rp) :: dqneut(ka,ia,ja,qa_lt)
345  real(rp) :: dqneut_real(ka,ia,ja,qa_lt)
346  real(rp) :: dqneut_real_tot(ka,ia,ja)
347  logical :: flg_chrged(qa_lt)
348  real(rp) :: emax, emax_old
349  logical :: output_step
350  integer :: flg_lt_neut
351  integer :: i, j, k, m, n, countbin, ip
352  real(rp) :: sw, zerosw, positive, negative
353  integer :: count_neut
354 
355  logical :: hist_sw(w_nmax)
356  real(rp) :: w3d(ka,ia,ja)
357 
358  real(rp) :: diff_qcrg(0:1), lack(0:1), sum_crg(0:1)
359  real(rp) :: crg_rate(qa_lt), qcrg_before(qa_lt)
360  integer :: int_sw
361 
362  num_end(:,:,:,:) = 0.0_rp
363 
364  !$omp parallel do &
365  !$omp private(RHOQ0)
366  do j = js, je
367  do i = is, ie
368 
369  ! calc total charge density
370  do k = ks, ke
371  qcrg(k,i,j) = 0.0_rp
372  do n = 1, qa_lt
373  qcrg(k,i,j) = qcrg(k,i,j) + qtrc(k,i,j,n)
374  end do
375  qcrg(k,i,j) = qcrg(k,i,j) * dens(k,i,j) * 1.e-6_rp ![fC/kg] -> [nc/m3]
376  enddo
377 
378  do k = ks, ke
379  qhyd_mass(k,i,j) = qhyd(k,i,j) * dens(k,i,j) ![kg/kg] -> [kg/m3]
380  enddo
381 
382  enddo
383  enddo
384 
385  !--- Calculate E field
386  call atmos_phy_lt_electric_field( ka, ks, ke, & ! [IN]
387  ia, is, ie, & ! [IN]
388  ja, js, je, & ! [IN]
389  qcrg(:,:,:), & ! [IN]
390  dens(:,:,:), & ! [IN]
391  rhot(:,:,:), & ! [IN]
392  epot(:,:,:), & ! [INOUT]
393  efield(:,:,:,i_lt_x:i_lt_z) ) ! [OUT]
394 
395  !$omp parallel do
396  do j = js, je
397  do i = is, ie
398  do k = ks, ke
399  efield(k,i,j,i_lt_abs) = sqrt( efield(k,i,j,i_lt_x)*efield(k,i,j,i_lt_x) &
400  + efield(k,i,j,i_lt_y)*efield(k,i,j,i_lt_y) &
401  + efield(k,i,j,i_lt_z)*efield(k,i,j,i_lt_z) )
402 
403 
404  lt_path(k,i,j) = 0.0_rp
405  enddo
406  enddo
407  enddo
408 
409 
410  if ( lt_do_lightning ) then
411 
412  call atmos_phy_lt_judge_abse( ka, ks, ke, & ! [IN]
413  ia, is, ie, & ! [IN]
414  ja, js, je, & ! [IN]
415  dens(:,:,:), & ! [IN]
416  efield(:,:,:,i_lt_abs), & ! [IN]
417  emax, & ! [OUT]
418  flg_lt_neut ) ! [OUT]
419 
420 
421  count_neut = 0
422  do while( flg_lt_neut > 0 )
423 
424  emax_old = emax
425 
426  call comm_vars8( efield(:,:,:,i_lt_x), 1 )
427  call comm_vars8( efield(:,:,:,i_lt_y), 2 )
428  call comm_vars8( efield(:,:,:,i_lt_z), 3 )
429  call comm_vars8( efield(:,:,:,i_lt_abs), 4 )
430  call comm_vars8( qhyd_mass(:,:,:), 5 )
431  call comm_vars8( qcrg(:,:,:), 6 )
432  call comm_vars8( epot(:,:,:), 7 )
433  call comm_wait ( efield(:,:,:,i_lt_x), 1 )
434  call comm_wait ( efield(:,:,:,i_lt_y), 2 )
435  call comm_wait ( efield(:,:,:,i_lt_z), 3 )
436  call comm_wait ( efield(:,:,:,i_lt_abs), 4 )
437  call comm_wait ( qhyd_mass(:,:,:), 5 )
438  call comm_wait ( qcrg(:,:,:), 6 )
439  call comm_wait ( epot(:,:,:), 7 )
440 
441  !--- Calculate lightning path and charge neutralization
442  if( nutr_type == 'MG2001' ) then
443  call atmos_phy_lt_neutralization_mg2001( &
444  ka, ks, ke, & ! [IN]
445  ia, is, ie, & ! [IN]
446  ja, js, je, & ! [IN]
447  kijmax, imax, jmax, & ! [IN]
448  efield(:,:,:,:), & ! [IN]
449  epot(:,:,:), & ! [IN]
450  dens(:,:,:), & ! [IN]
451  qcrg(:,:,:), & ! [IN]
452  qhyd_mass(:,:,:), & ! [IN]
453  num_end(:,:,:,:), & ! [INOUT]
454  lt_path(:,:,:), & ! [INOUT]
455  fls_int_p(:,:,:), & ! [OUT]
456  d_qcrg(:,:,:) ) ! [OUT]
457  elseif( nutr_type == 'F2013' ) then
458  call atmos_phy_lt_neutralization_f2013( &
459  ka, ks, ke, & ! [IN]
460  ia, is, ie, & ! [IN]
461  ja, js, je, & ! [IN]
462  kijmax, & ! [IN]
463  efield(:,:,:,:), & ! [IN]
464  epot(:,:,:), & ! [IN]
465  dens(:,:,:), & ! [IN]
466  qcrg(:,:,:), & ! [IN]
467  qhyd_mass(:,:,:), & ! [IN]
468  num_end(:,:,:,:), & ! [INOUT]
469  lt_path(:,:,:), & ! [INOUT]
470  fls_int_p(:,:,:), & ! [OUT]
471  d_qcrg(:,:,:) ) ! [OUT]
472  endif
473 
474  call comm_vars8( lt_path(:,:,:),1 )
475  call comm_vars8( d_qcrg(:,:,:),2 )
476  call comm_wait ( lt_path(:,:,:),1 )
477  call comm_wait ( d_qcrg(:,:,:),2 )
478 
479  dqneut(:,:,:,:) = 0.0_rp
480  dqneut_real(:,:,:,:) = 0.0_rp
481  !-- Calculate neutralization of each hydrometeor or each category
482  select case( nutr_qhyd )
483  case ( 'TOTAL' )
484 
485  !$omp parallel do &
486  !$omp private(Total_Sarea,r_totalSarea,zerosw)
487  do j = js, je
488  do i = is, ie
489  do k = ks, ke
490  if( abs( d_qcrg(k,i,j) ) > 0.0_rp ) then
491  total_sarea(1) = 0.0_rp
492  do n = 1, qa_lt
493  total_sarea(1) = total_sarea(1) + sarea(k,i,j,n)
494  enddo
495  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea(1)-small )
496  r_totalsarea(1) = 1.0_rp / ( total_sarea(1) + zerosw ) * ( 1.0_rp - zerosw )
497  do n = 1, qa_lt
498  dqneut(k,i,j,n) = d_qcrg(k,i,j)*1.0e+6_rp &
499  * sarea(k,i,j,n) * r_totalsarea(1) / dens(k,i,j)
500  qtrc(k,i,j,n) = qtrc(k,i,j,n) + dqneut(k,i,j,n)
501  dqneut_real(k,i,j,n) = dqneut(k,i,j,n)
502  enddo
503  endif
504  enddo
505  enddo
506  enddo
507 
508  case ( 'Each_POLARITY', 'Each_POLARITY2' )
509 
510  !$omp parallel do &
511  !$omp private(Total_Sarea,r_totalSarea,flg_chrged,pos_crg,neg_crg,frac, &
512  !$omp int_sw,lack, &
513  !$omp positive,negative,zerosw,sw)
514  do j = js, je
515  do i = is, ie
516  do k = ks, ke
517  lack(:) = 0.0_rp
518  if( abs( d_qcrg(k,i,j) ) > 0.0_rp ) then
519 
520  !--- flg whether the charged or not (0.0-> not charged, 1.0->charged)
521  do n = 1, qa_lt
522  flg_chrged(n) = abs(qtrc(k,i,j,n)) >= small
523  enddo
524 
525  total_sarea(:) = 0.0_rp
526  pos_crg = 0.0_rp
527  neg_crg = 0.0_rp
528  do n = 1, qa_lt
529  if ( flg_chrged(n) ) then
530  positive = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
531  negative = 1.0_rp - positive
532  !--- total of positive charge
533  pos_crg = pos_crg + qtrc(k,i,j,n) * positive
534  !--- total of negative charge
535  neg_crg = neg_crg + qtrc(k,i,j,n) * negative
536  !--- Sarea of positively charged hydrometeor
537  total_sarea(1) = total_sarea(1) + sarea(k,i,j,n) * positive
538  !--- Sarea of negatively charged hydrometeor
539  total_sarea(2) = total_sarea(2) + sarea(k,i,j,n) * negative
540  end if
541  end do
542 
543  zerosw = 0.5_rp - sign( 0.5_rp, abs( qcrg(k,i,j) ) - small )
544  frac = d_qcrg(k,i,j) / ( qcrg(k,i,j) + zerosw ) * ( 1.0_rp - zerosw )
545  pos_crg = frac * pos_crg
546  neg_crg = frac * neg_crg
547 
548  !--- remove 0 surface area ( no crg for each porality )
549  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea(1) - small )
550  r_totalsarea(1) = 1.0_rp / ( total_sarea(1) + zerosw ) * ( 1.0_rp - zerosw )
551  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea(2) - small )
552  r_totalsarea(2) = 1.0_rp / ( total_sarea(2) + zerosw ) * ( 1.0_rp - zerosw )
553 
554  diff_qcrg(:) = 0.0_rp
555  crg_rate(:) = 0.0_rp
556  sum_crg(:) = 0.0_rp
557  do n = 1, qa_lt
558  if ( flg_chrged(n) ) then
559  sw = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
560  dqneut(k,i,j,n) = &
561  + pos_crg * sarea(k,i,j,n) * r_totalsarea(1) &
562  * sw & ! sw = 1 positive
563  + neg_crg * sarea(k,i,j,n) * r_totalsarea(2) &
564  * ( 1.0_rp - sw ) ! sw = 0 negative
565  qcrg_before(n) = qtrc(k,i,j,n)
566 
567  if( sw == 1.0_rp ) then
568  qtrc(k,i,j,n) = max( qtrc(k,i,j,n) + dqneut(k,i,j,n), 0.0_rp ) !--- limiter
569  elseif( sw == 0.0_rp ) then
570  qtrc(k,i,j,n) = min( qtrc(k,i,j,n) + dqneut(k,i,j,n), 0.0_rp ) !--- limiter
571  endif
572 
573  int_sw = int( sw ) ! 0-> negative, 1-> positive
574  diff_qcrg(int_sw) = diff_qcrg(int_sw) &
575  + ( qtrc(k,i,j,n) - qcrg_before(n) )
576  sum_crg(int_sw) = sum_crg(int_sw) + qtrc(k,i,j,n)
577  end if
578  enddo
579 
580  if( nutr_qhyd == 'Each_POLARITY2' ) then ! Adjust
581  lack(0) = neg_crg - diff_qcrg(0) ! negative (Should be Positive)
582  lack(1) = pos_crg - diff_qcrg(1) ! positive (Should be Negative)
583 #ifdef DEBUG
584  if( lack(0) > 0.0_rp .and. abs(lack(0)/diff_qcrg(0)) > 1.0e-10_rp ) then
585  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,4E15.7)') &
586  "Large negative for lack(0) ", lack(0)/diff_qcrg(0), neg_crg, lack(0), diff_qcrg(0)
587  endif
588  if( lack(1) < 0.0_rp .and. abs(lack(1)/diff_qcrg(1)) > 1.0e-10_rp ) then
589  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,4E15.7)') &
590  "Large positive for lack(1) ", lack(1)/diff_qcrg(1), pos_crg, lack(1), diff_qcrg(1)
591  endif
592 #endif
593  lack(0) = max( lack(0), 0.0_rp ) ! negative (Should be Positive)
594  lack(1) = min( lack(1), 0.0_rp ) ! positive (Should be Negative)
595  do n = 1, qa_lt
596  if ( flg_chrged(n) ) then
597  sw = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
598  int_sw = int( sw ) ! 0-> negative, 1-> positive
599  if( sum_crg(int_sw) /= 0.0_rp ) then
600  crg_rate(n) = qtrc(k,i,j,n)/sum_crg(int_sw)
601  else
602  crg_rate(n) = 0.0_rp
603  endif
604  qtrc(k,i,j,n) = qtrc(k,i,j,n) + crg_rate(n) * lack(int_sw)
605  endif
606  enddo
607  endif
608 
609  do n = 1, qa_lt
610  if ( flg_chrged(n) ) then
611  dqneut_real(k,i,j,n) = qtrc(k,i,j,n) - qcrg_before(n)
612  endif
613  enddo
614 
615  endif
616 
617  enddo
618  enddo
619  enddo
620 
621  end select
622 
623 
624  do j = js, je
625  do i = is, ie
626 
627  ! calc total charge density
628  do k = ks, ke
629  qcrg(k,i,j) = 0.0_rp
630  dqneut_real_tot(k,i,j) = 0.0_rp
631  do n = 1, qa_lt
632  qcrg(k,i,j) = qcrg(k,i,j) + qtrc(k,i,j,n)
633  dqneut_real_tot(k,i,j) = dqneut_real_tot(k,i,j) + dqneut_real(k,i,j,n)
634  end do
635  qcrg(k,i,j) = qcrg(k,i,j) * dens(k,i,j) * 1.e-6_rp ![fC/kg] -> [nc/m3]
636  enddo
637 
638  enddo
639  enddo
640 
641  !--- Calculate E field
642  call atmos_phy_lt_electric_field( ka, ks, ke, & ! [IN]
643  ia, is, ie, & ! [IN]
644  ja, js, je, & ! [IN]
645  qcrg(:,:,:), & ! [IN]
646  dens(:,:,:), & ! [IN]
647  rhot(:,:,:), & ! [IN]
648  epot(:,:,:), & ! [INOUT]
649  efield(:,:,:,i_lt_x:i_lt_z) ) ! [OUT]
650 
651  !--- Add Total number of path
652  !$omp parallel do
653  do j = js, je
654  do i = is, ie
655  do k = ks, ke
656  efield(k,i,j,i_lt_abs) = sqrt( efield(k,i,j,i_lt_x)*efield(k,i,j,i_lt_x) &
657  + efield(k,i,j,i_lt_y)*efield(k,i,j,i_lt_y) &
658  + efield(k,i,j,i_lt_z)*efield(k,i,j,i_lt_z) )
659  end do
660  end do
661  end do
662 
663  !--- Add Total number of charge neutralization and flash point
664  if ( hist_id(i_qneut) > 0 ) then
665  !$omp parallel do
666  do j = js, je
667  do i = is, ie
668  do k = ks, ke
669  d_qcrg_tot(k,i,j) = d_qcrg_tot(k,i,j) + dqneut_real_tot(k,i,j)*1.e-6_rp ![fC/m3]->[nC/m3]
670  end do
671  end do
672  end do
673  end if
674  if ( hist_id(i_flashpoint) > 0 ) then
675  !$omp parallel do
676  do j = js, je
677  do i = is, ie
678  do k = ks, ke
679  fls_int_p_tot(k,i,j) = fls_int_p_tot(k,i,j) + fls_int_p(k,i,j)
680  end do
681  end do
682  end do
683  end if
684 
685  if ( hist_id(i_posflash) > 0 ) then
686  !$omp parallel do
687  do j = js, je
688  do i = is, ie
689  do k = ks, ke
690  lt_path_tot(k,i,j,1) = lt_path_tot(k,i,j,1) &
691  + 0.5_rp + sign( 0.5_rp,-dqneut_real_tot(k,i,j)-small )
692  end do
693  end do
694  end do
695  end if
696  if ( hist_id(i_negflash) > 0 ) then
697  !$omp parallel do
698  do j = js, je
699  do i = is, ie
700  do k = ks, ke
701  lt_path_tot(k,i,j,2) = lt_path_tot(k,i,j,2) &
702  + 0.5_rp + sign( 0.5_rp, dqneut_real_tot(k,i,j)-small )
703  end do
704  end do
705  end do
706  end if
707  if ( hist_id(i_ltpath) > 0 ) then
708  !$omp parallel do
709  do j = js, je
710  do i = is, ie
711  do k = ks, ke
712  lt_path_tot(k,i,j,3) = lt_path_tot(k,i,j,3) + lt_path(k,i,j)
713  end do
714  end do
715  end do
716  end if
717 
718  call atmos_phy_lt_judge_abse( ka, ks, ke, & ! [IN]
719  ia, is, ie, & ! [IN]
720  ja, js, je, & ! [IN]
721  dens(:,:,:), & ! [IN]
722  efield(:,:,:,i_lt_abs), & ! [IN]
723  emax, & ! [OUT]
724  flg_lt_neut ) ! [OUT]
725 
726  count_neut = count_neut + 1
727 #ifdef DEBUG
728  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,F15.7,A,F15.7,1X,I0)') &
729  'CHECK', &
730  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, count_neut
731 #endif
732 
733  if( flg_lt_neut == 1 .and. emax == emax_old ) then
734  flg_lt_neut = 0
735  if( prc_ismaster ) then
736  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,2E15.7,1X,I0)') &
737  'Eabs value after neutralization is same as previous value, Finish', &
738  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, count_neut
739  endif
740  elseif( flg_lt_neut == 1 .and. count_neut == max_nutr ) then
741  flg_lt_neut = 0
742  if( prc_ismaster ) then
743  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
744  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
745  ' [kV/m], reach maximum neutralization count, Finish', count_neut
746  endif
747  elseif( flg_lt_neut == 1 .and. emax > emax_old ) then
748  flg_lt_neut = 0
749  if( prc_ismaster ) then
750  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
751  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
752  ' [kV/m] larger than previous one by neutralization, back to previous step, Finish', &
753  count_neut
754  endif
755  !---- Back to Charge density as previous step
756  do j = js, je
757  do i = is, ie
758  do k = ks, ke
759  do n = 1, qa_lt
760  qtrc(k,i,j,n) = qtrc(k,i,j,n) - dqneut_real(k,i,j,n)
761  enddo
762  d_qcrg_tot(k,i,j) = d_qcrg_tot(k,i,j) - dqneut_real_tot(k,i,j)*1.e-6_rp
763  d_qcrg(k,i,j) = 0.0_rp
764  enddo
765  enddo
766  enddo
767  elseif( flg_lt_neut == 2 ) then
768  flg_lt_neut = 0
769  if( prc_ismaster ) then
770  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
771  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
772  '[kV/m] After neutralization, Finish' , count_neut
773  endif
774  elseif( flg_lt_neut == 0 ) then
775  if( prc_ismaster ) then
776  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(F15.7,A),1X,I0)') &
777  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
778  '[kV/m] After neutralization, Finish', count_neut
779  endif
780  endif
781 
782  enddo
783 
784  else
785 
786  !$omp parallel do
787  do j = js, je
788  do i = is, ie
789  do k = ks, ke
790  d_qcrg(k,i,j) = 0.0_rp
791  end do
792  end do
793  end do
794 
795  endif
796 
797  !--- For history output
798  do ip = 1, w_nmax
799  call file_history_query( hist_id(ip), hist_sw(ip) )
800  end do
801 
802  if ( hist_sw(i_ex ) ) then
803  !$omp parallel do
804  do j = js, je
805  do i = is, ie
806  do k = ks, ke
807  w3d(k,i,j) = efield(k,i,j,i_lt_x )*1.0e-3_rp ![kV/m]
808  enddo
809  enddo
810  enddo
811  call file_history_put( hist_id(i_ex ), w3d(:,:,:) )
812  endif
813  if ( hist_sw(i_ey ) ) then
814  !$omp parallel do
815  do j = js, je
816  do i = is, ie
817  do k = ks, ke
818  w3d(k,i,j) = efield(k,i,j,i_lt_y )*1.0e-3_rp ![kV/m]
819  enddo
820  enddo
821  enddo
822  call file_history_put( hist_id(i_ey ), w3d(:,:,:) )
823  endif
824  if ( hist_sw(i_ez ) ) then
825  !$omp parallel do
826  do j = js, je
827  do i = is, ie
828  do k = ks, ke
829  w3d(k,i,j) = efield(k,i,j,i_lt_z )*1.0e-3_rp ![kV/m]
830  enddo
831  enddo
832  enddo
833  call file_history_put( hist_id(i_ez ), w3d(:,:,:) )
834  endif
835  if ( hist_sw(i_eabs) ) then
836  !$omp parallel do
837  do j = js, je
838  do i = is, ie
839  do k = ks, ke
840  w3d(k,i,j) = efield(k,i,j,i_lt_abs)*1.0e-3_rp ![kV/m]
841  enddo
842  enddo
843  enddo
844  call file_history_put( hist_id(i_eabs), w3d(:,:,:) )
845  endif
846  if ( hist_sw(i_epot) ) then
847  call file_history_put( hist_id(i_epot), epot(:,:,:) )
848  end if
849  if ( hist_sw(i_qneut) ) then
850  call file_history_put( hist_id(i_qneut), d_qcrg_tot(:,:,:) )
851  !$omp parallel do
852  do j = js, je
853  do i = is, ie
854  do k = ks, ke
855  d_qcrg_tot(k,i,j) = 0.0_rp
856  end do
857  end do
858  end do
859  endif
860  if ( hist_sw(i_ltpath) ) then
861  call file_history_put( hist_id(i_ltpath), lt_path_tot(:,:,:,3) )
862  !$omp parallel do
863  do j = js, je
864  do i = is, ie
865  do k = ks, ke
866  lt_path_tot(k,i,j,3) = 0.0_rp
867  end do
868  end do
869  end do
870  endif
871  if ( hist_sw(i_posflash) ) then
872  call file_history_put( hist_id(i_posflash), lt_path_tot(:,:,:,1) )
873  !$omp parallel do
874  do j = js, je
875  do i = is, ie
876  do k = ks, ke
877  lt_path_tot(k,i,j,1) = 0.0_rp
878  end do
879  end do
880  end do
881  endif
882  if ( hist_sw(i_negflash) ) then
883  call file_history_put( hist_id(i_negflash), lt_path_tot(:,:,:,2) )
884  !$omp parallel do
885  do j = js, je
886  do i = is, ie
887  do k = ks, ke
888  lt_path_tot(k,i,j,2) = 0.0_rp
889  end do
890  end do
891  end do
892  endif
893  if ( hist_sw(i_flashpoint) ) then
894  call file_history_put( hist_id(i_flashpoint), fls_int_p_tot(:,:,:) )
895  !$omp parallel do
896  do j = js, je
897  do i = is, ie
898  do k = ks, ke
899  fls_int_p_tot(k,i,j) = 0.0_rp
900  end do
901  end do
902  end do
903  end if
904 
905 
906  return
907  end subroutine atmos_phy_lt_sato2019_adjustment
908  !-----------------------------------------------------------------------------
911  !-----------------------------------------------------------------------------
912  subroutine atmos_phy_lt_electric_field( &
913  KA, KS, KE, & ! [IN]
914  IA, IS, IE, & ! [IN]
915  JA, JS, JE, & ! [IN]
916  QCRG, & ! [IN]
917  DENS, & ! [IN]
918  RHOT, & ! [IN]
919  E_pot, & ! [OUT]
920  Efield ) ! [OUT]
921  use scale_prc, only: &
922  prc_abort, &
923  prc_ismaster, &
924  prc_myrank
925  use scale_const, only: &
926  eps => const_eps, &
927  epsvac => const_epsvac, &
928  epsair => const_epsair
929  use scale_atmos_grid_cartesc, only: &
930  rcdz => atmos_grid_cartesc_rcdz, &
931  rcdx => atmos_grid_cartesc_rcdx, &
932  rcdy => atmos_grid_cartesc_rcdy, &
933  rfdz => atmos_grid_cartesc_rfdz, &
934  rfdx => atmos_grid_cartesc_rfdx, &
935  rfdy => atmos_grid_cartesc_rfdy
936  use scale_atmos_grid_cartesc_metric, only: &
937  mapf => atmos_grid_cartesc_metric_mapf, &
938  gsqrt => atmos_grid_cartesc_metric_gsqrt, &
939  j13g => atmos_grid_cartesc_metric_j13g, &
940  j23g => atmos_grid_cartesc_metric_j23g, &
941  j33g => atmos_grid_cartesc_metric_j33g
942  use scale_atmos_grid_cartesc_index, only: &
943  i_xyz, &
944  i_xyw, &
945  i_uyw, &
946  i_xvw, &
947  i_uyz, &
948  i_xvz, &
949  i_uvz, &
950  i_xy, &
951  i_uy, &
952  i_xv, &
953  i_uv
954  use scale_prc_cartesc, only: &
955  prc_has_w, &
956  prc_has_e, &
957  prc_has_n, &
958  prc_has_s
959  use scale_comm_cartesc, only: &
960  comm_datatype, &
961  comm_world, &
962  comm_vars8, &
963  comm_wait
964  implicit none
965 
966  integer, intent(in) :: KA, KS, KE
967  integer, intent(in) :: IA, IS, IE
968  integer, intent(in) :: JA, JS, JE
969  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
970  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
971  real(RP), intent(in) :: RHOT (KA,IA,JA) !-- density weighted potential temperature [K kg/m3]
972  real(RP), intent(inout) :: E_pot (KA,IA,JA) !-- Electric potential [V]
973  real(RP), intent(out) :: Efield(KA,IA,JA,3) !-- Electric field [V/m]
974 
975  real(RP) :: eps_air(KA,IA,JA)
976  !--- A x E_pott = - QCRG/epsiron
977  real(RP) :: A(KA,15,IA,JA) !--- A : Laplasian
978  real(RP) :: B(KA,IA,JA) !--- B : -QCRG*DENS/epsiron
979  real(RP) :: E_pot_N(KA,IA,JA) !--- electrical potential calculated by Bi-CGSTAB
980 
981  integer :: i, j, k, ijk, ierror
982  real(RP) :: iprod, buf
983 
984  call prof_rapstart('LT_E_field', 1)
985 
986  iprod = 0.0_rp
987  !$omp parallel do reduction(+:iprod)
988  do j = js, je
989  do i = is, ie
990  do k = ks, ke
991  iprod = iprod + abs( qcrg(k,i,j) )
992  enddo
993  enddo
994  enddo
995 
996  call mpi_allreduce(iprod, buf, 1, comm_datatype, mpi_sum, comm_world, ierror)
997 
998  if( buf <= eps ) then
999  !$omp parallel do
1000  do j = js, je
1001  do i = is, ie
1002  do k = ks, ke
1003  e_pot(k,i,j) = 0.0_rp
1004  efield(k,i,j,1:3) = 0.0_rp
1005  end do
1006  end do
1007  end do
1008 
1009  return
1010  endif
1011 
1012  !$omp parallel do
1013  do j = js, je
1014  do i = is, ie
1015  do k = ks, ke
1016  eps_air(k,i,j) = epsvac * epsair !--- temporary, dependency of epsiron on P and T will be implemented
1017  b(k,i,j) = - qcrg(k,i,j)/eps_air(k,i,j) * 1.0e-9_rp ! [nC/m3] -> [C/m3 * m/F] = [V/m2]
1018  enddo
1019  enddo
1020  enddo
1021 
1022  !---- fill halo
1023  call comm_vars8( eps_air,1 )
1024  call comm_vars8( b, 2 )
1025 
1026  !---- input vector A
1027  !$omp parallel do
1028  do j = js, je
1029  do i = is, ie
1030  do k = ks, ke
1031  ! (k,i,j)
1032  a(k,1,i,j) = &
1033  - mapf(i,j ,1,i_xy)*mapf(i,j ,1,i_xy)*rcdx(i )*rfdx(i) &
1034  - mapf(i,j ,1,i_xy)*mapf(i,j ,1,i_xy)*rcdx(i-1)*rfdx(i) &
1035  + mapf(i,j ,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k ,i,j,2)*0.5_rp*rfdx(i)*rfdz(k) &
1036  - mapf(i,j ,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k-1,i,j,1)*0.5_rp*rfdx(i)*rfdz(k) &
1037  - mapf(i,j ,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k ,i,j,2)*0.5_rp*rfdx(i)*rfdz(k) &
1038  + mapf(i,j ,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k-1,i,j,1)*0.5_rp*rfdx(i)*rfdz(k) &
1039  - j13g(k,i,j,i_xyz)*j13g(k ,i,j,i_xyw)*rcdz(k)*rfdz(k) &
1040  - j13g(k,i,j,i_xyz)*j13g(k-1,i,j,i_xyw)*rcdz(k)*rfdz(k) &
1041  - mapf(i,j ,2,i_xy)*mapf(i,j ,2,i_xy)*rcdy(j )*rfdy(j) &
1042  - mapf(i,j ,2,i_xy)*mapf(i,j ,2,i_xy)*rcdy(j-1)*rfdy(j) &
1043  + mapf(i,j ,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k ,i,j,2)*0.5_rp*rfdy(j)*rfdz(k) &
1044  - mapf(i,j ,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k-1,i,j,1)*0.5_rp*rfdy(j)*rfdz(k) &
1045  - mapf(i,j ,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k ,i,j,2)*0.5_rp*rfdy(j)*rfdz(k) &
1046  + mapf(i,j ,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k-1,i,j,1)*0.5_rp*rfdy(j)*rfdz(k) &
1047  - j23g(k,i,j,i_xyz)*j23g(k ,i,j,i_xyw)*rcdz(k)*rfdz(k) &
1048  - j23g(k,i,j,i_xyz)*j23g(k-1,i,j,i_xyw)*rcdz(k)*rfdz(k) &
1049  - j33g*j33g*rfdz(k)*rcdz(k) &
1050  - j33g*j33g*rfdz(k)*rcdz(k-1)
1051 
1052  ! (k-1,i,j)
1053  a(k,2,i,j) = &
1054  - mapf(i,j,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k-1,i,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1055  + mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k-1,i,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1056  + j13g(k,i,j,i_xyz)*j13g(k-1,i,j,i_xyw)*rfdz(k)*rcdz(k-1) &
1057  - mapf(i,j,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k-1,i,j,2)*0.50_rp*rfdy(j)*rfdz(k) &
1058  + mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k-1,i,j,2)*0.50_rp*rfdy(j)*rfdz(k) &
1059  + j23g(k,i,j,i_xyz)*j23g(k-1,i,j,i_xyw)*rfdz(k)*rcdz(k-1) &
1060  + j33g*j33g*rfdz(k)*rcdz(k-1)
1061 
1062  ! (k+1,i,j)
1063  a(k,3,i,j) = &
1064  mapf(i,j,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k,i,j,1)*0.50_rp*rfdx(i)*rfdz(k) &
1065  - mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k,i,j,1)*0.50_rp*rfdx(i)*rfdz(k) &
1066  + j13g(k,i,j,i_xyz)*j13g(k,i,j,i_xyw)*rfdz(k)*rcdz(k) &
1067  + mapf(i,j,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k,i,j,1)*0.50_rp*rfdy(j)*rfdz(k) &
1068  - mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k,i,j,1)*0.50_rp*rfdy(j)*rfdz(k) &
1069  + j23g(k,i,j,i_xyz)*j23g(k,i,j,i_xyw)*rfdz(k)*rcdz(k) &
1070  + j33g*j33g*rfdz(k)*rcdz(k)
1071 
1072  ! (k,i-1,j)
1073  a(k,4,i,j) = &
1074  mapf(i,j,1,i_xy)*mapf(i-1,j,1,i_xy)*rfdx(i)*rcdx(i-1) &
1075  - mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k ,i-1,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1076  + mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k-1,i-1,j,1)*0.50_rp*rfdx(i)*rfdz(k) &
1077  - j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k ,i-1,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1078  + j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k-1,i-1,j,1)*0.50_rp*rfdx(i)*rfdz(k)
1079 
1080  ! (k,i+1,j)
1081  a(k,5,i,j) = &
1082  mapf(i,j,1,i_xy)*mapf(i+1,j,1,i_xy)*rfdx(i)*rcdx(i) &
1083  + mapf(i,j,1,i_xy)*j13g(k,i,j,i_uyz)*f2h(k ,i+1,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1084  - mapf(i,j,1,i_xy)*j13g(k,i,j,i_uyz)*f2h(k-1,i+1,j,1)*0.50_rp*rfdx(i)*rfdz(k) &
1085  + j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k ,i+1,j,2)*0.50_rp*rfdx(i)*rfdz(k) &
1086  - j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k-1,i+1,j,1)*0.50_rp*rfdx(i)*rfdz(k)
1087 
1088  ! (k,i,j-1)
1089  a(k,6,i,j) = &
1090  mapf(i,j,2,i_xy)*mapf(i,j-1,2,i_xy)*rfdy(j)*rcdy(j-1) &
1091  - mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k ,i,j-1,2)*0.50_rp*rfdy(j)*rfdz(k) &
1092  + mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k-1,i,j-1,1)*0.50_rp*rfdy(j)*rfdz(k) &
1093  - j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k ,i,j-1,2)*0.50_rp*rfdy(j)*rfdz(k) &
1094  + j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k-1,i,j-1,1)*0.50_rp*rfdy(j)*rfdz(k)
1095 
1096  ! (k,i,j+1)
1097  a(k,7,i,j) = &
1098  mapf(i,j,2,i_xy)*mapf(i,j+1,2,i_xy)*rfdy(j)*rcdy(j) &
1099  + mapf(i,j,2,i_xy)*j23g(k,i,j,i_xvz)*f2h(k ,i,j+1,2)*0.50_rp*rfdy(j)*rfdz(k) &
1100  - mapf(i,j,2,i_xy)*j23g(k,i,j,i_xvz)*f2h(k-1,i,j+1,1)*0.50_rp*rfdy(j)*rfdz(k) &
1101  + j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k ,i,j+1,2)*0.50_rp*rfdy(j)*rfdz(k) &
1102  - j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k-1,i,j+1,1)*0.50_rp*rfdy(j)*rfdz(k)
1103 
1104  ! (k-1,i-1,j)
1105  a(k,8,i,j) = &
1106  mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k-1,i-1,j,2)*0.5_rp*rfdx(i)*rfdz(k) &
1107  + j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k-1,i-1,j,2)*0.5_rp*rfdx(i)*rfdz(k)
1108 
1109  ! (k-1,i+1,j)
1110  a(k,9,i,j) = &
1111  - mapf(i,j,1,i_xy)*j13g(k,i,j,i_uyz)*f2h(k-1,i+1,j,2)*0.5_rp*rfdx(i)*rfdz(k) &
1112  - j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k-1,i+1,j,2)*0.5_rp*rfdx(i)*rfdz(k)
1113 
1114  ! (k-1,i,j-1)
1115  a(k,10,i,j) = &
1116  mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k-1,i,j-1,2)*0.5_rp*rfdy(j)*rfdz(k) &
1117  + j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k-1,i,j-1,2)*0.5_rp*rfdy(j)*rfdz(k)
1118 
1119  ! (k-1,i,j+1)
1120  a(k,11,i,j) = &
1121  - mapf(i,j,2,i_xy)*j23g(k,i,j,i_xvz)*f2h(k-1,i,j+1,2)*0.5_rp*rfdy(j)*rfdz(k) &
1122  - j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k-1,i,j+1,2)*0.5_rp*rfdy(j)*rfdz(k)
1123 
1124  ! (k+1,i-1,j)
1125  a(k,12,i,j) = &
1126  - mapf(i,j,1,i_xy)*j13g(k,i-1,j,i_uyz)*f2h(k,i-1,j,1)*0.5_rp*rfdx(i)*rfdz(k) &
1127  - j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k,i-1,j,1)*0.5_rp*rfdx(i)*rfdz(k)
1128 
1129  ! (k+1,i+1,j)
1130  a(k,13,i,j) = &
1131  mapf(i,j,1,i_xy)*j13g(k,i,j,i_uyz)*f2h(k,i+1,j,1)*0.5_rp*rfdx(i)*rfdz(k) &
1132  + j13g(k,i,j,i_xyz)*mapf(i,j,1,i_xy)*f2h(k,i+1,j,1)*0.5_rp*rfdx(i)*rfdz(k)
1133 
1134  ! (k+1,i,j-1)
1135  a(k,14,i,j) = &
1136  - mapf(i,j,2,i_xy)*j23g(k,i,j-1,i_xvz)*f2h(k,i,j-1,1)*0.5_rp*rfdy(j)*rfdz(k) &
1137  - j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k,i,j-1,1)*0.5_rp*rfdy(j)*rfdz(k)
1138 
1139  ! (k+1,i,j+1)
1140  a(k,15,i,j) = &
1141  mapf(i,j,2,i_xy)*j23g(k,i,j,i_xvz)*f2h(k,i,j+1,1)*0.5_rp*rfdy(j)*rfdz(k) &
1142  + j23g(k,i,j,i_xyz)*mapf(i,j,2,i_xy)*f2h(k,i,j+1,1)*0.5_rp*rfdy(j)*rfdz(k)
1143 
1144  enddo
1145  enddo
1146  enddo
1147 
1148 
1149  !$omp parallel do
1150  do j = js, je
1151  do i = is, ie
1152  do k = ks, ke
1153  e_pot_n(k,i,j) = e_pot(k,i,j) !-- initial value -> previous step value
1154  end do
1155  end do
1156  end do
1157  call comm_vars8( e_pot_n, 3 )
1158 
1159 
1160  call comm_wait ( eps_air, 1 )
1161  call comm_wait ( b, 2 )
1162  call comm_wait ( e_pot_n, 3 )
1163 
1164  !--- calcuclate counter matrix
1165  call atmos_phy_lt_solve_bicgstab( &
1166  ka, ks, ke, & ! (in)
1167  ia, is, ie, & ! (in)
1168  ja, js, je, & ! (in)
1169  e_pot, & ! (out)
1170  e_pot_n, & ! (in)
1171  a, b ) ! (in)
1172 
1173  call comm_vars8( e_pot, 1 )
1174  call comm_wait ( e_pot, 1, .true. )
1175 
1176  !$omp parallel do
1177  do j = 1, ja
1178  do i = 1, ia
1179  e_pot(1:ks-1,i,j) = 0.0_rp
1180  e_pot(ke+1:ka,i,j) = 0.0_rp
1181  enddo
1182  enddo
1183 
1184  !---- Calculate Electrical Field
1185  !$omp parallel do
1186  do j = js, je
1187  do i = is, ie
1188  do k = ks, ke
1189  efield(k,i,j,i_lt_x) = - mapf(i,j,1,i_xyz)/gsqrt(k,i,j,i_xyz) * &
1190  ( &
1191  ( gsqrt(k,i+1,j,i_xyz)*e_pot(k,i+1,j) - gsqrt(k,i-1,j,i_xyz)*e_pot(k,i-1,j) ) &
1192  * rcdx(i) * 0.5_rp &
1193  + ( j13g(k+1,i,j,i_uyw)*gsqrt(k+1,i,j,i_uyw)*e_pot(k+1,i,j) &
1194  - j13g(k-1,i,j,i_uyw)*gsqrt(k-1,i,j,i_uyw)*e_pot(k-1,i,j) &
1195  ) &
1196  * rfdz(k) * 0.5_rp &
1197  )
1198  efield(k,i,j,i_lt_y) = - mapf(i,j,2,i_xyz)/gsqrt(k,i,j,i_xyz) * &
1199  ( &
1200  ( gsqrt(k,i,j+1,i_xyz)*e_pot(k,i,j+1) - gsqrt(k,i,j-1,i_xyz)*e_pot(k,i,j-1) ) &
1201  * rcdy(j) * 0.5_rp &
1202  + ( j23g(k+1,i,j,i_uyw)*gsqrt(k+1,i,j,i_uyw)*e_pot(k+1,i,j) &
1203  - j23g(k-1,i,j,i_uyw)*gsqrt(k-1,i,j,i_uyw)*e_pot(k-1,i,j) &
1204  ) &
1205  * rfdz(k) * 0.5_rp &
1206  )
1207  efield(k,i,j,i_lt_z) = - 1.0_rp/gsqrt(k,i,j,i_xyz) * &
1208  ( j33g*e_pot(k+1,i ,j )*gsqrt(k+1,i,j,i_xyz) &
1209  - j33g*e_pot(k-1,i ,j )*gsqrt(k-1,i,j,i_xyz) &
1210  ) &
1211  * rfdz(k) * 0.5_rp
1212  enddo
1213  enddo
1214  enddo
1215 
1216  call prof_rapend('LT_E_field', 1)
1217 
1218  return
1219  end subroutine atmos_phy_lt_electric_field
1220 
1221  subroutine atmos_phy_lt_solve_bicgstab( &
1222  KA, KS, KE, & ! [IN]
1223  IA, IS, IE, & ! [IN]
1224  JA, JS, JE, & ! [IN]
1225  PHI_N, & ! [OUT]
1226  PHI, & ! [IN]
1227  M, B ) ! [IN]
1228  use scale_prc, only: &
1229  prc_abort, &
1230  prc_ismaster
1231  use scale_comm_cartesc, only: &
1232  comm_datatype, &
1233  comm_world, &
1234  comm_vars8, &
1235  comm_wait
1236  implicit none
1237  integer, intent(in) :: KA, KS, KE
1238  integer, intent(in) :: IA, IS, IE
1239  integer, intent(in) :: JA, JS, JE
1240  real(RP), intent(out) :: PHI_N(KA,IA,JA)
1241  real(RP), intent(in) :: PHI(KA,IA,JA)
1242  real(RP), intent(in) :: M(KA,15,IA,JA)
1243  real(RP), intent(in) :: B(KA,IA,JA)
1244 
1245  real(RP) :: r0(KA,IA,JA)
1246 
1247  real(RP) :: p(KA,IA,JA)
1248  real(RP) :: Mp(KA,IA,JA)
1249  real(RP) :: s(KA,IA,JA)
1250  real(RP) :: Ms(KA,IA,JA)
1251  real(RP) :: al, be, w
1252 
1253  real(RP), pointer :: r(:,:,:)
1254  real(RP), pointer :: rn(:,:,:)
1255  real(RP), pointer :: swap(:,:,:)
1256  real(RP), target :: v0(KA,IA,JA)
1257  real(RP), target :: v1(KA,IA,JA)
1258  real(RP) :: r0r
1259  real(RP) :: norm, error, error2
1260 
1261  real(RP) :: iprod(2)
1262  real(RP) :: buf(2)
1263 
1264  integer :: k, i, j
1265  integer :: iis, iie, jjs, jje
1266  integer :: iter
1267  integer :: ierror
1268 
1269  r => v0
1270  rn => v1
1271 
1272  call mul_matrix( ka, ks, ke, & ! (in)
1273  ia, is, ie, & ! (in)
1274  ja, js, je, & ! (in)
1275  v1, m, phi ) ! v1 = M x0
1276 
1277  norm = 0.0_rp
1278  !$omp parallel do reduction(+:norm)
1279  do j = js, je
1280  do i = is, ie
1281  do k = ks, ke
1282  norm = norm + b(k,i,j)**2
1283  enddo
1284  enddo
1285  enddo
1286 
1287  ! r = b - M x0
1288  !$omp parallel do
1289  do j = js, je
1290  do i = is, ie
1291  do k = ks, ke
1292  r(k,i,j) = b(k,i,j) - v1(k,i,j)
1293  enddo
1294  enddo
1295  enddo
1296 
1297  !$omp parallel do
1298  do j = js, je
1299  do i = is, ie
1300  do k = ks, ke
1301  r0(k,i,j) = r(k,i,j)
1302  p(k,i,j) = r(k,i,j)
1303  enddo
1304  enddo
1305  enddo
1306 
1307  r0r = 0.0_rp
1308  !$omp parallel do reduction(+:r0r)
1309  do j = js, je
1310  do i = is, ie
1311  do k = ks, ke
1312  r0r = r0r + r0(k,i,j) * r(k,i,j)
1313  enddo
1314  enddo
1315  enddo
1316 
1317  !$omp parallel do
1318  do j = js-1, je+1
1319  do i = is-1, ie+1
1320  do k = ks, ke
1321  phi_n(k,i,j) = phi(k,i,j)
1322  end do
1323  end do
1324  end do
1325 
1326 
1327  iprod(1) = r0r
1328  iprod(2) = norm
1329  call mpi_allreduce(iprod, buf, 2, comm_datatype, mpi_sum, comm_world, ierror)
1330  r0r = buf(1)
1331  norm = buf(2)
1332 
1333  error2 = norm
1334 
1335  do iter = 1, itmax
1336 
1337  error = 0.0_rp
1338  !$omp parallel do reduction(+:error)
1339  do j = js, je
1340  do i = is, ie
1341  do k = ks, ke
1342  error = error + r(k,i,j)**2
1343  enddo
1344  enddo
1345  enddo
1346  call mpi_allreduce(error, buf, 1, comm_datatype, mpi_sum, comm_world, ierror)
1347  error = buf(1)
1348 
1349  if ( sqrt(error/norm) < epsilon ) then
1350  log_info("ATMOS_PHY_LT_Efield",'(a,1x,i0,1x,2e15.7)') "Bi-CGSTAB converged:", iter, sqrt(error/norm),norm
1351  exit
1352  endif
1353  error2 = error
1354 
1355  call comm_vars8( p, 1 )
1356  call comm_wait ( p, 1 )
1357  call mul_matrix( ka, ks, ke, & ! (in)
1358  ia, is, ie, & ! (in)
1359  ja, js, je, & ! (in)
1360  mp, m, p )
1361 
1362  iprod(1) = 0.0_rp
1363  !$omp parallel do reduction(+:iprod)
1364  do j = js, je
1365  do i = is, ie
1366  do k = ks, ke
1367  iprod(1) = iprod(1) + r0(k,i,j) * mp(k,i,j)
1368  enddo
1369  enddo
1370  enddo
1371  call mpi_allreduce(iprod, buf, 1, comm_datatype, mpi_sum, comm_world, ierror)
1372  if ( buf(1) == 0.0_rp ) then
1373  log_info("ATMOS_PHY_LT_Efield",'(a,1x,e15.7,1x,i10)') 'Buf(1) is zero(Bi-CGSTAB) skip:', buf(1), iter
1374  exit
1375  endif
1376  al = r0r / buf(1) ! (r0,r) / (r0,Mp)
1377 
1378  !$omp parallel do
1379  do j = js, je
1380  do i = is, ie
1381  do k = ks, ke
1382  s(k,i,j) = r(k,i,j) - al*mp(k,i,j)
1383  enddo
1384  enddo
1385  enddo
1386 
1387  call comm_vars8( s, 1 )
1388  call comm_wait ( s, 1 )
1389  call mul_matrix( ka, ks, ke, & ! (in)
1390  ia, is, ie, & ! (in)
1391  ja, js, je, & ! (in)
1392  ms, m, s )
1393  iprod(1) = 0.0_rp
1394  iprod(2) = 0.0_rp
1395  !$omp parallel do reduction(+:iprod)
1396  do j = js, je
1397  do i = is, ie
1398  do k = ks, ke
1399  iprod(1) = iprod(1) + ms(k,i,j) * s(k,i,j)
1400  iprod(2) = iprod(2) + ms(k,i,j) * ms(k,i,j)
1401  enddo
1402  enddo
1403  enddo
1404  call mpi_allreduce(iprod, buf, 2, comm_datatype, mpi_sum, comm_world, ierror)
1405  if ( buf(2) == 0.0_rp ) then
1406  log_info("ATMOS_PHY_LT_Efield",'(a,1x,e15.7,1x,i10)') 'Buf(2) is zero(Bi-CGSTAB) skip:', buf(2), iter
1407  exit
1408  endif
1409  w = buf(1) / buf(2) ! (Ms,s) / (Ms,Ms)
1410 
1411  !$omp parallel do
1412  do j = js, je
1413  do i = is, ie
1414  do k = ks, ke
1415  phi_n(k,i,j) = phi_n(k,i,j) + al*p(k,i,j) + w*s(k,i,j)
1416  enddo
1417  enddo
1418  enddo
1419 
1420  !$omp parallel do
1421  do j = js, je
1422  do i = is, ie
1423  do k = ks, ke
1424  rn(k,i,j) = s(k,i,j) - w*ms(k,i,j)
1425  enddo
1426  enddo
1427  enddo
1428 
1429  iprod(1) = 0.0_rp
1430  !$omp parallel do reduction(+:iprod)
1431  do j = js, je
1432  do i = is, ie
1433  do k = ks, ke
1434  iprod(1) = iprod(1) + r0(k,i,j) * rn(k,i,j)
1435  enddo
1436  enddo
1437  enddo
1438 
1439  be = al/w / r0r
1440 
1441  call mpi_allreduce(iprod, r0r, 1, comm_datatype, mpi_sum, comm_world, ierror)
1442 
1443  be = be * r0r ! al/w * (r0,rn)/(r0,r)
1444 
1445  !$omp parallel do
1446  do j = js, je
1447  do i = is, ie
1448  do k = ks, ke
1449  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mp(k,i,j) )
1450  enddo
1451  enddo
1452  enddo
1453 
1454  swap => rn
1455  rn => r
1456  r => swap
1457 
1458  if ( r0r == 0.0_rp ) then
1459  log_info("ATMOS_PHY_LT_Efield",'(a,1x,i0,1x,3e15.7)') "Inner product of r0 and r_itr is zero(Bi-CGSTAB) :", &
1460  iter, r0r, sqrt(error/norm), norm
1461  exit
1462  endif
1463 
1464  enddo
1465 
1466  if ( iter >= itmax ) then
1467  if( prc_ismaster ) then
1468  log_warn("ATMOS_PHY_LT_solve_bicgstab",'(a,1x,2e15.7)') 'Bi-CGSTAB not converged:', error, norm
1469  log_warn_cont('(a,1x,2e15.7)') 'Bi-CGSTAB not converged:', epsilon, sqrt(error/norm)
1470  log_warn_cont('(a,1x,2e15.7)') 'xxx epsilon(set,last)=', epsilon, sqrt(error/norm)
1471  if( error /= error ) then
1472  log_error("ATMOS_PHY_LT_solve_bicgstab",*) 'xxx error or norm is NaN Stop!'
1473  call prc_abort
1474  endif
1475  endif
1476  endif
1477 
1478  return
1479  end subroutine atmos_phy_lt_solve_bicgstab
1480 
1481  subroutine mul_matrix( KA,KS,KE, &
1482  IA,IS,IE, &
1483  JA,JS,JE, &
1484  V, M, C)
1485  implicit none
1486  integer, intent(in) :: KA, KS, KE
1487  integer, intent(in) :: IA, IS, IE
1488  integer, intent(in) :: JA, JS, JE
1489  real(RP), intent(out) :: V(KA,IA,JA)
1490  real(RP), intent(in) :: M(KA,15,IA,JA)
1491  real(RP), intent(in) :: C(KA,IA,JA)
1492 
1493  integer :: k, i, j
1494 
1495  !$omp parallel do
1496  do j = js, je
1497  do i = is, ie
1498  do k = ks+1, ke-1
1499  v(k,i,j) = m(k,1,i,j) * c(k ,i ,j ) &
1500  + m(k,2,i,j) * c(k-1,i ,j ) &
1501  + m(k,3,i,j) * c(k+1,i ,j ) &
1502  + m(k,4,i,j) * c(k ,i-1,j ) &
1503  + m(k,5,i,j) * c(k ,i+1,j ) &
1504  + m(k,6,i,j) * c(k ,i ,j-1) &
1505  + m(k,7,i,j) * c(k ,i ,j+1) &
1506  + m(k,8,i,j) * c(k-1,i-1,j ) &
1507  + m(k,9,i,j) * c(k-1,i+1,j ) &
1508  + m(k,10,i,j)* c(k-1,i ,j-1) &
1509  + m(k,11,i,j)* c(k-1,i ,j+1) &
1510  + m(k,12,i,j)* c(k+1,i-1,j ) &
1511  + m(k,13,i,j)* c(k+1,i+1,j ) &
1512  + m(k,14,i,j)* c(k+1,i ,j-1) &
1513  + m(k,15,i,j)* c(k+1,i ,j+1)
1514  enddo
1515  v(ks,i,j) = m(ks,1,i,j) * c(ks ,i ,j ) &
1516  + m(ks,3,i,j) * c(ks+1,i ,j ) &
1517  + m(ks,4,i,j) * c(ks ,i-1,j ) &
1518  + m(ks,5,i,j) * c(ks ,i+1,j ) &
1519  + m(ks,6,i,j) * c(ks ,i ,j-1) &
1520  + m(ks,7,i,j) * c(ks ,i ,j+1) &
1521  + m(ks,12,i,j)* c(ks+1,i-1,j ) &
1522  + m(ks,13,i,j)* c(ks+1,i+1,j ) &
1523  + m(ks,14,i,j)* c(ks+1,i ,j-1) &
1524  + m(ks,15,i,j)* c(ks+1,i ,j+1)
1525  v(ke,i,j) = m(ke,1,i,j) * c(ke ,i ,j ) &
1526  + m(ke,2,i,j) * c(ke-1,i ,j ) &
1527  + m(ke,4,i,j) * c(ke ,i-1,j ) &
1528  + m(ke,5,i,j) * c(ke ,i+1,j ) &
1529  + m(ke,6,i,j) * c(ke ,i ,j-1) &
1530  + m(ke,7,i,j) * c(ke ,i ,j+1) &
1531  + m(ke,8,i,j) * c(ke-1,i-1,j ) &
1532  + m(ke,9,i,j) * c(ke-1,i+1,j ) &
1533  + m(ke,10,i,j)* c(ke-1,i ,j-1) &
1534  + m(ke,11,i,j)* c(ke-1,i ,j+1)
1535  enddo
1536  enddo
1537 
1538  return
1539  end subroutine mul_matrix
1540  !-----------------------------------------------------------------------------
1541  function f2h( k,i,j,p )
1542 
1543  use scale_atmos_grid_cartesc, only: &
1544  cdz => atmos_grid_cartesc_cdz
1545  use scale_atmos_grid_cartesc_metric, only: &
1546  gsqrt => atmos_grid_cartesc_metric_gsqrt
1547  use scale_atmos_grid_cartesc_index, only: &
1548  i_xyz
1549 
1550  implicit none
1551 
1552  real(RP) :: f2h
1553  integer, intent(in) :: k, i, j, p
1554 
1555  f2h = (cdz(k+p-1)*gsqrt(k+p-1,i,j,i_xyz) &
1556  / (cdz(k)*gsqrt(k,i,j,i_xyz)+cdz(k+1)*gsqrt(k+1,i,j,i_xyz)))
1557 
1558  end function
1559  !-----------------------------------------------------------------------------
1561  !-----------------------------------------------------------------------------
1562  subroutine atmos_phy_lt_neutralization_mg2001( &
1563  KA, KS, KE, & ! [IN]
1564  IA, IS, IE, & ! [IN]
1565  JA, JS, JE, & ! [IN]
1566  KIJMAX, & ! [IN]
1567  IMAX, & ! [IN]
1568  JMAX, & ! [IN]
1569  Efield, & ! [IN]
1570  E_pot, & ! [IN]
1571  DENS, & ! [IN]
1572  QCRG, & ! [IN]
1573  QHYD, & ! [IN]
1574  NUM_end, & ! [INOUT]
1575  LT_path, & ! [INOUT]
1576  fls_int_p, & ! [OUT]
1577  d_QCRG ) ! [OUT]
1578  use scale_const, only: &
1579  eps => const_eps, &
1580  epsvac => const_epsvac, &
1581  epsair => const_epsair, &
1582  large_num => const_huge
1583  use scale_atmos_grid_cartesc, only: &
1584  cz => atmos_grid_cartesc_cz, &
1585  cx => atmos_grid_cartesc_cx, &
1586  cy => atmos_grid_cartesc_cy, &
1587  cdz => atmos_grid_cartesc_cdz, &
1588  cdx => atmos_grid_cartesc_cdx, &
1589  cdy => atmos_grid_cartesc_cdy, &
1590  rfdz => atmos_grid_cartesc_rfdz, &
1591  rfdx => atmos_grid_cartesc_rfdx, &
1592  rfdy => atmos_grid_cartesc_rfdy
1593  use scale_prc_cartesc, only: &
1594  prc_has_w, &
1595  prc_has_e, &
1596  prc_has_n, &
1597  prc_has_s, &
1598  prc_w, &
1599  prc_n, &
1600  prc_e, &
1601  prc_s, &
1602  prc_nw, &
1603  prc_ne, &
1604  prc_sw, &
1605  prc_se, &
1606  prc_next, &
1607  prc_num_x, &
1608  prc_num_y
1609  use scale_prc, only: &
1610  prc_abort, &
1611  prc_ismaster, &
1612  prc_nprocs, &
1613  prc_masterrank, &
1614  prc_mpibarrier, &
1615  prc_myrank
1616  use scale_comm_cartesc, only: &
1617  comm_datatype, &
1618  comm_world, &
1619  comm_vars8, &
1620  comm_wait, &
1621  comm_bcast
1622  use scale_random, only: &
1623  random_uniform
1624  implicit none
1625 
1626  integer, intent(in) :: KA, KS, KE
1627  integer, intent(in) :: IA, IS, IE
1628  integer, intent(in) :: JA, JS, JE
1629  integer, intent(in) :: KIJMAX, IMAX, JMAX
1630  real(RP), intent(in) :: Efield (KA,IA,JA,4) !-- Electric field [V/m] (4->|E|)
1631  real(RP), intent(in) :: E_pot (KA,IA,JA) !-- Electric potential [V]
1632  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
1633  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
1634  real(RP), intent(in) :: QHYD (KA,IA,JA) !-- Total hydrometeor mixing ratio
1635  real(RP), intent(inout) :: NUM_end (KA,IA,JA,3) !-- Number of Lightning stop grid
1636  real(RP), intent(inout) :: LT_path (KA,IA,JA) !-- Number of path
1637  real(RP), intent(out) :: fls_int_p(KA,IA,JA) !-- Flash initiation point (0->no flash, 1->flash start at the point)
1638  real(RP), intent(out) :: d_QCRG (KA,IA,JA) !-- Charge density [nC/m3]
1639 
1640  real(RP) :: E_det, E_x, E_y, E_z, E_max
1641  real(RP) :: dqrho_pls(KA,IA,JA), dqrho_mns(KA,IA,JA)
1642  real(RP) :: L_path(KA,IA,JA)
1643  real(RP) :: sumdqrho_pls, sumdqrho_mns
1644  real(RP) :: Ncrit, dqrho_cor
1645  real(RP) :: randnum(1)
1646  real(RP) :: rprod1, rprod2, rbuf, rbuf2, phi_end(2)
1647  integer :: Eovid(4,KIJMAX)
1648  integer :: Npls, Nmns, Ntot
1649  integer :: count1(PRC_nprocs)
1650  integer :: countindx(PRC_nprocs+1)
1651  integer :: own_prc_total !--- number of grid whose |E| > Eint-dEint for each process
1652  integer :: iprod(2), ibuf(6), status(MPI_STATUS_SIZE)
1653  integer :: init_point, rank_initpoint, grid_initpoint
1654  integer :: current_prc !--- 1->lightning flash path is include, 0-> nopath in the node
1655  integer :: comm_flag, comm_target, stop_flag, corr_flag(2)
1656  integer :: end_grid(4), wild_grid(6)
1657  real(RP) :: pm
1658  integer :: k, i, j, ipp, iq, direction, ierr
1659  integer :: k1, i1, j1, k2, i2, j2, sign_flash
1660  integer :: wild_flag, count_path
1661  integer :: flg_num_end(KA,IA,JA,3)
1662  real(RP) :: Eint_hgt(KA,IA,JA)
1663 
1664  call prof_rapstart('LT_neut_MG2001', 1)
1665 
1666  !$omp parallel do
1667  do j = js, je
1668  do i = is, ie
1669  do k = ks, ke
1670  fls_int_p(k,i,j) = 0.0_rp
1671  end do
1672  end do
1673  end do
1674 
1675  do count_path = 1, nutr_itmax
1676 
1677 
1678  !--- search grid whose Efield is over threshold of flash inititation
1679  own_prc_total = 0
1680  do j = js, je
1681  do i = is, ie
1682  do k = ks, ke
1683  eint_hgt(k,i,j) = min( eint*dens(k,i,j)/rho0,eint )*flg_eint_hgt &
1684  + ( eint-deleint )*( 1.0_rp-flg_eint_hgt )
1685  if( flg_eint_hgt == 1.0_rp .and. eint_hgt(k,i,j) < estop ) then
1686  eint_hgt(k,i,j) = large_num !--- ignore when Eint < Estop
1687  endif
1688  e_det = efield(k,i,j,i_lt_abs) - eint_hgt(k,i,j)
1689  if( e_det > 0.0_rp ) then
1690  own_prc_total = own_prc_total + 1
1691  eovid(1,own_prc_total) = i
1692  eovid(2,own_prc_total) = j
1693  eovid(3,own_prc_total) = k
1694  eovid(4,own_prc_total) = prc_myrank
1695  endif
1696  enddo
1697  enddo
1698  enddo
1699 
1700  call mpi_allgather( own_prc_total, 1, mpi_integer, &
1701  count1, 1, mpi_integer, &
1702  comm_world, ierr )
1703  countindx(1) = 0
1704 
1705  !**** countindx(PROC_nprocs+1) -> total number of grid with |E|>E_threthold
1706  do ipp = 1, prc_nprocs
1707  countindx(ipp+1) = countindx(ipp) + count1(ipp)
1708  enddo
1709 
1710  !---- select initial point of flash by random select
1711  !**** rank_initpoint -> rank number including initpoint
1712  !**** grid_initpoint -> grid number of init point in rank (rank_initpoint)
1713  if( prc_ismaster ) then
1714  call random_uniform(randnum)
1715  init_point = int( randnum(1) * countindx(prc_nprocs+1) ) + 1
1716  do ipp = 1, prc_nprocs
1717  if ( countindx(ipp+1) >= init_point ) then
1718  rank_initpoint = ipp-1
1719  exit
1720  end if
1721  end do
1722  grid_initpoint = init_point - countindx(rank_initpoint+1)
1723  ibuf(1) = rank_initpoint
1724  ibuf(2) = grid_initpoint
1725  endif
1726 
1727  call comm_bcast( ibuf, 2 )
1728 
1729  rank_initpoint = ibuf(1)
1730  grid_initpoint = ibuf(2)
1731 
1732 
1733  !--- Propagate lightning
1734 
1735  !$omp parallel do
1736  do j = js, je
1737  do i = is, ie
1738  do k = ks, ke
1739  l_path(k,i,j) = 0.0_rp !--- +-1 -> path already passed, +-2 -> path calculate current step
1740  flg_num_end(k,i,j,:) = 0
1741  end do
1742  end do
1743  end do
1744  wild_flag = 0
1745  end_grid(:) = 0
1746  comm_target = -999
1747  do iq = 1, 2 !--- loop for direction (1-> parallel, 2-> anti-parallel)
1748  if( iq == 1 ) then
1749  pm = 1.0_rp
1750  elseif( iq == 2 ) then
1751  pm = - 1.0_rp
1752  endif
1753  stop_flag = 0
1754  current_prc = rank_initpoint
1755 
1756  !---- determine initiation point
1757  if( rank_initpoint == prc_myrank ) then
1758  i = eovid(1,grid_initpoint)
1759  j = eovid(2,grid_initpoint)
1760  k = eovid(3,grid_initpoint)
1761  l_path(k,i,j) = pm
1762  if( iq == 1 ) fls_int_p(k,i,j) = fls_int_p(k,i,j) + 1.0_rp
1763  sign_flash = 2 + int( sign( 1.0_rp, qcrg(k,i,j) ) )
1764  else
1765  i = 0
1766  j = 0
1767  k = 0
1768  endif
1769 
1770  loop_path : do ipp = 1, kijmaxg !--- loop for path
1771  comm_flag = 0
1772  !---- calculate path of lightning
1773  if( current_prc == prc_myrank ) then
1774 
1775  !--- determine direction of path
1776  e_x = abs( efield(k,i,j,i_lt_x) )/cdx(i)
1777  e_y = abs( efield(k,i,j,i_lt_y) )/cdy(j)
1778  e_z = abs( efield(k,i,j,i_lt_z) )/cdz(k)
1779  e_det = max( e_x,e_y,e_z )
1780 
1781  i1 = i
1782  j1 = j
1783  k1 = k
1784  if( e_det == e_x ) then
1785  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_x) ) * pm )
1786  i1 = i+direction
1787  elseif( e_det == e_y ) then
1788  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_y) ) * pm )
1789  j1 = j+direction
1790  elseif( e_det == e_z ) then
1791  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_z) ) * pm )
1792  k1 = k+direction
1793  endif
1794 
1795 
1796  if( efield(k1,i1,j1,i_lt_abs) <= estop ) then
1797  !--- stop if |E| < Estop
1798  phi_end(iq) = qcrg(k,i,j)
1799  wild_flag = 1
1800  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1801  flg_num_end(k,i,j,sign_flash) = 1
1802  l_path(k,i,j) = pm
1803  stop_flag = 1
1804  end_grid(1) = i
1805  end_grid(2) = j
1806  end_grid(3) = k
1807  end_grid(4) = prc_myrank
1808  if( qhyd(k,i,j) < eps ) then
1809  corr_flag(iq) = 0
1810  else
1811  corr_flag(iq) = 1
1812  endif
1813  elseif( efield(k1,i1,j1,i_lt_abs) > estop ) then
1814  !--- propagate lightning path
1815  l_path(k, i ,j ) = pm
1816  l_path(k1,i1,j1) = pm
1817  endif
1818 
1819  !--- check whether lightning path reach top or bottom layer
1820  if( k1 == ke ) then
1821  !--- reach at top
1822  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1823  flg_num_end(k,i,j,sign_flash) = 1
1824  l_path(k,i,j) = pm
1825  stop_flag = 1
1826  end_grid(1) = i
1827  end_grid(2) = j
1828  end_grid(3) = k
1829  end_grid(4) = prc_myrank
1830  if( qhyd(k,i,j) < eps ) then
1831  corr_flag(iq) = 0
1832  else
1833  corr_flag(iq) = 1
1834  endif
1835  elseif( cz(k1) < zcg ) then
1836  !--- reach at groud
1837  num_end(k,i,j,2) = num_end(k,i,j,2) + 1.0_rp
1838  flg_num_end(k,i,j,2) = 1
1839  l_path(k,i,j) = pm
1840  stop_flag = 1
1841  end_grid(1) = i
1842  end_grid(2) = j
1843  end_grid(3) = k
1844  end_grid(4) = prc_myrank
1845  if( qhyd(k,i,j) < eps ) then
1846  corr_flag(iq) = 0
1847  else
1848  corr_flag(iq) = 1
1849  endif
1850  endif
1851 
1852  if( stop_flag /= 1 ) then
1853  !--- check whether lightning path reachs boundary in i-direction
1854  if( i1 == is-1 .and. .not. prc_has_w ) then
1855  !--- reach west boundary of global domain(stop)
1856  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1857  flg_num_end(k,i,j,sign_flash) = 1
1858  l_path(k,i,j) = pm
1859  l_path(k1,i1,j1) = pm
1860  stop_flag = 1
1861  end_grid(1) = i
1862  end_grid(2) = j
1863  end_grid(3) = k
1864  end_grid(4) = prc_myrank
1865  if( qhyd(k,i,j) < eps ) then
1866  corr_flag(iq) = 0
1867  else
1868  corr_flag(iq) = 1
1869  endif
1870  elseif( i1 == is-1 .and. prc_has_w ) then
1871  !--- reach west boundary of local domain(propagate)
1872  l_path(k,i,j) = pm
1873  l_path(k1,i1,j1) = pm
1874  comm_flag = 1
1875  comm_target = prc_next(prc_w)
1876  i1 = ie
1877  elseif( i1 == ie+1 .and. .not. prc_has_e ) then
1878  !--- reach east boundary of global domain(stop)
1879  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1880  flg_num_end(k,i,j,sign_flash) = 1
1881  l_path(k,i,j) = pm
1882  l_path(k1,i1,j1) = pm
1883  stop_flag = 1
1884  end_grid(1) = i
1885  end_grid(2) = j
1886  end_grid(3) = k
1887  end_grid(4) = prc_myrank
1888  if( qhyd(k,i,j) < eps ) then
1889  corr_flag(iq) = 0
1890  else
1891  corr_flag(iq) = 1
1892  endif
1893  elseif( i1 == ie+1 .and. prc_has_e ) then
1894  !--- reach east boundary of local domain(propagate)
1895  l_path(k,i,j) = pm
1896  l_path(k1,i1,j1) = pm
1897  comm_flag = 1
1898  comm_target = prc_next(prc_e)
1899  i1 = is
1900  endif
1901 
1902  !--- check whether lightning path reachs boundary in i-direction
1903  if( j1 == js-1 .and. .not. prc_has_s ) then
1904  !--- reach south boundary of global domain(stop)
1905  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1906  flg_num_end(k,i,j,sign_flash) = 1
1907  l_path(k,i,j) = pm
1908  l_path(k1,i1,j1) = pm
1909  stop_flag = 1
1910  end_grid(1) = i
1911  end_grid(2) = j
1912  end_grid(3) = k
1913  end_grid(4) = prc_myrank
1914  if( qhyd(k,i,j) < eps ) then
1915  corr_flag(iq) = 0
1916  else
1917  corr_flag(iq) = 1
1918  endif
1919  elseif( j1 == js-1 .and. prc_has_s ) then
1920  !--- reach south boundary of local domain(propagate)
1921  l_path(k,i,j) = pm
1922  l_path(k1,i1,j1) = pm
1923  comm_flag = 1
1924  comm_target = prc_next(prc_s)
1925  j1 = je
1926  elseif( j1 == je+1 .and. .not. prc_has_n ) then
1927  !--- reach north boundary of global domain(stop)
1928  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
1929  flg_num_end(k,i,j,sign_flash) = 1
1930  l_path(k,i,j) = pm
1931  l_path(k1,i1,j1) = pm
1932  comm_flag = 1
1933  stop_flag = 1
1934  end_grid(1) = i
1935  end_grid(2) = j
1936  end_grid(3) = k
1937  end_grid(4) = prc_myrank
1938  if( qhyd(k,i,j) < eps ) then
1939  corr_flag(iq) = 0
1940  else
1941  corr_flag(iq) = 1
1942  endif
1943  elseif( j1 == je+1 .and. prc_has_n ) then
1944  !--- reach north boundary of local domain(propagate)
1945  l_path(k,i,j) = pm
1946  l_path(k1,i1,j1) = pm
1947  comm_flag = 1
1948  comm_target = prc_next(prc_n)
1949  j1 = js
1950  endif
1951  endif
1952 
1953  endif
1954 
1955  !---- determine stop or not
1956  call mpi_bcast(stop_flag, 1, mpi_integer, current_prc, comm_world, ierr)
1957 
1958  if( stop_flag == 1 ) then
1959  !---- send flag wildfire
1960  call mpi_bcast(wild_flag, 1, mpi_integer, current_prc, comm_world, ierr)
1961 
1962  if ( prc_myrank == current_prc ) then
1963  ibuf(1:4) = end_grid(:)
1964  ibuf(5:6) = corr_flag(:)
1965  end if
1966  call mpi_bcast(ibuf, 6, mpi_integer, current_prc, comm_world, ierr)
1967  end_grid(:) = ibuf(1:4)
1968  corr_flag(:) = ibuf(5:6)
1969 
1970  stop_flag = 0
1971  !---- If lightning path reaches end_grid stop
1972  exit loop_path
1973  endif
1974 
1975  !---- determine wether process change occurs or not
1976  call mpi_bcast(comm_flag, 1, mpi_integer, current_prc, comm_world, ierr)
1977 
1978  !--- process change occurs
1979  if( comm_flag == 1 ) then
1980  call mpi_bcast(comm_target, 1, mpi_integer, current_prc, comm_world, ierr)
1981 
1982  if ( prc_myrank == comm_target ) then
1983  call mpi_recv(ibuf, 6, mpi_integer, current_prc, 1, comm_world, status, ierr)
1984  k1 = ibuf(1)
1985  i1 = ibuf(2)
1986  j1 = ibuf(3)
1987  corr_flag(:) = ibuf(4:5)
1988  sign_flash = ibuf(6)
1989  end if
1990 
1991  if ( prc_myrank == current_prc ) then
1992  ibuf(1) = k1
1993  ibuf(2) = i1
1994  ibuf(3) = j1
1995  ibuf(4:5) = corr_flag(:)
1996  ibuf(6) = sign_flash
1997  call mpi_send(ibuf, 6, mpi_integer, comm_target, 1, comm_world, ierr)
1998  end if
1999 
2000  !--- change current proc
2001  current_prc = comm_target
2002  endif
2003 
2004  if( current_prc == prc_myrank ) then
2005  k = k1
2006  j = j1
2007  i = i1
2008  endif
2009 
2010  enddo loop_path
2011 
2012 ! call PRC_MPIbarrier
2013 
2014  !---- Wildfire method
2015  if( wild_flag == 1 .and. end_grid(4) == prc_myrank ) then
2016  i1 = end_grid(1)
2017  j1 = end_grid(2)
2018  k1 = end_grid(3)
2019  i2 = end_grid(1)+1
2020  j2 = end_grid(2)
2021  k2 = end_grid(3)
2022  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2023  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2024  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2025  endif
2026  i1 = end_grid(1)
2027  j1 = end_grid(2)
2028  k1 = end_grid(3)
2029  i2 = end_grid(1)-1
2030  j2 = end_grid(2)
2031  k2 = end_grid(3)
2032  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2033  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2034  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2035  endif
2036  i1 = end_grid(1)
2037  j1 = end_grid(2)
2038  k1 = end_grid(3)
2039  i2 = end_grid(1)
2040  j2 = end_grid(2)+1
2041  k2 = end_grid(3)
2042  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2043  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2044  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2045  endif
2046  i1 = end_grid(1)
2047  j1 = end_grid(2)
2048  k1 = end_grid(3)
2049  i2 = end_grid(1)
2050  j2 = end_grid(2)-1
2051  k2 = end_grid(3)
2052  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2053  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2054  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2055  endif
2056  i1 = end_grid(1)
2057  j1 = end_grid(2)
2058  k1 = end_grid(3)
2059  i2 = end_grid(1)
2060  j2 = end_grid(2)
2061  k2 = end_grid(3)+1
2062  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2063  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2064  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2065  endif
2066  i1 = end_grid(1)
2067  j1 = end_grid(2)
2068  k1 = end_grid(3)
2069  i2 = end_grid(1)
2070  j2 = end_grid(2)
2071  k2 = end_grid(3)-1
2072  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
2073  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
2074  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
2075  endif
2076  endif
2077 
2078 ! call PRC_MPIbarrier
2079 
2080  enddo !--- loop for direction
2081 
2082 ! call COMM_vars8( L_path,1 )
2083 ! call COMM_wait ( L_path,1 )
2084 
2085  !---- Neutralization
2086  npls = 0
2087  nmns = 0
2088  sumdqrho_pls = 0.0_rp
2089  sumdqrho_mns = 0.0_rp
2090  !$omp parallel do reduction(+:Npls,Nmns,sumdqrho_pls,sumdqrho_mns)
2091  do j = js, je
2092  do i = is, ie
2093  do k = ks, ke
2094 
2095  !---- whether the grid is on lightning path or not
2096  if( l_path(k,i,j) /= 0.0_rp .and. abs( qcrg(k,i,j) ) > qrho_chan ) then
2097  if ( qcrg(k,i,j) >= 0.0_rp ) then
2098  npls = npls + 1
2099  dqrho_pls(k,i,j) = fp * ( abs( qcrg(k,i,j) )-qrho_neut )
2100  sumdqrho_pls = sumdqrho_pls &
2101  + fp * ( abs( qcrg(k,i,j) )-qrho_neut )
2102  dqrho_mns(k,i,j) = 0.0_rp
2103  else
2104  nmns = nmns + 1
2105  dqrho_mns(k,i,j) = fp * ( abs( qcrg(k,i,j) )-qrho_neut )
2106  sumdqrho_mns = sumdqrho_mns &
2107  + fp * ( abs( qcrg(k,i,j) )-qrho_neut )
2108  dqrho_pls(k,i,j) = 0.0_rp
2109  endif
2110  else
2111  dqrho_pls(k,i,j) = 0.0_rp
2112  dqrho_mns(k,i,j) = 0.0_rp
2113  endif
2114 
2115  enddo
2116  enddo
2117  enddo
2118 
2119  iprod(1) = npls
2120  iprod(2) = nmns
2121  call mpi_allreduce(iprod, ibuf, 2, mpi_integer, mpi_sum, comm_world, ierr)
2122  npls = ibuf(1)
2123  nmns = ibuf(2)
2124 
2125  ntot = npls + nmns
2126 
2127  if ( ntot > 0 ) exit
2128 
2129  end do
2130 
2131  if ( count_path > nutr_itmax ) then
2132  log_info("ATMOS_PHY_LT_Efield",*) "Reach limit iteration for searching flash path", count_path, npls, nmns, current_prc
2133  endif
2134 
2135  if( corr_flag(1) == 1 .and. corr_flag(2) == 1 ) then
2136  dqrho_cor = sumdqrho_pls - sumdqrho_mns
2137  rprod2 = dqrho_cor
2138  call mpi_allreduce(rprod2, rbuf, 1, comm_datatype, mpi_sum, comm_world, ierr)
2139  dqrho_cor = rbuf
2140 
2141  dqrho_cor = dqrho_cor / ntot !---- Eq.(2) of MacGorman et al. (2001)
2142  else
2143  dqrho_cor = 0.0_rp !---- No correlation if flash to ground and out of cloud
2144  endif
2145 
2146  !$omp parallel do
2147  do j = js, je
2148  do i = is, ie
2149  do k = ks, ke
2150  d_qcrg(k,i,j) = 0.0_rp
2151  end do
2152  end do
2153  end do
2154  rprod1 = 0.0_rp
2155  rprod2 = 0.0_rp
2156  !--- pseud-2D experiment for x-direction
2157  if( prc_num_x == 1 .and. imax == 2 ) then
2158 
2159  !$omp parallel do
2160  do j = js, je
2161  do k = ks, ke
2162  if( l_path(k,is,j) /= 0.0_rp ) then
2163  l_path(k,ie,j) = l_path(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
2164  if( dqrho_pls(k,is,j) /= 0.0_rp ) then
2165  d_qcrg(k,is,j) = - dqrho_pls(k,is,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2166  d_qcrg(k,ie,j) = d_qcrg(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
2167  elseif( dqrho_mns(k,is,j) /= 0.0_rp ) then
2168  d_qcrg(k,is,j) = dqrho_mns(k,is,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2169  d_qcrg(k,ie,j) = d_qcrg(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
2170  endif
2171  elseif( l_path(k,ie,j) /= 0.0_rp ) then
2172  l_path(k,is,j) = l_path(k,ie,j) !--- copy IE value to IS (which is regarded as same grid in psued-2D)
2173  if( dqrho_pls(k,ie,j) /= 0.0_rp ) then
2174  d_qcrg(k,ie,j) = - dqrho_pls(k,ie,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2175  d_qcrg(k,is,j) = d_qcrg(k,ie,j) !--- copy IS value to IS (which is regarded as same grid in psued-2D)
2176  elseif( dqrho_mns(k,ie,j) /= 0.0_rp ) then
2177  d_qcrg(k,ie,j) = dqrho_mns(k,ie,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2178  d_qcrg(k,is,j) = d_qcrg(k,ie,j) !--- copy IS value to IS (which is regarded as same grid in psued-2D)
2179  endif
2180  endif
2181 
2182  do iq = 1, 3
2183  if( flg_num_end(k,is,j,iq) == 1 ) then
2184  num_end(k,ie,j,iq) = num_end(k,ie,j,iq) + 1.0_rp
2185  endif
2186  if( flg_num_end(k,ie,j,iq) == 1 ) then
2187  num_end(k,is,j,iq) = num_end(k,is,j,iq) + 1.0_rp
2188  endif
2189  enddo
2190 
2191  enddo
2192  enddo
2193 
2194  !--- pseud-2D experiment for y-direction
2195  elseif( prc_num_y == 1 .and. jmax == 2 ) then
2196 
2197  !$omp parallel do
2198  do i = is, ie
2199  do k = ks, ke
2200  if( l_path(k,i,js) /= 0.0_rp ) then
2201  l_path(k,i,je) = l_path(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
2202  if( dqrho_pls(k,i,js) /= 0.0_rp ) then
2203  d_qcrg(k,i,js) = - dqrho_pls(k,i,js) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2204  d_qcrg(k,i,je) = d_qcrg(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
2205  elseif( dqrho_mns(k,i,js) /= 0.0_rp ) then
2206  d_qcrg(k,i,js) = dqrho_mns(k,i,js) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2207  d_qcrg(k,i,je) = d_qcrg(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
2208  endif
2209  elseif( l_path(k,i,je) /= 0.0_rp ) then
2210  l_path(k,i,js) = l_path(k,i,je) !--- copy JE value to JS (which is regarded as same grid in psued-2D)
2211  if( dqrho_pls(k,i,je) /= 0.0_rp ) then
2212  d_qcrg(k,i,je) = - dqrho_pls(k,i,je) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2213  d_qcrg(k,i,js) = d_qcrg(k,i,je) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
2214  elseif( dqrho_mns(k,i,je) /= 0.0_rp ) then
2215  d_qcrg(k,i,je) = dqrho_mns(k,i,je) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2216  d_qcrg(k,i,js) = d_qcrg(k,i,je) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
2217  endif
2218  endif
2219 
2220  do iq = 1, 3
2221  if( flg_num_end(k,i,js,iq) == 1 ) then
2222  num_end(k,i,je,iq) = num_end(k,i,je,iq) + 1.0_rp
2223  endif
2224  if( flg_num_end(k,i,je,iq) == 1 ) then
2225  num_end(k,i,js,iq) = num_end(k,i,js,iq) + 1.0_rp
2226  endif
2227  enddo
2228 
2229  enddo
2230  enddo
2231 
2232  else
2233 
2234  !$omp parallel do
2235  do j = js, je
2236  do i = is, ie
2237  do k = ks, ke
2238  if( l_path(k,i,j) /= 0.0_rp .and. dqrho_pls(k,i,j) /= 0.0_rp ) then
2239  d_qcrg(k,i,j) = - dqrho_pls(k,i,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2240  elseif( l_path(k,i,j) /= 0.0_rp .and. dqrho_mns(k,i,j) /= 0.0_rp ) then
2241  d_qcrg(k,i,j) = dqrho_mns(k,i,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
2242  endif
2243  enddo
2244  enddo
2245  enddo
2246 
2247  endif
2248 
2249  !$omp parallel do
2250  do j = js, je
2251  do i = is, ie
2252  do k = ks, ke
2253  lt_path(k,i,j) = lt_path(k,i,j) + l_path(k,i,j)
2254  enddo
2255  enddo
2256  enddo
2257 
2258  call prof_rapend('LT_neut_MG2001', 1)
2259 
2260  return
2261 
2262  end subroutine atmos_phy_lt_neutralization_mg2001
2263  !-----------------------------------------------------------------------------
2265  !-----------------------------------------------------------------------------
2266  subroutine atmos_phy_lt_neutralization_f2013( &
2267  KA, KS, KE, & ! [IN]
2268  IA, IS, IE, & ! [IN]
2269  JA, JS, JE, & ! [IN]
2270  KIJMAX, & ! [IN]
2271  Efield, & ! [IN]
2272  E_pot, & ! [IN]
2273  DENS, & ! [IN]
2274  QCRG, & ! [IN]
2275  QHYD, & ! [IN]
2276  NUM_end, & ! [INOUT]
2277  LT_path, & ! [INOUT]
2278  fls_int_p, & ! [OUT]
2279  d_QCRG ) ! [OUT]
2280  use scale_const, only: &
2281  eps => const_eps, &
2282  epsvac => const_epsvac, &
2283  epsair => const_epsair
2284  use scale_atmos_grid_cartesc, only: &
2285  cz => atmos_grid_cartesc_cz, &
2286  cx => atmos_grid_cartesc_cx, &
2287  cy => atmos_grid_cartesc_cy, &
2288  cdz => atmos_grid_cartesc_cdz, &
2289  cdx => atmos_grid_cartesc_cdx, &
2290  cdy => atmos_grid_cartesc_cdy, &
2291  rfdz => atmos_grid_cartesc_rfdz, &
2292  rfdx => atmos_grid_cartesc_rfdx, &
2293  rfdy => atmos_grid_cartesc_rfdy
2294  use scale_prc_cartesc, only: &
2295  prc_has_w, &
2296  prc_has_e, &
2297  prc_has_n, &
2298  prc_has_s, &
2299  prc_w, &
2300  prc_n, &
2301  prc_e, &
2302  prc_s, &
2303  prc_nw, &
2304  prc_ne, &
2305  prc_sw, &
2306  prc_se, &
2307  prc_next, &
2308  prc_num_x,&
2309  prc_num_y
2310  use scale_prc, only: &
2311  prc_abort, &
2312  prc_ismaster, &
2313  prc_nprocs, &
2314  prc_masterrank, &
2315  prc_mpibarrier, &
2316  prc_myrank
2317  use scale_comm_cartesc, only: &
2318  comm_datatype, &
2319  comm_world
2320 ! use scale_random, only: &
2321 ! RANDOM_reset, &
2322 ! RANDOM_uniform
2323  implicit none
2324 
2325  integer, intent(in) :: KA, KS, KE
2326  integer, intent(in) :: IA, IS, IE
2327  integer, intent(in) :: JA, JS, JE
2328  integer, intent(in) :: KIJMAX
2329  real(RP), intent(in) :: Efield (KA,IA,JA,4) !-- Electric field [V/m] (4->|E|)
2330  real(RP), intent(in) :: E_pot (KA,IA,JA) !-- Electric potential [V]
2331  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
2332  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
2333  real(RP), intent(in) :: QHYD (KA,IA,JA) !-- Total hydrometeor mixing ratio
2334  real(RP), intent(inout) :: NUM_end (KA,IA,JA,3) !-- Number of Lightning stop grid
2335  real(RP), intent(inout) :: LT_path (KA,IA,JA) !-- Number of path
2336  real(RP), intent(out) :: fls_int_p(KA,IA,JA) !-- Flash initiation point (0->no flash, 1->flash start at the point)
2337  real(RP), intent(out) :: d_QCRG (KA,IA,JA) !-- Charge density [nC/m3]
2338 
2339  real(RP), parameter :: q_thre = 0.1_rp ! threshold of discharge zone (Fierro et al. 2013) [nC/m3]
2340 
2341  real(RP) :: B(IA,JA) !--- B in Fierro et al. 2013
2342  integer :: C(IA,JA) !--- flg for cylinder (in cylinder -> C=1)
2343  real(RP) :: Q_d, Fpls, Fmns
2344  real(RP) :: Spls, Smns
2345  real(RP) :: Spls_g, Smns_g
2346  real(RP) :: distance
2347  real(RP) :: Edif(KS:KE), Edif_max, abs_qcrg_max
2348  real(RP) :: sw
2349 
2350  integer, allocatable :: proc_num(:), proc_numg(:)
2351  real(RP),allocatable :: E_exce_x(:), E_exce_x_g(:) !--- x point of column in which |E|>Eint is included [m] (_g means global attribute)
2352  real(RP),allocatable :: E_exce_y(:), E_exce_y_g(:) !--- y point of column in which |E|>Eint is included [m] (_g means global attribute)
2353  real(RP) :: exce_grid(2,KIJMAX)
2354  integer :: count1(PRC_nprocs)
2355  integer :: countindx(PRC_nprocs+1)
2356  integer :: num_own !--- number of column whose |E| > Eint-dEint for each process
2357  integer :: num_total !--- total number of column whose |E| > Eint-dEint
2358  real(RP) :: rbuf1(2), rbuf2(2)
2359  integer :: k, i, j, ipp, iq, ierr
2360 
2361  call prof_rapstart('LT_neut_F2013', 1)
2362 
2363  !$omp parallel do
2364  do j = js, je
2365  do i = is, ie
2366  do k = ks, ie
2367  num_end(k,i,j,:) = 0.0_rp
2368  end do
2369  end do
2370  end do
2371  !--- search grid whose Efield is over threshold of flash inititation
2372  num_own = 0
2373  do j = js, je
2374  do i = is, ie
2375  edif_max = -1.0_rp
2376  do k = ks, ke
2377  edif(k) = efield(k,i,j,i_lt_abs) - ( eint-deleint )
2378  edif_max = max( edif_max, edif(k) )
2379  enddo
2380  if ( edif_max > 0.0_rp ) then
2381  num_own = num_own + 1
2382  exce_grid(1,num_own) = cx(i)
2383  exce_grid(2,num_own) = cy(j)
2384  do k = ks, ke
2385  fls_int_p(k,i,j) = 0.5_rp + sign( 0.5_rp, edif(k)-eps )
2386  enddo
2387  else
2388  do k = ks, ke
2389  fls_int_p(k,i,j) = 0.0_rp
2390  enddo
2391  endif
2392  enddo
2393  enddo
2394 
2395  !**** proc_num(0~) -> process number of each grid with |E|> E_threthold (local)
2396 ! allocate(proc_num(own_prc_total))
2397  allocate(e_exce_x(num_own))
2398  allocate(e_exce_y(num_own))
2399 
2400 ! proc_num(:) = PRC_myrank
2401  e_exce_x(1:num_own) = exce_grid(1,1:num_own)
2402  e_exce_y(1:num_own) = exce_grid(2,1:num_own)
2403 
2404  call mpi_allgather( num_own, 1, mpi_integer, &
2405  count1, 1, mpi_integer, &
2406  comm_world, ierr )
2407 
2408  countindx(1) = 0
2409  do ipp = 1, prc_nprocs
2410  countindx(ipp+1) = countindx(ipp) + count1(ipp)
2411  enddo
2412 
2413  !---- Create global version of proc_num(proc_numg)
2414  !**** countindx(PROC_nprocs) -> total number of grid with |E|>E_threthold
2415  !**** proc_numg(0~) -> process number of each grid with |E|> E_threthold (global)
2416  num_total = countindx(prc_nprocs+1)
2417  allocate(e_exce_x_g(num_total))
2418  allocate(e_exce_y_g(num_total))
2419 
2420  call mpi_allgatherv( e_exce_x, num_own, comm_datatype, &
2421  e_exce_x_g, count1, countindx, comm_datatype, &
2422  comm_world, ierr )
2423 
2424  call mpi_allgatherv( e_exce_y, num_own, comm_datatype, &
2425  e_exce_y_g, count1, countindx, comm_datatype, &
2426  comm_world, ierr )
2427 
2428  c(:,:) = 0
2429  do ipp = 1, num_total
2430  do j = js, je
2431  do i = is, ie
2432  distance = sqrt( ( cx(i)-e_exce_x_g(ipp) )**2 &
2433  + ( cy(j)-e_exce_y_g(ipp) )**2 )
2434  if( distance <= r_neut ) then
2435  c(i,j) = 1
2436  endif
2437  enddo
2438  enddo
2439  enddo
2440 
2441  spls = 0.0_rp
2442  smns = 0.0_rp
2443  !$omp parallel do reduction(+:Spls,Smns) &
2444  !$omp private(abs_qcrg_max,sw)
2445  do j = js, je
2446  do i = is, ie
2447 
2448  b(i,j) = 0.0_rp
2449  if ( c(i,j) == 1 ) then
2450  abs_qcrg_max = abs( qcrg(ks,i,j) )
2451  do k = ks+1, ke
2452  abs_qcrg_max = max( abs_qcrg_max, abs( qcrg(k,i,j) ) )
2453  end do
2454  if ( abs_qcrg_max >= q_thre ) then
2455  b(i,j) = 1.0_rp
2456  do k = ks, ke
2457  sw = 0.5_rp + sign( 0.5_rp, qcrg(k,i,j) )
2458  spls = spls + qcrg(k,i,j) * sw
2459  smns = smns - qcrg(k,i,j) * ( 1.0_rp - sw )
2460  enddo
2461  end if
2462  end if
2463 
2464  enddo
2465  enddo
2466 
2467  rbuf1(1) = spls
2468  call mpi_allreduce( rbuf1, rbuf2, 1, comm_datatype, &
2469  mpi_sum, comm_world, ierr )
2470  spls_g = rbuf2(1)
2471 
2472  rbuf1(1) = smns
2473  call mpi_allreduce( rbuf1, rbuf2, 1, comm_datatype, &
2474  mpi_sum, comm_world, ierr )
2475  smns_g = rbuf2(1)
2476 
2477  if( max( spls_g,smns_g )*0.3_rp < min( spls_g,smns_g ) ) then
2478  q_d = 0.3_rp * max( spls_g,smns_g )
2479  else
2480  q_d = min( spls_g,smns_g )
2481  endif
2482 
2483  if( spls_g /= 0.0_rp ) then
2484  fpls = q_d / spls_g
2485  else
2486  fpls = 0.0_rp
2487  endif
2488 
2489  if( smns_g /= 0.0_rp ) then
2490  fmns = q_d / smns_g
2491  else
2492  fmns = 0.0_rp
2493  endif
2494 
2495  !---- select initial point of flash by random select
2496  !$omp parallel do
2497  do j = js, je
2498  do i = is, ie
2499  do k = ks, ke
2500  if( qcrg(k,i,j) > q_thre ) then
2501  d_qcrg(k,i,j) = - fpls * ( qcrg(k,i,j) - q_thre ) * b(i,j)
2502  lt_path(k,i,j) = lt_path(k,i,j) + b(i,j)
2503  elseif( qcrg(k,i,j) < -q_thre ) then
2504  d_qcrg(k,i,j) = + fmns * ( - qcrg(k,i,j) - q_thre ) * b(i,j)
2505  lt_path(k,i,j) = lt_path(k,i,j) - b(i,j)
2506  else
2507  d_qcrg(k,i,j) = 0.0_rp
2508  endif
2509  enddo
2510  enddo
2511  enddo
2512 
2513 
2514  deallocate(e_exce_x)
2515  deallocate(e_exce_y)
2516  deallocate(e_exce_x_g)
2517  deallocate(e_exce_y_g)
2518 
2519 
2520  call prof_rapend('LT_neut_F2013', 1)
2521 
2522  return
2523  end subroutine atmos_phy_lt_neutralization_f2013
2524  !-----------------------------------------------------------------------------
2526  !-----------------------------------------------------------------------------
2527  subroutine atmos_phy_lt_judge_abse( &
2528  KA, KS, KE, & ! [IN]
2529  IA, IS, IE, & ! [IN]
2530  JA, JS, JE, & ! [IN]
2531  DENS, & ! [IN]
2532  Efield, & ! [IN]
2533  Emax, & ! [OUT]
2534  flg_lt_neut ) ! [OUT]
2535  use scale_comm_cartesc, only: &
2536  comm_datatype, &
2537  comm_world
2538  use scale_const, only: &
2539  large_num => const_huge
2540  implicit none
2541 
2542  integer, intent(in) :: KA, KS, KE
2543  integer, intent(in) :: IA, IS, IE
2544  integer, intent(in) :: JA, JS, JE
2545  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
2546  real(RP), intent(in) :: Efield (KA,IA,JA) !-- Absolute value of Electric field [V/m]
2547  real(RP), intent(out) :: Emax !-- Maximum value of Electric field [V/m]
2548  integer, intent(out) :: flg_lt_neut !-- flg 1-> neutralization, 0->no neutralization
2549 
2550  real(RP) :: Eint_hgt(KA,IA,JA)
2551  real(RP) :: E_det, rbuf, rprod1
2552  integer :: own_prc_total, iprod1, buf
2553  integer :: k, i, j, ierr
2554 
2555  !--- search grid whose Efield is over threshold of flash inititation
2556  own_prc_total = 0
2557  if( flg_eint_hgt == 1.0_rp ) then
2558  !$omp parallel do &
2559  !$omp reduction(+:own_prc_total) &
2560  !$omp private(E_det)
2561  do j = js, je
2562  do i = is, ie
2563  do k = ks, ke
2564  eint_hgt(k,i,j) = min( eint*dens(k,i,j)/rho0,eint )
2565  if( eint_hgt(k,i,j) < estop ) then
2566  eint_hgt(k,i,j) = large_num !--- ignore when Eint < Estop
2567  endif
2568  e_det = efield(k,i,j) - eint_hgt(k,i,j)
2569  if( e_det > 0.0_rp ) then
2570  own_prc_total = own_prc_total + 1
2571  endif
2572  enddo
2573  enddo
2574  enddo
2575  else
2576  !$omp parallel do &
2577  !$omp reduction(+:own_prc_total) &
2578  !$omp private(E_det)
2579  do j = js, je
2580  do i = is, ie
2581  do k = ks, ke
2582  e_det = efield(k,i,j) - ( eint-deleint )
2583  if( e_det > 0.0_rp ) then
2584  own_prc_total = own_prc_total + 1
2585  endif
2586  enddo
2587  enddo
2588  enddo
2589  endif
2590 
2591  !--- Add number of grids, whose |E| are over threshold of flash initiaion, for all process
2592  iprod1 = own_prc_total
2593  call mpi_allreduce(iprod1, buf, 1, mpi_integer, mpi_sum, comm_world, ierr)
2594  rprod1 = maxval(efield(ks:ke,is:ie,js:je))
2595  call mpi_allreduce(rprod1, rbuf, 1, comm_datatype, mpi_max, comm_world, ierr)
2596  !--- exit when no grid point with |E| over threshold of flash initiation exist
2597  emax = rbuf
2598  if( buf == 0 ) then
2599  flg_lt_neut = 0
2600  else
2601  flg_lt_neut = 1
2602  if( emax < eint .and. emax >= eint-deleint .and. flg_eint_hgt == 0.0_rp ) then
2603  flg_lt_neut = 2
2604  endif
2605  endif
2606 
2607  end subroutine atmos_phy_lt_judge_abse
2608  !-----------------------------------------------------------------------------
2610  !-----------------------------------------------------------------------------
2611  subroutine atmos_phy_lt_sato2019_select_dqcrg_from_lut( &
2612  KA, KS, KE, & ! [IN]
2613  IA, IS, IE, & ! [IN]
2614  JA, JS, JE, & ! [IN]
2615  NLIQ, & ! [IN]
2616  TEMP, & ! [IN]
2617  DENS, & ! [IN]
2618  QLIQ, & ! [IN]
2619  dqcrg, & ! [OUT]
2620  beta_crg ) ! [OUT]
2621  use scale_const, only: &
2622  t00 => const_tem00
2623  implicit none
2624 
2625  integer, intent(in) :: ka, ks, ke
2626  integer, intent(in) :: ia, is, ie
2627  integer, intent(in) :: ja, js, je
2628  integer, intent(in) :: nliq
2629  real(rp), intent(in) :: temp(ka,ia,ja)
2630  real(rp), intent(in) :: dens(ka,ia,ja)
2631  real(rp), intent(in) :: qliq(ka,ia,ja,nliq)
2632  real(rp), intent(out) :: dqcrg(ka,ia,ja)
2633  real(rp), intent(out) :: beta_crg(ka,ia,ja)
2634 
2635  integer :: i, j, k, pp, qq, iq
2636  integer :: grid(2)
2637  real(rp) :: cwc
2638  real(rp) :: diffx(nxlut_lt), diffy(nylut_lt)
2639 
2640  !$omp parallel do &
2641  !$omp private(cwc,diffx,diffy,grid)
2642  do j = js, je
2643  do i = is, ie
2644  do k = ks, ke
2645  if( temp(k,i,j) <= t00 .and. temp(k,i,j) >= tcrglimit ) then
2646  cwc = 0.0_rp
2647  do iq = 1, nliq
2648  cwc = cwc + qliq(k,i,j,iq) * dens(k,i,j) * 1.0e+3_rp ![g/m3]
2649  enddo
2650  do pp = 1, nxlut_lt
2651  diffx(pp) = abs( grid_lut_t(pp)-temp(k,i,j) )
2652  enddo
2653  grid(1) = minloc( diffx,1 )
2654  do qq = 1, nylut_lt
2655  diffy(qq) = abs( grid_lut_l(qq)-cwc )
2656  enddo
2657  grid(2) = minloc( diffy,1 )
2658  dqcrg(k,i,j) = dq_chrg( grid(1), grid(2) ) &
2659  *( 0.5_rp + sign( 0.5_rp,cwc-1.0e-2_rp ) ) !--- no charge separation when cwc < 0.01 [g/m3]
2660  else
2661  dqcrg(k,i,j) = 0.0_rp
2662  endif
2663  if( temp(k,i,j) >= -30.0_rp+t00 ) then
2664  beta_crg(k,i,j) = 1.0_rp
2665  elseif( temp(k,i,j) < -30.0_rp+t00 .and. temp(k,i,j) >= -43.0_rp+t00 ) then
2666  beta_crg(k,i,j) = 1.0_rp - ( ( temp(k,i,j)-t00+30.0_rp )/13.0_rp )**2
2667 ! elseif( TEMP(k,i,j) < -43.0_RP+T00 ) then
2668  else
2669  beta_crg(k,i,j) = 0.0_rp
2670  endif
2671  enddo
2672  enddo
2673  enddo
2674 
2675  return
2676 
2677  end subroutine atmos_phy_lt_sato2019_select_dqcrg_from_lut
2678  !-----------------------------------------------------------------------------
2679 end module scale_atmos_phy_lt_sato2019
scale_atmos_grid_cartesc_index::i_uy
integer, public i_uy
Definition: scale_atmos_grid_cartesC_index.F90:99
scale_comm_cartesc::comm_datatype
integer, public comm_datatype
datatype of variable
Definition: scale_comm_cartesC.F90:98
scale_atmos_grid_cartesc_index::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:100
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_prc_cartesc::prc_n
integer, parameter, public prc_n
[node direction] north
Definition: scale_prc_cartesC.F90:33
scale_atmos_phy_lt_sato2019::atmos_phy_lt_sato2019_setup
subroutine, public atmos_phy_lt_sato2019_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, IMAXG, JMAXG, KMAX, MP_TYPE, CDX, CDY)
Setup.
Definition: scale_atmos_phy_lt_sato2019.F90:150
scale_atmos_grid_cartesc::atmos_grid_cartesc_cdz
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cdz
z-length of control volume [m]
Definition: scale_atmos_grid_cartesC.F90:42
scale_atmos_phy_lt_sato2019::atmos_phy_lt_electric_field
subroutine atmos_phy_lt_electric_field(KA, KS, KE, IA, IS, IE, JA, JS, JE, QCRG, DENS, RHOT, E_pot, Efield)
calculate electric field from charge density of each grid temporaly Bi-CGSTAB is used in this compone...
Definition: scale_atmos_phy_lt_sato2019.F90:921
scale_prc_cartesc::prc_sw
integer, parameter, public prc_sw
[node direction] southwest
Definition: scale_prc_cartesC.F90:38
scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdx
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rfdx
reciprocal of face-dx
Definition: scale_atmos_grid_cartesC.F90:67
scale_prc_cartesc::prc_next
integer, dimension(8), public prc_next
node ID of 8 neighbour process
Definition: scale_prc_cartesC.F90:45
scale_prc_cartesc::prc_has_s
logical, public prc_has_s
Definition: scale_prc_cartesC.F90:50
scale_const::const_epsvac
real(rp), parameter, public const_epsvac
parts par million
Definition: scale_const.F90:93
scale_atmos_phy_lt_sato2019::atmos_phy_lt_sato2019_adjustment
subroutine, public atmos_phy_lt_sato2019_adjustment(KA, KS, KE, IA, IS, IE, JA, JS, JE, KIJMAX, IMAX, JMAX, QA_LT, DENS, RHOT, QHYD, Sarea, dt_LT, QTRC, Epot)
Update of charge density.
Definition: scale_atmos_phy_lt_sato2019.F90:304
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_phy_lt_sato2019::atmos_phy_lt_judge_abse
subroutine atmos_phy_lt_judge_abse(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, Efield, Emax, flg_lt_neut)
judge max of |E| exceeds E_init
Definition: scale_atmos_phy_lt_sato2019.F90:2535
scale_atmos_grid_cartesc_index::i_uv
integer, public i_uv
Definition: scale_atmos_grid_cartesC_index.F90:101
scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdx
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rcdx
reciprocal of center-dx
Definition: scale_atmos_grid_cartesC.F90:65
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
scale_atmos_grid_cartesc_metric
module Atmosphere Grid CartesianC metirc
Definition: scale_atmos_grid_cartesC_metric.F90:12
scale_atmos_grid_cartesc::atmos_grid_cartesc_cdy
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cdy
y-length of control volume [m]
Definition: scale_atmos_grid_cartesC.F90:61
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:90
scale_random
module RANDOM
Definition: scale_random.F90:11
scale_io::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_io.F90:321
scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdz
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rcdz
reciprocal of center-dz
Definition: scale_atmos_grid_cartesC.F90:44
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_prc_cartesc::prc_se
integer, parameter, public prc_se
[node direction] southeast
Definition: scale_prc_cartesC.F90:39
scale_prc_cartesc::prc_has_n
logical, public prc_has_n
Definition: scale_prc_cartesC.F90:48
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:90
scale_prc_cartesc::prc_has_e
logical, public prc_has_e
Definition: scale_prc_cartesC.F90:49
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_mapf
real(rp), dimension(:,:,:,:), allocatable, public atmos_grid_cartesc_metric_mapf
map factor
Definition: scale_atmos_grid_cartesC_metric.F90:34
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:98
scale_prc::prc_mpibarrier
subroutine, public prc_mpibarrier
Barrier MPI.
Definition: scale_prc.F90:831
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdy
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rcdy
reciprocal of center-dy
Definition: scale_atmos_grid_cartesC.F90:66
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_prc::prc_masterrank
integer, parameter, public prc_masterrank
master process in each communicator
Definition: scale_prc.F90:66
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_gsqrt
real(rp), dimension(:,:,:,:), allocatable, public atmos_grid_cartesc_metric_gsqrt
transformation metrics from Z to Xi, {G}^1/2
Definition: scale_atmos_grid_cartesC_metric.F90:37
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_prc_cartesc::prc_w
integer, parameter, public prc_w
[node direction] west
Definition: scale_prc_cartesC.F90:32
scale_prc_cartesc
module process / cartesC
Definition: scale_prc_cartesC.F90:11
scale_const::const_huge
real(rp), public const_huge
huge number
Definition: scale_const.F90:35
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:92
scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_j23g
real(rp), dimension(:,:,:,:), allocatable, public atmos_grid_cartesc_metric_j23g
(2,3) element of Jacobian matrix * {G}^1/2
Definition: scale_atmos_grid_cartesC_metric.F90:39
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_comm_cartesc::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm_cartesC.F90:99
scale_prc_cartesc::prc_s
integer, parameter, public prc_s
[node direction] south
Definition: scale_prc_cartesC.F90:35
scale_atmos_phy_lt_sato2019::mul_matrix
subroutine mul_matrix(KA, KS, KE, IA, IS, IE, JA, JS, JE, V, M, C)
Definition: scale_atmos_phy_lt_sato2019.F90:1485
scale_prc::prc_nprocs
integer, public prc_nprocs
myrank in local communicator
Definition: scale_prc.F90:89
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:90
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_prc_cartesc::prc_nw
integer, parameter, public prc_nw
[node direction] northwest
Definition: scale_prc_cartesC.F90:36
scale_prc_cartesc::prc_num_y
integer, public prc_num_y
y length of 2D processor topology
Definition: scale_prc_cartesC.F90:42
scale_prc_cartesc::prc_e
integer, parameter, public prc_e
[node direction] east
Definition: scale_prc_cartesC.F90:34
scale_atmos_grid_cartesc::atmos_grid_cartesc_cy
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cy
center coordinate [m]: y, local
Definition: scale_atmos_grid_cartesC.F90:56
scale_const::const_epsair
real(rp), public const_epsair
parts par million
Definition: scale_const.F90:94
scale_atmos_phy_lt_sato2019::atmos_phy_lt_neutralization_f2013
subroutine atmos_phy_lt_neutralization_f2013(KA, KS, KE, IA, IS, IE, JA, JS, JE, KIJMAX, Efield, E_pot, DENS, QCRG, QHYD, NUM_end, LT_path, fls_int_p, d_QCRG)
calculate charge neutralization Fierro et al. (2013)
Definition: scale_atmos_phy_lt_sato2019.F90:2280
scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdz
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rfdz
reciprocal of face-dz
Definition: scale_atmos_grid_cartesC.F90:45
scale_atmos_phy_lt_sato2019
module atmosphere / physics / lightninh / SATO2019
Definition: scale_atmos_phy_lt_sato2019.F90:12
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:650
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_j33g
real(rp), public atmos_grid_cartesc_metric_j33g
(3,3) element of Jacobian matrix * {G}^1/2
Definition: scale_atmos_grid_cartesC_metric.F90:40
scale_prc_cartesc::prc_num_x
integer, public prc_num_x
x length of 2D processor topology
Definition: scale_prc_cartesC.F90:41
scale_atmos_grid_cartesc::atmos_grid_cartesc_cdx
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cdx
x-length of control volume [m]
Definition: scale_atmos_grid_cartesC.F90:60
scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_j13g
real(rp), dimension(:,:,:,:), allocatable, public atmos_grid_cartesc_metric_j13g
(1,3) element of Jacobian matrix * {G}^1/2
Definition: scale_atmos_grid_cartesC_metric.F90:38
scale_atmos_phy_lt_sato2019::atmos_phy_lt_sato2019_select_dqcrg_from_lut
subroutine, public atmos_phy_lt_sato2019_select_dqcrg_from_lut(KA, KS, KE, IA, IS, IE, JA, JS, JE, NLIQ, TEMP, DENS, QLIQ, dqcrg, beta_crg)
Select cwc-temp point on LUT.
Definition: scale_atmos_phy_lt_sato2019.F90:2621
scale_atmos_grid_cartesc
module atmosphere / grid / cartesC
Definition: scale_atmos_grid_cartesC.F90:12
scale_atmos_grid_cartesc::atmos_grid_cartesc_cz
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cz
center coordinate [m]: z, local
Definition: scale_atmos_grid_cartesC.F90:40
scale_atmos_grid_cartesc_index::i_uvz
integer, public i_uvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_atmos_phy_lt_sato2019::atmos_phy_lt_neutralization_mg2001
subroutine atmos_phy_lt_neutralization_mg2001(KA, KS, KE, IA, IS, IE, JA, JS, JE, KIJMAX, IMAX, JMAX, Efield, E_pot, DENS, QCRG, QHYD, NUM_end, LT_path, fls_int_p, d_QCRG)
calculate charge neutralization MacGorman et al. (2001)
Definition: scale_atmos_phy_lt_sato2019.F90:1578
scale_prc_cartesc::prc_ne
integer, parameter, public prc_ne
[node direction] northeast
Definition: scale_prc_cartesC.F90:37
scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdy
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rfdy
reciprocal of face-dy
Definition: scale_atmos_grid_cartesC.F90:68
scale_atmos_phy_lt_sato2019::atmos_phy_lt_solve_bicgstab
subroutine atmos_phy_lt_solve_bicgstab(KA, KS, KE, IA, IS, IE, JA, JS, JE, PHI_N, PHI, M, B)
Definition: scale_atmos_phy_lt_sato2019.F90:1228
scale_prc_cartesc::prc_has_w
logical, public prc_has_w
Definition: scale_prc_cartesC.F90:47
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_atmos_grid_cartesc::atmos_grid_cartesc_cx
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cx
center coordinate [m]: x, local
Definition: scale_atmos_grid_cartesC.F90:55
scale_prc::prc_ismaster
logical, public prc_ismaster
master process in local communicator?
Definition: scale_prc.F90:91
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:91
Generated by   doxygen 1.8.17