SCALE-RM
scale_atmos_phy_lt_sato2019.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
15 !-------------------------------------------------------------------------------
16 #ifndef COLORING
17 #ifdef _OPENACC
18 #define COLORING 2
19 #elif defined(_OPENMP)
20 #define COLORING 1
21 #else
22 #define COLORING 0
23 #endif
24 #endif
25 
26 #include "scalelib.h"
27 module scale_atmos_phy_lt_sato2019
28  !-----------------------------------------------------------------------------
29  !
30  !++ used modules
31  !
32  use mpi
33  use scale_precision
34  use scale_io
35  use scale_prof
36  !-----------------------------------------------------------------------------
37  implicit none
38  private
39  !-----------------------------------------------------------------------------
40  !
41  !++ Public procedure
42  !
43  public :: atmos_phy_lt_sato2019_setup
44  public :: atmos_phy_lt_sato2019_finalize
45  public :: atmos_phy_lt_sato2019_adjustment
46  public :: atmos_phy_lt_sato2019_select_dqcrg_from_lut
47 ! public :: ATMOS_PHY_LT_sato2019_mkinit
48 ! public :: ATMOS_PHY_LT_electric_field
49 ! public :: ATMOS_PHY_LT_neutralization_MG2001
50 ! public :: ATMOS_PHY_LT_neutralization_F2013
51 ! public :: ATMOS_PHY_LT_judge_absE
52  !-----------------------------------------------------------------------------
53  !
54  !++ Public parameters & variables
55  !
56  !-----------------------------------------------------------------------------
57  !
58  !++ Private procedure
59  !
60  !-----------------------------------------------------------------------------
61  !
62  !++ Private parameters & variables
63  !
64  !--- for Bi-CGSTAB to solve the poisson equation
65  integer, private :: ITMAX = 3000
66  real(RP), private :: epsilon = 0.1_rp ** (rp*2)
67  character(len=64),private :: NUTR_TYPE='F2013'
68  character(len=64),private,save :: NUTR_qhyd='Each_POLARITY2'
69  integer, private :: MAX_NUTR = 100
70  real(RP), private :: Eint = 150.0e+3_rp ! [V/m]
71  real(RP), private :: delEint = 10.0e+3_rp ! [V/m]
72  real(RP), private :: Estop = 15.0e+3_rp ! [V/m]
73  real(RP), private :: qrho_chan = 0.5_rp ! [nC/m3]
74  real(RP), private :: qrho_neut = 0.5_rp ! [nC/m3]
75  real(RP), private :: fp = 0.3_rp ! [-]
76  real(RP), private :: zcg = 0.0_rp ! lowest grid point of lightning path for MG2001 [m]
77  integer, private :: NUTR_ITMAX = 1000
78  integer, private :: FLAG_preprocessing = 2
82  integer, private :: KIJMAXG
83  character(len=H_LONG) :: ATMOS_PHY_LT_LUT_FILENAME !--- LUT file name
84  character(len=H_LONG) :: fname_lut_lt="LUT_TK1978_v.txt" !--- LUT file name
85  integer, save :: fid_lut_lt
86  real(RP), private :: R_neut = 2000.0_rp ! [m]
87  real(RP), private :: rho0 = 1.225_rp ! [kg/m3]
88  real(RP), save :: flg_eint_hgt
89  logical, private :: LT_DO_Lightning = .true.
90 
91  !--- for history output
92  real(RP), allocatable, private :: d_QCRG_TOT(:,:,:)
93  real(RP), allocatable, private :: LT_PATH_TOT(:,:,:,:)
94  real(RP), allocatable, private :: fls_int_p_tot(:,:,:)
95  real(RP), allocatable, private :: B_F2013_TOT(:,:)
96  real(RP), allocatable, private :: G_F2013(:,:)
97  real(RP), private :: C_F2013
98 
99  real(RP), allocatable, private :: A(:,:,:,:) !--- A : Laplasian (Coefficient matrix)
100  !$acc declare create(d_QCRG_TOT,LT_PATH_TOT,fls_int_p_tot,B_F2013_TOT,G_F2013,A)
101 
102  !---
103  integer, parameter, private :: nxlut_lt = 200, nylut_lt = 200
104  real(RP), private :: dq_chrg( nxlut_lt,nylut_lt ) !--- charge separation [fC]
105  real(RP), private :: grid_lut_t( nxlut_lt )
106  real(RP), private :: grid_lut_l( nylut_lt )
107 ! !$acc declare create(dq_chrg,grid_lut_t,grid_lut_l)
108  real(RP), private :: tcrglimit
109  logical, private :: Hgt_dependency_Eint = .false.
110 
111  !--- Index for local variable
112  integer, private, parameter :: I_lt_x = 1
113  integer, private, parameter :: I_lt_y = 2
114  integer, private, parameter :: I_lt_z = 3
115  integer, private, parameter :: I_lt_abs = 4
116 
117  !--- For history output
118  integer, private, parameter :: w_nmax = 11
119  integer, private, parameter :: I_Ex = 1
120  integer, private, parameter :: I_Ey = 2
121  integer, private, parameter :: I_Ez = 3
122  integer, private, parameter :: I_Eabs = 4
123  integer, private, parameter :: I_Epot = 5
124  integer, private, parameter :: I_Qneut = 6
125  integer, private, parameter :: I_LTpath = 7
126  integer, private, parameter :: I_PosFLASH = 8
127  integer, private, parameter :: I_NegFLASH = 9
128  integer, private, parameter :: I_FlashPoint = 10
129  integer, private, parameter :: I_FOD = 11
130  integer, private :: HIST_id(w_nmax)
131  character(len=H_SHORT), private :: w_name(w_nmax)
132  character(len=H_MID), private :: w_longname(w_nmax)
133  character(len=H_SHORT), private :: w_unit(w_nmax)
134  data w_name / 'Ex', &
135  'Ey', &
136  'Ez', &
137  'Eabs', &
138  'Epot', &
139  'Qneut', &
140  'LTpath', &
141  'PosFLASH', &
142  'NegFLASH', &
143  'FlashPoint', &
144  'FOD' /
145  data w_longname / &
146  'X component of Electrical Field', &
147  'Y component of Electrical Field', &
148  'Z component of Electrical Field', &
149  'Absolute value of Electrical Field', &
150  'Electric Potential', &
151  'Cumulative Neutralizated charge', &
152  'Cumulative Number of flash path', &
153  'Cumulative Number of Positive flash', &
154  'Cumulative Number of Negative flash', &
155  'Cumulative Number of Flash point', &
156  'Flash Origin Density' /
157  data w_unit / &
158  'kV/m', &
159  'kV/m', &
160  'kV/m', &
161  'kV/m', &
162  'V', &
163  'nC/m3/s', &
164  'num/grid/s', &
165  'num/grid/s', &
166  'num/grid/s', &
167  'num/grid/s', &
168  'num/grid/s' /
169  !-----------------------------------------------------------------------------
170 contains
171  !-----------------------------------------------------------------------------
173  subroutine atmos_phy_lt_sato2019_setup( KA, KS, KE, &
174  IA, IS, IE, &
175  JA, JS, JE, &
176  IMAXG, &
177  JMAXG, &
178  KMAX, &
179  MP_TYPE, &
180  CDX, CDY )
181  use scale_prc, only: &
182  prc_abort, &
183  prc_ismaster, &
184  prc_myrank
185  use scale_comm_cartesc, only: &
186  comm_bcast
187  use scale_atmos_grid_cartesc, only: &
188  cz => atmos_grid_cartesc_cz
189  use scale_const, only: &
190  t00 => const_tem00, &
191  pi => const_pi
192  use scale_file_history, only: &
193  file_history_reg
194  use scale_atmos_grid_cartesc, only: &
195  rcdz => atmos_grid_cartesc_rcdz, &
196  rcdx => atmos_grid_cartesc_rcdx, &
197  rcdy => atmos_grid_cartesc_rcdy, &
198  rfdz => atmos_grid_cartesc_rfdz, &
199  rfdx => atmos_grid_cartesc_rfdx, &
200  rfdy => atmos_grid_cartesc_rfdy
201  use scale_atmos_grid_cartesc_metric, only: &
202  mapf => atmos_grid_cartesc_metric_mapf, &
203  j13g => atmos_grid_cartesc_metric_j13g, &
204  j23g => atmos_grid_cartesc_metric_j23g, &
205  j33g => atmos_grid_cartesc_metric_j33g
206  use scale_atmos_grid_cartesc_index, only: &
207  i_xyz, &
208  i_xyw, &
209  i_uyw, &
210  i_xvw, &
211  i_uyz, &
212  i_xvz, &
213  i_uvz, &
214  i_xy, &
215  i_uy, &
216  i_xv, &
217  i_uv
218  implicit none
219  integer, intent(in) :: ka, ks, ke
220  integer, intent(in) :: ia, is, ie
221  integer, intent(in) :: ja, js, je
222  integer, intent(in) :: imaxg
223  integer, intent(in) :: jmaxg
224  integer, intent(in) :: kmax
225 
226  character(len=*), intent(in) :: mp_type
227  real(rp), intent(in) :: cdx(ia)
228  real(rp), intent(in) :: cdy(ja)
229 
230  integer :: n, myu, ip
231  integer :: ierr
232  integer :: i, j, k
233 
234  namelist / param_atmos_phy_lt_sato2019 / &
235  nutr_type, &
236  atmos_phy_lt_lut_filename, &
237  itmax, &
238  epsilon, &
239  eint, &
240  deleint, &
241  estop, &
242  qrho_chan, &
243  qrho_neut, &
244  fp, &
245  nutr_itmax, &
246  max_nutr, &
247  nutr_qhyd, &
248  zcg, &
249  r_neut, &
250  hgt_dependency_eint, &
251  lt_do_lightning, &
252  flag_preprocessing
253 
254  !--- read namelist
255  rewind(io_fid_conf)
256  read(io_fid_conf,nml=param_atmos_phy_lt_sato2019,iostat=ierr)
257 
258  if( ierr < 0 ) then !--- missing
259  log_info("ATMOS_PHY_LT_sato2019_setup",*) '*** Not found namelist. Default used.'
260  elseif( ierr > 0 ) then !--- fatal error
261  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_LT_SATO2019. Check!'
262  call prc_abort
263  endif
264  log_nml(param_atmos_phy_lt_sato2019)
265 
266  !--- If zcg is lower than lowest grid point, zcg set lowest grid point
267  if( zcg <= cz(ks) ) then
268  zcg = cz(ks)
269  endif
270 
271  if( r_neut <= minval( cdx,1 ) .or. r_neut <= minval( cdy,1 ) ) then
272  r_neut = 2.d0 * min( minval( cdx,1 ), minval( cdy,1 ) )
273  endif
274 
275  !$acc enter data create(dq_chrg,grid_lut_t,grid_lut_l)
276  if ( prc_ismaster ) then
277  fname_lut_lt = atmos_phy_lt_lut_filename
278  fid_lut_lt = io_get_available_fid()
279  !--- open parameter of cloud microphysics
280  open ( fid_lut_lt, file = fname_lut_lt, form = 'formatted', status = 'old', iostat=ierr )
281 
282  if ( ierr == 0 ) then
283  log_info("ATMOS_PHY_LT_sato2019_setup",'(2A)') 'Read LUT of TK78 table from ', trim(fname_lut_lt)
284  read( fid_lut_lt,* )
285  do n = 1, nylut_lt
286  do myu = 1, nxlut_lt
287  read( fid_lut_lt,* ) grid_lut_t( myu ), grid_lut_l( n ), dq_chrg( myu,n )
288  enddo
289  enddo
290  !$acc update device(dq_chrg,grid_lut_t,grid_lut_l)
291  !--- LUT file does not exist
292  else
293  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx LUT for LT is requied when ATMOS_PHY_LT_TYPE = SATO2019, stop!'
294  call prc_abort
295  endif
296 
297  if( nutr_type /= 'MG2001' .and. nutr_type /= 'F2013' ) then
298  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx NUTR_TYPE should be MG2001 or F2013 stop!'
299  call prc_abort
300  endif
301  endif
302 
303  call comm_bcast( nxlut_lt, nylut_lt, dq_chrg )
304  call comm_bcast( nxlut_lt, grid_lut_t )
305  call comm_bcast( nylut_lt, grid_lut_l )
306 
307 ! KIJMAXG = (IEG-ISG+1)*(JEG-JSG+1)*(KE-KS+1)
308  kijmaxg = imaxg*jmaxg*kmax
309 
310  ! for history output
311  allocate( d_qcrg_tot(ka,ia,ja) )
312  allocate( lt_path_tot(ka,ia,ja,3) )
313  allocate( fls_int_p_tot(ka,ia,ja) )
314  d_qcrg_tot(:,:,:) = 0.0_rp
315  lt_path_tot(:,:,:,:) = 0.0_rp
316  fls_int_p_tot(:,:,:) = 0.0_rp
317  !$acc update device(d_QCRG_TOT,LT_PATH_TOT,fls_int_p_tot)
318 
319  tcrglimit = -60.0_rp+t00
320 
321  if( nutr_type == 'F2013' ) then
322  log_info("ATMOS_PHY_LT_sato2019_setup",'(A,F15.7,A)') 'Radius of neutralization is ', r_neut, "[m]"
323  endif
324  allocate( b_f2013_tot(ia,ja) )
325  allocate( g_f2013(ia,ja) )
326  b_f2013_tot(:,:) = 0.0_rp
327  do i = 1, ia
328  do j = 1, ja
329  g_f2013(i,j) = cdx(i)*cdy(j)*1.0e-6_rp ! [m2] -> [km2]
330  enddo
331  enddo
332  !$acc update device(B_F2013_TOT,G_F2013)
333  c_f2013 = pi*r_neut*r_neut*1.0e-6_rp ! [m2] -> [km2]
334 
335  flg_eint_hgt = 0.0_rp
336  if( hgt_dependency_eint .and. nutr_type == 'MG2001') then
337  flg_eint_hgt = 1.0_rp
338  endif
339 
340  if( nutr_qhyd /= 'TOTAL' .and. nutr_qhyd /= 'Each_POLARITY' .and. nutr_qhyd /= 'Each_POLARITY2' ) then
341  log_error("ATMOS_PHY_LT_sato2019_setup",*) 'xxx NUTR_qhyd should be TOTAL, or Each_POLARITY, or Each_POLARITY2, stop!'
342  call prc_abort
343  endif
344 
345 ! if( ( NUTR_qhyd == 'Each_POLARITY' .or. NUTR_qhyd /= 'Each_POLARITY2' ) .and. MP_TYPE == 'SUZUKI10' ) then
346 ! LOG_ERROR("ATMOS_PHY_LT_sato2019_setup",*) 'NUTR_qhyd = Each_POLARITY, Each_POLARITY2 is not supported for MP_TYPE SUZUKI10'
347 ! call PRC_abort
348 ! endif
349 
350  do ip = 1, w_nmax
351  if( ip /= i_fod ) then
352  call file_history_reg( w_name(ip), w_longname(ip), w_unit(ip), & ! [IN]
353  hist_id(ip) ) ! [OUT]
354  elseif( ip == i_fod ) then
355  call file_history_reg( w_name(ip), w_longname(ip), w_unit(ip), & ! [IN]
356  hist_id(ip), dim_type='XY' ) ! [OUT]
357  endif
358  end do
359 
360  allocate( a(ka,15,ia,ja) )
361  !---- input vector A
362  !$omp parallel do
363  do j = js, je
364  do i = is, ie
365  do k = ks, ke
366  ! (k,i,j)
367  a(k,1,i,j) = &
368  - mapf(i,j ,1,i_xy)*mapf(i,j ,1,i_xy)*rcdx(i )*rfdx(i) &
369  - mapf(i,j ,1,i_xy)*mapf(i,j ,1,i_xy)*rcdx(i-1)*rfdx(i) &
370  + mapf(i,j ,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k ,i,j,2)*0.5_rp*rfdx(i)*rfdz(k) &
371  - 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) &
372  - 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) &
373  + 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) &
374  - j13g(k,i,j,i_xyz)*j13g(k ,i,j,i_xyw)*rcdz(k)*rfdz(k) &
375  - j13g(k,i,j,i_xyz)*j13g(k-1,i,j,i_xyw)*rcdz(k)*rfdz(k) &
376  - mapf(i,j ,2,i_xy)*mapf(i,j ,2,i_xy)*rcdy(j )*rfdy(j) &
377  - mapf(i,j ,2,i_xy)*mapf(i,j ,2,i_xy)*rcdy(j-1)*rfdy(j) &
378  + mapf(i,j ,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k ,i,j,2)*0.5_rp*rfdy(j)*rfdz(k) &
379  - 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) &
380  - 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) &
381  + 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) &
382  - j23g(k,i,j,i_xyz)*j23g(k ,i,j,i_xyw)*rcdz(k)*rfdz(k) &
383  - j23g(k,i,j,i_xyz)*j23g(k-1,i,j,i_xyw)*rcdz(k)*rfdz(k) &
384  - j33g*j33g*rfdz(k)*rcdz(k) &
385  - j33g*j33g*rfdz(k)*rcdz(k-1)
386 
387  ! (k-1,i,j)
388  a(k,2,i,j) = &
389  - 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) &
390  + 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) &
391  + j13g(k,i,j,i_xyz)*j13g(k-1,i,j,i_xyw)*rfdz(k)*rcdz(k-1) &
392  - 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) &
393  + 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) &
394  + j23g(k,i,j,i_xyz)*j23g(k-1,i,j,i_xyw)*rfdz(k)*rcdz(k-1) &
395  + j33g*j33g*rfdz(k)*rcdz(k-1)
396 
397  ! (k+1,i,j)
398  a(k,3,i,j) = &
399  mapf(i,j,1,i_xy)*j13g(k,i ,j,i_uyz)*f2h(k,i,j,1)*0.50_rp*rfdx(i)*rfdz(k) &
400  - 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) &
401  + j13g(k,i,j,i_xyz)*j13g(k,i,j,i_xyw)*rfdz(k)*rcdz(k) &
402  + mapf(i,j,2,i_xy)*j23g(k,i,j ,i_xvz)*f2h(k,i,j,1)*0.50_rp*rfdy(j)*rfdz(k) &
403  - 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) &
404  + j23g(k,i,j,i_xyz)*j23g(k,i,j,i_xyw)*rfdz(k)*rcdz(k) &
405  + j33g*j33g*rfdz(k)*rcdz(k)
406 
407  ! (k,i-1,j)
408  a(k,4,i,j) = &
409  mapf(i,j,1,i_xy)*mapf(i-1,j,1,i_xy)*rfdx(i)*rcdx(i-1) &
410  - 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) &
411  + 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) &
412  - 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) &
413  + 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)
414 
415  ! (k,i+1,j)
416  a(k,5,i,j) = &
417  mapf(i,j,1,i_xy)*mapf(i+1,j,1,i_xy)*rfdx(i)*rcdx(i) &
418  + 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) &
419  - 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) &
420  + 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) &
421  - 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)
422 
423  ! (k,i,j-1)
424  a(k,6,i,j) = &
425  mapf(i,j,2,i_xy)*mapf(i,j-1,2,i_xy)*rfdy(j)*rcdy(j-1) &
426  - 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) &
427  + 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) &
428  - 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) &
429  + 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)
430 
431  ! (k,i,j+1)
432  a(k,7,i,j) = &
433  mapf(i,j,2,i_xy)*mapf(i,j+1,2,i_xy)*rfdy(j)*rcdy(j) &
434  + 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) &
435  - 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) &
436  + 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) &
437  - 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)
438 
439  ! (k-1,i-1,j)
440  a(k,8,i,j) = &
441  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) &
442  + 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)
443 
444  ! (k-1,i+1,j)
445  a(k,9,i,j) = &
446  - 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) &
447  - 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)
448 
449  ! (k-1,i,j-1)
450  a(k,10,i,j) = &
451  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) &
452  + 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)
453 
454  ! (k-1,i,j+1)
455  a(k,11,i,j) = &
456  - 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) &
457  - 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)
458 
459  ! (k+1,i-1,j)
460  a(k,12,i,j) = &
461  - 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) &
462  - 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)
463 
464  ! (k+1,i+1,j)
465  a(k,13,i,j) = &
466  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) &
467  + 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)
468 
469  ! (k+1,i,j-1)
470  a(k,14,i,j) = &
471  - 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) &
472  - 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)
473 
474  ! (k+1,i,j+1)
475  a(k,15,i,j) = &
476  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) &
477  + 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)
478 
479  enddo
480  enddo
481  enddo
482  !$acc update device(A)
483 
484  return
485  end subroutine atmos_phy_lt_sato2019_setup
486 
487  !-----------------------------------------------------------------------------
489  subroutine atmos_phy_lt_sato2019_finalize
490 
491  !$acc exit data delete(dq_chrg,grid_lut_t,grid_lut_l)
492 
493  deallocate( d_qcrg_tot )
494  deallocate( lt_path_tot )
495  deallocate( fls_int_p_tot )
496 
497  deallocate( b_f2013_tot )
498  deallocate( g_f2013 )
499 
500  return
501  end subroutine atmos_phy_lt_sato2019_finalize
502 
503  !-----------------------------------------------------------------------------
505  subroutine atmos_phy_lt_sato2019_adjustment( &
506  KA, KS, KE, &
507  IA, IS, IE, &
508  JA, JS, JE, &
509  KIJMAX, &
510  IMAX, &
511  JMAX, &
512  QA_LT, &
513  DENS, &
514  RHOT, &
515  QHYD, &
516  Sarea, &
517  dt_LT, &
518  QTRC, &
519  Epot )
520  use scale_const, only: &
521  small => const_eps
522  use scale_prc, only: &
523  prc_ismaster, &
524  prc_abort
525  use scale_comm_cartesc, only: &
526  comm_datatype, &
527  comm_world, &
528  comm_wait, &
529  comm_vars8
530  use scale_file_history, only: &
531  file_history_query, &
532  file_history_put
533  implicit none
534 
535  integer, intent(in) :: ka, ks, ke
536  integer, intent(in) :: ia, is, ie
537  integer, intent(in) :: ja, js, je
538  integer, intent(in) :: kijmax
539  integer, intent(in) :: imax
540  integer, intent(in) :: jmax
541 ! character(len=H_SHORT)
542  integer, intent(in) :: qa_lt
543  real(rp), intent(in) :: dens(ka,ia,ja)
544  real(rp), intent(in) :: rhot(ka,ia,ja)
545  real(rp), intent(in) :: qhyd(ka,ia,ja)
546  real(rp), intent(in) :: sarea(ka,ia,ja,qa_lt) !--- Surface area and that of each catergory [m2]
547  real(dp), intent(in) :: dt_lt
548  real(rp), intent(inout) :: qtrc(ka,ia,ja,qa_lt)
549  real(rp), intent(inout) :: epot(ka,ia,ja) !--- Electrical potential at previous time step [V]
550 
551  real(rp) :: qhyd_mass(ka,ia,ja) !--- Mass of total hydrometeor [kg/m3]
552  real(rp) :: qcrg(ka,ia,ja) !--- Total charge density [nC/m3]
553  real(rp) :: efield(ka,ia,ja,i_lt_abs) !--- Electrical field (1-3 ->x,y,z, 4->abs. )
554  real(rp) :: num_end(ka,ia,ja,3) !--- Number of each flash type (1->negative, 2->ground, 3->positive)
555  real(rp) :: d_qcrg(ka,ia,ja) !--- Change of charge by charge neutralization [fC/m3]
556  real(rp) :: lt_path(ka,ia,ja) !--- Lightning path (0-> no path, 1-> path)
557  real(rp) :: fls_int_p(ka,ia,ja)
558  real(rp) :: total_sarea1, total_sarea2 !--- Sum of surface area and that of each category [m2]
559  real(rp) :: r_totalsarea1, r_totalsarea2
560  real(rp) :: neg_crg, pos_crg
561  real(rp) :: frac
562  real(rp) :: dqneut(ka,ia,ja,qa_lt)
563  real(rp) :: dqneut_real(ka,ia,ja,qa_lt)
564  real(rp) :: dqneut_real_tot(ka,ia,ja)
565  logical :: flg_chrged(qa_lt)
566  real(rp) :: emax, emax_old
567  logical :: output_step
568  integer :: flg_lt_neut
569  integer :: i, j, k, m, n, countbin, ip
570  real(rp) :: sw, zerosw, positive, negative
571  integer :: count_neut
572 
573  logical :: hist_sw(w_nmax)
574  real(rp) :: w3d(ka,ia,ja)
575 
576  real(rp) :: diff_qcrg(0:1), lack(0:1), sum_crg(0:1)
577  real(rp) :: crg_rate(qa_lt), qcrg_before(qa_lt)
578  real(rp) :: b_f2013(ia,ja)
579  integer :: int_sw
580  real(rp) :: iprod, buf, tmp_qcrg, tmp_dqcrg
581  integer :: ierror
582 
583  !$acc data &
584  !$acc copyin(DENS,RHOT,QHYD,Sarea) &
585  !$acc copy(QTRC,Epot) &
586  !$acc create(QHYD_mass,QCRG,Efield,NUM_end,d_QCRG,LT_PATH,fls_int_p, &
587  !$acc dqneut,dqneut_real,dqneut_real_tot,flg_chrged,HIST_sw,w3d, &
588  !$acc diff_qcrg,lack,sum_crg,crg_rate,qcrg_before,B_F2013)
589 
590  !$acc kernels
591  num_end(:,:,:,:) = 0.0_rp
592  b_f2013(:,:) = 0.0_rp
593  !$acc end kernels
594 
595  !$omp parallel do private(tmp_qcrg)
596  !$acc kernels
597  !$acc loop collapse(2)
598  do j = js, je
599  do i = is, ie
600 
601  ! calc total charge density
602  do k = ks, ke
603  tmp_qcrg = 0.0_rp
604  !$acc loop seq
605  do n = 1, qa_lt
606  tmp_qcrg = tmp_qcrg + qtrc(k,i,j,n)
607  end do
608  qcrg(k,i,j) = tmp_qcrg * dens(k,i,j) * 1.e-6_rp ![fC/kg] -> [nc/m3]
609  enddo
610 
611  do k = ks, ke
612  qhyd_mass(k,i,j) = qhyd(k,i,j) * dens(k,i,j) ![kg/kg] -> [kg/m3]
613  enddo
614 
615  enddo
616  enddo
617  !$acc end kernels
618 
619  iprod = 0.0_rp
620  !$omp parallel do reduction(+:iprod) private(i,j,k)
621  !$acc kernels
622  !$acc loop collapse(3) reduction(+:iprod)
623  do j = js, je
624  do i = is, ie
625  do k = ks, ke
626  iprod = iprod + abs( qcrg(k,i,j) )
627  enddo
628  enddo
629  enddo
630  !$acc end kernels
631 
632  call mpi_allreduce(iprod, buf, 1, comm_datatype, mpi_sum, comm_world, ierror)
633 
634  if( buf <= small ) then
635 
636  !$omp parallel do
637  !$acc kernels
638  !$acc loop collapse(3)
639  do j = js, je
640  do i = is, ie
641  do k = ks, ke
642  epot(k,i,j) = 0.0_rp
643  efield(k,i,j,i_lt_x) = 0.0_rp
644  efield(k,i,j,i_lt_y) = 0.0_rp
645  efield(k,i,j,i_lt_z) = 0.0_rp
646  end do
647  end do
648  end do
649  !$acc end kernels
650 
651  else
652 
653  !--- Calculate E field
654  call atmos_phy_lt_electric_field( ka, ks, ke, & ! [IN]
655  ia, is, ie, & ! [IN]
656  ja, js, je, & ! [IN]
657  qcrg(:,:,:), & ! [IN]
658  dens(:,:,:), & ! [IN]
659  rhot(:,:,:), & ! [IN]
660  epot(:,:,:), & ! [INOUT]
661  efield(:,:,:,i_lt_x:i_lt_z) ) ! [INOUT]
662 
663  endif
664 
665  !$omp parallel do
666  !$acc kernels
667  !$acc loop collapse(3)
668  do j = js, je
669  do i = is, ie
670  do k = ks, ke
671  efield(k,i,j,i_lt_abs) = sqrt( efield(k,i,j,i_lt_x)*efield(k,i,j,i_lt_x) &
672  + efield(k,i,j,i_lt_y)*efield(k,i,j,i_lt_y) &
673  + efield(k,i,j,i_lt_z)*efield(k,i,j,i_lt_z) )
674 
675  lt_path(k,i,j) = 0.0_rp
676 
677  enddo
678  enddo
679  enddo
680  !$acc end kernels
681 
682  if ( lt_do_lightning ) then
683 
684  !$acc kernels
685  d_qcrg_tot(:,:,:) = 0.0_rp
686  fls_int_p_tot(:,:,:) = 0.0_rp
687  lt_path_tot(:,:,:,:) = 0.0_rp
688  b_f2013_tot(:,:) = 0.0_rp
689  !$acc end kernels
690 
691  call atmos_phy_lt_judge_abse( ka, ks, ke, & ! [IN]
692  ia, is, ie, & ! [IN]
693  ja, js, je, & ! [IN]
694  dens(:,:,:), & ! [IN]
695  efield(:,:,:,i_lt_abs), & ! [IN]
696  emax, & ! [OUT]
697  flg_lt_neut ) ! [OUT]
698 
699  count_neut = 0
700  do while( flg_lt_neut > 0 )
701 
702  emax_old = emax
703 
704  call comm_vars8( efield(:,:,:,i_lt_x), 1 )
705  call comm_vars8( efield(:,:,:,i_lt_y), 2 )
706  call comm_vars8( efield(:,:,:,i_lt_z), 3 )
707  call comm_vars8( efield(:,:,:,i_lt_abs), 4 )
708  call comm_vars8( qhyd_mass(:,:,:), 5 )
709  call comm_vars8( qcrg(:,:,:), 6 )
710  call comm_vars8( epot(:,:,:), 7 )
711  call comm_wait ( efield(:,:,:,i_lt_x), 1 )
712  call comm_wait ( efield(:,:,:,i_lt_y), 2 )
713  call comm_wait ( efield(:,:,:,i_lt_z), 3 )
714  call comm_wait ( efield(:,:,:,i_lt_abs), 4 )
715  call comm_wait ( qhyd_mass(:,:,:), 5 )
716  call comm_wait ( qcrg(:,:,:), 6 )
717  call comm_wait ( epot(:,:,:), 7 )
718 
719  !--- Calculate lightning path and charge neutralization
720  if( nutr_type == 'MG2001' ) then
721  !$acc update host(Efield,Epot,QCRG,QHYD_mass,NUM_end,LT_PATH)
722  call atmos_phy_lt_neutralization_mg2001( &
723  ka, ks, ke, & ! [IN]
724  ia, is, ie, & ! [IN]
725  ja, js, je, & ! [IN]
726  kijmax, imax, jmax, & ! [IN]
727  efield(:,:,:,:), & ! [IN]
728  epot(:,:,:), & ! [IN]
729  dens(:,:,:), & ! [IN]
730  qcrg(:,:,:), & ! [IN]
731  qhyd_mass(:,:,:), & ! [IN]
732  num_end(:,:,:,:), & ! [INOUT]
733  lt_path(:,:,:), & ! [INOUT]
734  fls_int_p(:,:,:), & ! [OUT]
735  d_qcrg(:,:,:) ) ! [OUT]
736  !$acc update device(NUM_end,LT_PATH,fls_int_p,d_QCRG)
737  elseif( nutr_type == 'F2013' ) then
738  call atmos_phy_lt_neutralization_f2013( &
739  ka, ks, ke, & ! [IN]
740  ia, is, ie, & ! [IN]
741  ja, js, je, & ! [IN]
742  kijmax, & ! [IN]
743  efield(:,:,:,:), & ! [IN]
744  epot(:,:,:), & ! [IN]
745  dens(:,:,:), & ! [IN]
746  qcrg(:,:,:), & ! [IN]
747  qhyd_mass(:,:,:), & ! [IN]
748  num_end(:,:,:,:), & ! [INOUT]
749  lt_path(:,:,:), & ! [INOUT]
750  fls_int_p(:,:,:), & ! [OUT]
751  d_qcrg(:,:,:), & ! [OUT]
752  b_f2013(:,:) ) ! [OUT]
753  endif
754 
755  call comm_vars8( lt_path(:,:,:),1 )
756  call comm_vars8( d_qcrg(:,:,:),2 )
757  call comm_wait ( lt_path(:,:,:),1 )
758  call comm_wait ( d_qcrg(:,:,:),2 )
759 
760  !$acc kernels
761  dqneut(:,:,:,:) = 0.0_rp
762  dqneut_real(:,:,:,:) = 0.0_rp
763  !$acc end kernels
764 
765  !-- Calculate neutralization of each hydrometeor or each category
766  select case( nutr_qhyd )
767  case ( 'TOTAL' )
768 
769  !$omp parallel do &
770  !$omp private(Total_Sarea1,Total_Sarea2,r_totalSarea1,r_totalSarea2,zerosw)
771  !$acc kernels
772  !$acc loop collapse(3)
773  do j = js, je
774  do i = is, ie
775  do k = ks, ke
776  if( abs( d_qcrg(k,i,j) ) > 0.0_rp ) then
777  total_sarea1 = 0.0_rp
778  !$acc loop seq
779  do n = 1, qa_lt
780  total_sarea1 = total_sarea1 + sarea(k,i,j,n)
781  enddo
782  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea1-small )
783  r_totalsarea1 = 1.0_rp / ( total_sarea1 + zerosw ) * ( 1.0_rp - zerosw )
784  !$acc loop seq
785  do n = 1, qa_lt
786  dqneut(k,i,j,n) = d_qcrg(k,i,j)*1.0e+6_rp &
787  * sarea(k,i,j,n) * r_totalsarea1 / dens(k,i,j)
788  qtrc(k,i,j,n) = qtrc(k,i,j,n) + dqneut(k,i,j,n)
789  dqneut_real(k,i,j,n) = dqneut(k,i,j,n)
790  enddo
791  endif
792  enddo
793  enddo
794  enddo
795  !$acc end kernels
796 
797  case ( 'Each_POLARITY', 'Each_POLARITY2' )
798 
799  !$omp parallel do &
800  !$omp private(Total_Sarea1,Total_Sarea2,r_totalSarea1,r_totalSarea2,flg_chrged,pos_crg,neg_crg,frac,lack, &
801  !$omp diff_qcrg,crg_rate,sum_crg,qcrg_before, &
802  !$omp positive,negative,zerosw,sw,int_sw)
803  !$acc parallel
804  !$acc loop collapse(2) gang
805  do j = js, je
806  do i = is, ie
807  !$acc loop vector private(lack,flg_chrged,diff_qcrg,crg_rate,sum_crg)
808  do k = ks, ke
809  lack(:) = 0.0_rp
810  if( abs( d_qcrg(k,i,j) ) > 0.0_rp ) then
811 
812  !--- flg whether the charged or not (0.0-> not charged, 1.0->charged)
813  !$acc loop seq
814  do n = 1, qa_lt
815  flg_chrged(n) = abs(qtrc(k,i,j,n)) >= small
816  enddo
817 
818  total_sarea1 = 0.0_rp
819  total_sarea2 = 0.0_rp
820  pos_crg = 0.0_rp
821  neg_crg = 0.0_rp
822  !$acc loop seq
823  do n = 1, qa_lt
824  if ( flg_chrged(n) ) then
825  positive = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
826  negative = 1.0_rp - positive
827  !--- total of positive charge
828  pos_crg = pos_crg + qtrc(k,i,j,n) * positive
829  !--- total of negative charge
830  neg_crg = neg_crg + qtrc(k,i,j,n) * negative
831  !--- Sarea of positively charged hydrometeor
832  total_sarea1 = total_sarea1 + sarea(k,i,j,n) * positive
833  !--- Sarea of negatively charged hydrometeor
834  total_sarea2 = total_sarea2 + sarea(k,i,j,n) * negative
835  end if
836  end do
837 
838  zerosw = 0.5_rp - sign( 0.5_rp, abs( qcrg(k,i,j) ) - small )
839  frac = d_qcrg(k,i,j) / ( qcrg(k,i,j) + zerosw ) * ( 1.0_rp - zerosw )
840  pos_crg = frac * pos_crg
841  neg_crg = frac * neg_crg
842 
843  !--- remove 0 surface area ( no crg for each porality )
844  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea1 - small )
845  r_totalsarea1 = 1.0_rp / ( total_sarea1 + zerosw ) * ( 1.0_rp - zerosw )
846  zerosw = 0.5_rp - sign( 0.5_rp, total_sarea2 - small )
847  r_totalsarea2 = 1.0_rp / ( total_sarea2 + zerosw ) * ( 1.0_rp - zerosw )
848 
849  diff_qcrg(:) = 0.0_rp
850  !$acc loop seq
851  do n = 1, qa_lt
852  crg_rate(n) = 0.0_rp
853  end do
854  sum_crg(:) = 0.0_rp
855  !$acc loop seq
856  do n = 1, qa_lt
857  if ( flg_chrged(n) ) then
858  sw = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
859  dqneut(k,i,j,n) = &
860  + pos_crg * sarea(k,i,j,n) * r_totalsarea1 &
861  * sw & ! sw = 1 positive
862  + neg_crg * sarea(k,i,j,n) * r_totalsarea2 &
863  * ( 1.0_rp - sw ) ! sw = 0 negative
864  qcrg_before(n) = qtrc(k,i,j,n)
865 
866  if( sw == 1.0_rp ) then
867  qtrc(k,i,j,n) = max( qtrc(k,i,j,n) + dqneut(k,i,j,n), 0.0_rp ) !--- limiter
868  elseif( sw == 0.0_rp ) then
869  qtrc(k,i,j,n) = min( qtrc(k,i,j,n) + dqneut(k,i,j,n), 0.0_rp ) !--- limiter
870  endif
871 
872  int_sw = int( sw ) ! 0-> negative, 1-> positive
873  diff_qcrg(int_sw) = diff_qcrg(int_sw) &
874  + ( qtrc(k,i,j,n) - qcrg_before(n) )
875  sum_crg(int_sw) = sum_crg(int_sw) + qtrc(k,i,j,n)
876  end if
877  enddo
878 
879  if( nutr_qhyd == 'Each_POLARITY2' ) then ! Adjust
880  lack(0) = neg_crg - diff_qcrg(0) ! negative (Should be Positive)
881  lack(1) = pos_crg - diff_qcrg(1) ! positive (Should be Negative)
882 #ifdef DEBUG
883  if( lack(0) > 0.0_rp .and. abs(lack(0)/diff_qcrg(0)) > 1.0e-10_rp ) then
884  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,4E15.7)') &
885  "Large negative for lack(0) ", lack(0)/diff_qcrg(0), neg_crg, lack(0), diff_qcrg(0)
886  endif
887  if( lack(1) < 0.0_rp .and. abs(lack(1)/diff_qcrg(1)) > 1.0e-10_rp ) then
888  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,4E15.7)') &
889  "Large positive for lack(1) ", lack(1)/diff_qcrg(1), pos_crg, lack(1), diff_qcrg(1)
890  endif
891 #endif
892  lack(0) = max( lack(0), 0.0_rp ) ! negative (Should be Positive)
893  lack(1) = min( lack(1), 0.0_rp ) ! positive (Should be Negative)
894  !$acc loop seq
895  do n = 1, qa_lt
896  if ( flg_chrged(n) ) then
897  sw = 0.5_rp + sign( 0.5_rp, qtrc(k,i,j,n) )
898  int_sw = int( sw ) ! 0-> negative, 1-> positive
899  if( sum_crg(int_sw) /= 0.0_rp ) then
900  crg_rate(n) = qtrc(k,i,j,n)/sum_crg(int_sw)
901  else
902  crg_rate(n) = 0.0_rp
903  endif
904  qtrc(k,i,j,n) = qtrc(k,i,j,n) + crg_rate(n) * lack(int_sw)
905  endif
906  enddo
907  endif
908 
909  !$acc loop seq
910  do n = 1, qa_lt
911  if ( flg_chrged(n) ) then
912  dqneut_real(k,i,j,n) = qtrc(k,i,j,n) - qcrg_before(n)
913  endif
914  enddo
915 
916  endif
917 
918  enddo
919  enddo
920  enddo
921  !$acc end parallel
922 
923  end select
924 
925  !$acc kernels
926  !$acc loop collapse(3)
927  do j = js, je
928  do i = is, ie
929 
930  ! calc total charge density
931  do k = ks, ke
932  tmp_qcrg = 0.0_rp
933  tmp_dqcrg = 0.0_rp
934  !$acc loop seq
935  do n = 1, qa_lt
936  tmp_qcrg = tmp_qcrg + qtrc(k,i,j,n)
937  tmp_dqcrg = tmp_dqcrg + dqneut_real(k,i,j,n)
938  end do
939  qcrg(k,i,j) = tmp_qcrg * dens(k,i,j) * 1.e-6_rp ![fC/kg] -> [nc/m3]
940  dqneut_real_tot(k,i,j) = tmp_dqcrg
941  enddo
942 
943  enddo
944  enddo
945  !$acc end kernels
946 
947  iprod = 0.0_rp
948  !$omp parallel do reduction(+:iprod) private(i,j,k)
949  !$acc kernels
950  !$acc loop collapse(3) reduction(+:iprod)
951  do j = js, je
952  do i = is, ie
953  do k = ks, ke
954  iprod = iprod + abs( qcrg(k,i,j) )
955  enddo
956  enddo
957  enddo
958  !$acc end kernels
959 
960  call mpi_allreduce(iprod, buf, 1, comm_datatype, mpi_sum, comm_world, ierror)
961 
962  if( buf <= small ) then
963 
964  !$omp parallel do
965  !$acc kernels
966  !$acc loop collapse(3)
967  do j = js, je
968  do i = is, ie
969  do k = ks, ke
970  epot(k,i,j) = 0.0_rp
971  efield(k,i,j,i_lt_x) = 0.0_rp
972  efield(k,i,j,i_lt_y) = 0.0_rp
973  efield(k,i,j,i_lt_z) = 0.0_rp
974  end do
975  end do
976  end do
977  !$acc end kernels
978 
979  else
980 
981  !--- Calculate E field
982  call atmos_phy_lt_electric_field( ka, ks, ke, & ! [IN]
983  ia, is, ie, & ! [IN]
984  ja, js, je, & ! [IN]
985  qcrg(:,:,:), & ! [IN]
986  dens(:,:,:), & ! [IN]
987  rhot(:,:,:), & ! [IN]
988  epot(:,:,:), & ! [INOUT]
989  efield(:,:,:,i_lt_x:i_lt_z) ) ! [INOUT]
990 
991  endif
992 
993  !--- Add Total number of path
994  !$omp parallel do
995  !$acc kernels
996  !$acc loop collapse(3)
997  do j = js, je
998  do i = is, ie
999  do k = ks, ke
1000  efield(k,i,j,i_lt_abs) = sqrt( efield(k,i,j,i_lt_x)*efield(k,i,j,i_lt_x) &
1001  + efield(k,i,j,i_lt_y)*efield(k,i,j,i_lt_y) &
1002  + efield(k,i,j,i_lt_z)*efield(k,i,j,i_lt_z) )
1003  end do
1004  end do
1005  end do
1006  !$acc end kernels
1007 
1008  !--- Add Total number of charge neutralization and flash point
1009  if ( hist_id(i_qneut) > 0 ) then
1010  !$omp parallel do
1011  !$acc kernels
1012  do j = js, je
1013  do i = is, ie
1014  do k = ks, ke
1015  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]
1016  end do
1017  end do
1018  end do
1019  !$acc end kernels
1020  end if
1021  if ( hist_id(i_flashpoint) > 0 ) then
1022  !$omp parallel do
1023  !$acc kernels
1024  do j = js, je
1025  do i = is, ie
1026  do k = ks, ke
1027  fls_int_p_tot(k,i,j) = fls_int_p_tot(k,i,j) + fls_int_p(k,i,j)
1028  end do
1029  end do
1030  end do
1031  !$acc end kernels
1032  end if
1033 
1034  if ( hist_id(i_posflash) > 0 ) then
1035  !$omp parallel do
1036  !$acc kernels
1037  do j = js, je
1038  do i = is, ie
1039  do k = ks, ke
1040  lt_path_tot(k,i,j,1) = lt_path_tot(k,i,j,1) &
1041  + 0.5_rp + sign( 0.5_rp,-dqneut_real_tot(k,i,j)-small )
1042  end do
1043  end do
1044  end do
1045  !$acc end kernels
1046  end if
1047  if ( hist_id(i_negflash) > 0 ) then
1048  !$omp parallel do
1049  !$acc kernels
1050  do j = js, je
1051  do i = is, ie
1052  do k = ks, ke
1053  lt_path_tot(k,i,j,2) = lt_path_tot(k,i,j,2) &
1054  + 0.5_rp + sign( 0.5_rp, dqneut_real_tot(k,i,j)-small )
1055  end do
1056  end do
1057  end do
1058  !$acc end kernels
1059  end if
1060  if ( hist_id(i_ltpath) > 0 ) then
1061  !$omp parallel do
1062  !$acc kernels
1063  do j = js, je
1064  do i = is, ie
1065  do k = ks, ke
1066  lt_path_tot(k,i,j,3) = lt_path_tot(k,i,j,3) + lt_path(k,i,j)
1067  end do
1068  end do
1069  end do
1070  !$acc end kernels
1071  end if
1072  if ( hist_id(i_fod) > 0 ) then
1073  !$acc kernels
1074  do j = js, je
1075  do i = is, ie
1076  b_f2013_tot(i,j) = b_f2013_tot(i,j) + g_f2013(i,j)/c_f2013*b_f2013(i,j)
1077  enddo
1078  enddo
1079  !$acc end kernels
1080  endif
1081 
1082  call atmos_phy_lt_judge_abse( ka, ks, ke, & ! [IN]
1083  ia, is, ie, & ! [IN]
1084  ja, js, je, & ! [IN]
1085  dens(:,:,:), & ! [IN]
1086  efield(:,:,:,i_lt_abs), & ! [IN]
1087  emax, & ! [OUT]
1088  flg_lt_neut ) ! [OUT]
1089 
1090  count_neut = count_neut + 1
1091 #ifdef DEBUG
1092  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,F15.7,A,F15.7,1X,I0)') &
1093  'CHECK', &
1094  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, count_neut
1095 #endif
1096 
1097  if( flg_lt_neut == 1 .and. emax == emax_old ) then
1098  flg_lt_neut = 0
1099  if( prc_ismaster ) then
1100  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(A,2(E15.7,A),1X,I0)') &
1101  'Eabs value after neutralization is same as previous value, Finish', &
1102  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, count_neut
1103  endif
1104  elseif( flg_lt_neut == 1 .and. count_neut == max_nutr ) then
1105  flg_lt_neut = 0
1106  if( prc_ismaster ) then
1107  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
1108  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
1109  ' [kV/m], reach maximum neutralization count, Finish', count_neut
1110  endif
1111  elseif( flg_lt_neut == 1 .and. emax > emax_old ) then
1112  flg_lt_neut = 0
1113  if( prc_ismaster ) then
1114  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
1115  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
1116  ' [kV/m] larger than previous one by neutralization, back to previous step, Finish', &
1117  count_neut
1118  endif
1119  !---- Back to Charge density as previous step
1120  !$acc kernels
1121  !$acc loop collapse(3)
1122  do j = js, je
1123  do i = is, ie
1124  do k = ks, ke
1125  !$acc loop seq
1126  do n = 1, qa_lt
1127  qtrc(k,i,j,n) = qtrc(k,i,j,n) - dqneut_real(k,i,j,n)
1128  enddo
1129  d_qcrg_tot(k,i,j) = d_qcrg_tot(k,i,j) - dqneut_real_tot(k,i,j)*1.e-6_rp
1130  d_qcrg(k,i,j) = 0.0_rp
1131  enddo
1132  enddo
1133  enddo
1134  !$acc end kernels
1135  elseif( flg_lt_neut == 2 ) then
1136  flg_lt_neut = 0
1137  if( prc_ismaster ) then
1138  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(E15.7,A),1X,I0)') &
1139  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
1140  '[kV/m] After neutralization, Finish' , count_neut
1141  endif
1142  elseif( flg_lt_neut == 0 ) then
1143  if( prc_ismaster ) then
1144  log_info("ATMOS_PHY_LT_sato2019_adjustment",'(2(F15.7,A),1X,I0)') &
1145  emax_old*1.e-3_rp, ' [kV/m] -> ', emax*1.e-3_rp, &
1146  '[kV/m] After neutralization, Finish', count_neut
1147  endif
1148  endif
1149 
1150  enddo
1151 
1152  else
1153 
1154  !$omp parallel do private(i,j,k)
1155  !$acc kernels
1156  !$acc loop collapse(3)
1157  do j = js, je
1158  do i = is, ie
1159  do k = ks, ke
1160  d_qcrg(k,i,j) = 0.0_rp
1161  end do
1162  end do
1163  end do
1164  !$acc end kernels
1165 
1166  endif
1167 
1168  !--- For history output
1169  do ip = 1, w_nmax
1170  call file_history_query( hist_id(ip), hist_sw(ip) )
1171  end do
1172 
1173  if ( hist_sw(i_ex ) ) then
1174  !$omp parallel do private(i,j,k)
1175  !$acc kernels
1176  !$acc loop collapse(3)
1177  do j = js, je
1178  do i = is, ie
1179  do k = ks, ke
1180  w3d(k,i,j) = efield(k,i,j,i_lt_x )*1.0e-3_rp ![kV/m]
1181  enddo
1182  enddo
1183  enddo
1184  !$acc end kernels
1185  call file_history_put( hist_id(i_ex ), w3d(:,:,:) )
1186  endif
1187  if ( hist_sw(i_ey ) ) then
1188  !$omp parallel do private(i,j,k)
1189  !$acc kernels
1190  !$acc loop collapse(3)
1191  do j = js, je
1192  do i = is, ie
1193  do k = ks, ke
1194  w3d(k,i,j) = efield(k,i,j,i_lt_y )*1.0e-3_rp ![kV/m]
1195  enddo
1196  enddo
1197  enddo
1198  !$acc end kernels
1199  call file_history_put( hist_id(i_ey ), w3d(:,:,:) )
1200  endif
1201  if ( hist_sw(i_ez ) ) then
1202  !$omp parallel do private(i,j,k)
1203  !$acc kernels
1204  !$acc loop collapse(3)
1205  do j = js, je
1206  do i = is, ie
1207  do k = ks, ke
1208  w3d(k,i,j) = efield(k,i,j,i_lt_z )*1.0e-3_rp ![kV/m]
1209  enddo
1210  enddo
1211  enddo
1212  !$acc end kernels
1213  call file_history_put( hist_id(i_ez ), w3d(:,:,:) )
1214  endif
1215  if ( hist_sw(i_eabs) ) then
1216  !$omp parallel do private(i,j,k)
1217  !$acc kernels
1218  !$acc loop collapse(3)
1219  do j = js, je
1220  do i = is, ie
1221  do k = ks, ke
1222  w3d(k,i,j) = efield(k,i,j,i_lt_abs)*1.0e-3_rp ![kV/m]
1223  enddo
1224  enddo
1225  enddo
1226  !$acc end kernels
1227  call file_history_put( hist_id(i_eabs), w3d(:,:,:) )
1228  endif
1229  if ( hist_sw(i_epot) ) then
1230  call file_history_put( hist_id(i_epot), epot(:,:,:) )
1231  end if
1232  if ( hist_sw(i_qneut) ) then
1233  !$omp parallel do private(i,j,k)
1234  !$acc kernels
1235  do j = js, je
1236  do i = is, ie
1237  do k = ks, ke
1238  w3d(k,i,j) = d_qcrg_tot(k,i,j)/dt_lt ![nC/m3/s]
1239  enddo
1240  enddo
1241  enddo
1242  !$acc end kernels
1243  call file_history_put( hist_id(i_qneut), w3d(:,:,:) )
1244  endif
1245  if ( hist_sw(i_ltpath) ) then
1246  !$omp parallel do private(i,j,k)
1247  !$acc kernels
1248  do j = js, je
1249  do i = is, ie
1250  do k = ks, ke
1251  w3d(k,i,j) = lt_path_tot(k,i,j,3)/dt_lt ![num/grid/s]
1252  enddo
1253  enddo
1254  enddo
1255  !$acc end kernels
1256  call file_history_put( hist_id(i_ltpath), w3d(:,:,:) )
1257  endif
1258  if ( hist_sw(i_posflash) ) then
1259  !$omp parallel do private(i,j,k)
1260  !$acc kernels
1261  do j = js, je
1262  do i = is, ie
1263  do k = ks, ke
1264  w3d(k,i,j) = lt_path_tot(k,i,j,1)/dt_lt ![num/grid/s]
1265  enddo
1266  enddo
1267  enddo
1268  !$acc end kernels
1269  call file_history_put( hist_id(i_posflash), w3d(:,:,:) )
1270  endif
1271  if ( hist_sw(i_negflash) ) then
1272  !$omp parallel do private(i,j,k)
1273  !$acc kernels
1274  do j = js, je
1275  do i = is, ie
1276  do k = ks, ke
1277  w3d(k,i,j) = lt_path_tot(k,i,j,2)/dt_lt ![num/grid/s]
1278  enddo
1279  enddo
1280  enddo
1281  !$acc end kernels
1282  call file_history_put( hist_id(i_negflash), w3d(:,:,:) )
1283  endif
1284  if ( hist_sw(i_flashpoint) ) then
1285  !$omp parallel do private(i,j,k)
1286  !$acc kernels
1287  do j = js, je
1288  do i = is, ie
1289  do k = ks, ke
1290  w3d(k,i,j) = fls_int_p_tot(k,i,j)/dt_lt ![num/grid/s]
1291  enddo
1292  enddo
1293  enddo
1294  !$acc end kernels
1295  call file_history_put( hist_id(i_flashpoint), w3d(:,:,:) )
1296  end if
1297  if ( hist_sw(i_fod) ) then
1298  !$omp parallel do private(i,j)
1299  !$acc kernels
1300  do j = js, je
1301  do i = is, ie
1302  w3d(1,i,j) = b_f2013_tot(i,j)/dt_lt ![num/grid/s]
1303  enddo
1304  enddo
1305  !$acc end kernels
1306  call file_history_put( hist_id(i_fod), w3d(1,:,:) )
1307  end if
1308 
1309  !$acc end data
1310 
1311  return
1312  end subroutine atmos_phy_lt_sato2019_adjustment
1313  !-----------------------------------------------------------------------------
1316  !-----------------------------------------------------------------------------
1317  subroutine atmos_phy_lt_electric_field( &
1318  KA, KS, KE, & ! [IN]
1319  IA, IS, IE, & ! [IN]
1320  JA, JS, JE, & ! [IN]
1321  QCRG, & ! [IN]
1322  DENS, & ! [IN]
1323  RHOT, & ! [IN]
1324  E_pot, & ! [INOUT]
1325  Efield ) ! [INOUT]
1326  use scale_prc, only: &
1327  prc_abort, &
1328  prc_ismaster, &
1329  prc_myrank
1330  use scale_const, only: &
1331  eps => const_eps, &
1332  epsvac => const_epsvac, &
1333  epsair => const_epsair
1334  use scale_atmos_grid_cartesc, only: &
1335  rcdz => atmos_grid_cartesc_rcdz, &
1336  rcdx => atmos_grid_cartesc_rcdx, &
1337  rcdy => atmos_grid_cartesc_rcdy, &
1338  rfdz => atmos_grid_cartesc_rfdz, &
1339  rfdx => atmos_grid_cartesc_rfdx, &
1340  rfdy => atmos_grid_cartesc_rfdy
1341  use scale_atmos_grid_cartesc_metric, only: &
1342  mapf => atmos_grid_cartesc_metric_mapf, &
1343  gsqrt => atmos_grid_cartesc_metric_gsqrt, &
1344  j13g => atmos_grid_cartesc_metric_j13g, &
1345  j23g => atmos_grid_cartesc_metric_j23g, &
1346  j33g => atmos_grid_cartesc_metric_j33g
1347  use scale_atmos_grid_cartesc_index, only: &
1348  i_xyz, &
1349  i_xyw, &
1350  i_uyw, &
1351  i_xvw, &
1352  i_uyz, &
1353  i_xvz, &
1354  i_uvz, &
1355  i_xy, &
1356  i_uy, &
1357  i_xv, &
1358  i_uv
1359  use scale_prc_cartesc, only: &
1360  prc_has_w, &
1361  prc_has_e, &
1362  prc_has_n, &
1363  prc_has_s
1364  use scale_comm_cartesc, only: &
1365  comm_datatype, &
1366  comm_world, &
1367  comm_vars8, &
1368  comm_wait
1369  implicit none
1370 
1371  integer, intent(in) :: KA, KS, KE
1372  integer, intent(in) :: IA, IS, IE
1373  integer, intent(in) :: JA, JS, JE
1374  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
1375  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
1376  real(RP), intent(in) :: RHOT (KA,IA,JA) !-- density weighted potential temperature [K kg/m3]
1377  real(RP), intent(inout) :: E_pot (KA,IA,JA) !-- Electric potential [V]
1378  real(RP), intent(inout) :: Efield(KA,IA,JA,3) !-- Electric field [V/m]
1379 
1380  real(RP) :: eps_air
1381  !--- A x E_pott = - QCRG/epsiron
1382  real(RP) :: B(KA,IA,JA) !--- B : -QCRG*DENS/epsiron
1383  real(RP) :: E_pot_N(KA,IA,JA) !--- electrical potential calculated by Bi-CGSTAB
1384 
1385  integer :: i, j, k, ijk, ierror
1386  real(RP) :: iprod, buf
1387 
1388  call prof_rapstart('LT_E_field', 2)
1389 
1390 
1391  !$acc data &
1392  !$acc copy(E_pot,Efield) &
1393  !$acc copyin(QCRG,DENS,RHOT) &
1394  !$acc create(B,E_pot_N) &
1395  !$acc copyin(GSQRT,MAPF,J13G,J23G,RCDX,RCDY,RFDZ)
1396 
1397  !$omp parallel do private(i,j,k)
1398  !$acc kernels
1399  do j = js, je
1400  do i = is, ie
1401  do k = ks, ke
1402  e_pot_n(k,i,j) = e_pot(k,i,j) !-- initial value -> previous step value
1403  end do
1404  end do
1405  end do
1406  !$acc end kernels
1407 
1408  call comm_vars8( e_pot_n, 1 )
1409 
1410 
1411  !$omp parallel do private(i,j,k,eps_air)
1412  !$acc kernels
1413  do j = js, je
1414  do i = is, ie
1415  do k = ks, ke
1416  eps_air = epsvac * epsair !--- temporary, dependency of epsiron on P and T will be implemented
1417  b(k,i,j) = - qcrg(k,i,j)/eps_air * 1.0e-9_rp ! [nC/m3] -> [C/m3 * m/F] = [V/m2]
1418  enddo
1419  enddo
1420  enddo
1421  !$acc end kernels
1422 
1423  call comm_wait ( e_pot_n, 1 )
1424 
1425  !--- calcuclate counter matrix
1426  call atmos_phy_lt_solve_bicgstab( &
1427  ka, ks, ke, & ! (in)
1428  ia, is, ie, & ! (in)
1429  ja, js, je, & ! (in)
1430  e_pot, & ! (out)
1431  e_pot_n, & ! (in)
1432  a, b ) ! (in)
1433 
1434  call comm_vars8( e_pot, 1 )
1435  call comm_wait ( e_pot, 1, .true. )
1436 
1437  !$omp parallel do private(i,j)
1438  !$acc parallel vector_length(32)
1439  !$acc loop collapse(2)
1440  do j = 1, ja
1441  do i = 1, ia
1442  do k = 1, ks-1
1443  e_pot(k,i,j) = 0.0_rp
1444  enddo
1445  do k = ke+1, ka
1446  e_pot(k,i,j) = 0.0_rp
1447  enddo
1448  enddo
1449  enddo
1450  !$acc end parallel
1451 
1452  !---- Calculate Electrical Field
1453  !$omp parallel do private(i,j,k)
1454  !$acc kernels
1455  do j = js, je
1456  do i = is, ie
1457  do k = ks, ke
1458  efield(k,i,j,i_lt_x) = - mapf(i,j,1,i_xyz)/gsqrt(k,i,j,i_xyz) * &
1459  ( &
1460  ( 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) ) &
1461  * rcdx(i) * 0.5_rp &
1462  + ( j13g(k+1,i,j,i_xyz)*gsqrt(k+1,i,j,i_xyz)*e_pot(k+1,i,j) &
1463  - j13g(k-1,i,j,i_xyz)*gsqrt(k-1,i,j,i_xyz)*e_pot(k-1,i,j) &
1464  ) &
1465  * rfdz(k) * 0.5_rp &
1466  )
1467  efield(k,i,j,i_lt_y) = - mapf(i,j,2,i_xyz)/gsqrt(k,i,j,i_xyz) * &
1468  ( &
1469  ( 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) ) &
1470  * rcdy(j) * 0.5_rp &
1471  + ( j23g(k+1,i,j,i_xyz)*gsqrt(k+1,i,j,i_xyz)*e_pot(k+1,i,j) &
1472  - j23g(k-1,i,j,i_xyz)*gsqrt(k-1,i,j,i_xyz)*e_pot(k-1,i,j) &
1473  ) &
1474  * rfdz(k) * 0.5_rp &
1475  )
1476  efield(k,i,j,i_lt_z) = - 1.0_rp/gsqrt(k,i,j,i_xyz) * &
1477  ( j33g*e_pot(k+1,i ,j )*gsqrt(k+1,i,j,i_xyz) &
1478  - j33g*e_pot(k-1,i ,j )*gsqrt(k-1,i,j,i_xyz) &
1479  ) &
1480  * rfdz(k) * 0.5_rp
1481  enddo
1482  enddo
1483  enddo
1484  !$acc end kernels
1485 
1486  !$acc end data
1487 
1488  call prof_rapend('LT_E_field', 2)
1489 
1490  return
1491  end subroutine atmos_phy_lt_electric_field
1492 
1493  subroutine atmos_phy_lt_solve_bicgstab( &
1494  KA, KS, KE, & ! [IN]
1495  IA, IS, IE, & ! [IN]
1496  JA, JS, JE, & ! [IN]
1497  PHI_N, & ! [OUT]
1498  PHI, & ! [IN]
1499  M, B ) ! [IN]
1500  use scale_prc, only: &
1501  prc_abort, &
1502  prc_ismaster
1503  use scale_comm_cartesc, only: &
1504  comm_datatype, &
1505  comm_world, &
1506  comm_vars8, &
1507  comm_wait
1508  implicit none
1509  integer, intent(in) :: KA, KS, KE
1510  integer, intent(in) :: IA, IS, IE
1511  integer, intent(in) :: JA, JS, JE
1512  real(RP), intent(out) :: PHI_N(KA,IA,JA)
1513  real(RP), intent(in) :: PHI(KA,IA,JA)
1514  real(RP), intent(in) :: M(KA,15,IA,JA)
1515  real(RP), intent(in) :: B(KA,IA,JA)
1516 
1517  real(RP) :: r0(KA,IA,JA)
1518 
1519  real(RP) :: p(KA,IA,JA)
1520  real(RP) :: Mp(KA,IA,JA)
1521  real(RP) :: s(KA,IA,JA)
1522  real(RP) :: Ms(KA,IA,JA)
1523  real(RP) :: al, be, w
1524 
1525  real(RP), pointer :: r(:,:,:), rn(:,:,:), swap(:,:,:)
1526  real(RP), target :: rbuf1(KA,IA,JA)
1527  real(RP), target :: rbuf2(KA,IA,JA)
1528  real(RP) :: v1(KA,IA,JA)
1529 
1530  real(RP) :: r0r
1531  real(RP) :: norm, error, error2
1532 
1533  real(RP) :: iprod1, iprod2
1534  real(RP) :: buf(2)
1535 
1536  real(RP):: diag(KA,IA,JA)
1537  real(RP):: z1(KA,IA,JA)
1538  real(RP):: z2(KA,IA,JA)
1539  real(RP):: Mz1(KA,IA,JA)
1540  real(RP):: Mz2(KA,IA,JA)
1541 
1542  integer :: k, i, j
1543  integer :: iis, iie, jjs, jje
1544  integer :: iter
1545  integer :: ierror
1546 
1547  !$acc data &
1548  !$acc copyout(PHI_N) &
1549  !$acc copyin(PHI,M,B) &
1550  !$acc create(r0,p,Mp,s,Ms,rbuf1,rbuf2,v1,diag,z1,z2,Mz1,Mz2)
1551 
1552  r => rbuf1
1553  rn => rbuf2
1554 
1555  if( flag_preprocessing == 3 ) then
1556  !$acc update host(M)
1557  call ilu_decomp(ka, ks, ke, &
1558  ia, is, ie, &
1559  ja, js, je, &
1560  m, diag)
1561  !$acc update device(diag)
1562  endif
1563 
1564  call mul_matrix( ka, ks, ke, & ! (in)
1565  ia, is, ie, & ! (in)
1566  ja, js, je, & ! (in)
1567  v1, m, phi ) ! v1 = M x0
1568 
1569  norm = 0.0_rp
1570  !$omp parallel do reduction(+:norm) private(i, j, k)
1571  !$acc kernels
1572  !$acc loop collapse(3) reduction(+:norm)
1573  do j = js, je
1574  do i = is, ie
1575  do k = ks, ke
1576  norm = norm + b(k,i,j)**2
1577  enddo
1578  enddo
1579  enddo
1580  !$acc end kernels
1581 
1582  ! r = b - M x0
1583  !$omp parallel do private(i, j, k)
1584  !$acc kernels
1585  do j = js, je
1586  do i = is, ie
1587  do k = ks, ke
1588  r(k,i,j) = b(k,i,j) - v1(k,i,j)
1589  enddo
1590  enddo
1591  enddo
1592  !$acc end kernels
1593 
1594  !$omp parallel do private(i, j, k)
1595  !$acc kernels
1596  do j = js, je
1597  do i = is, ie
1598  do k = ks, ke
1599  r0(k,i,j) = r(k,i,j)
1600  p(k,i,j) = r(k,i,j)
1601  enddo
1602  enddo
1603  enddo
1604  !$acc end kernels
1605 
1606  r0r = 0.0_rp
1607  !$omp parallel do reduction(+:r0r) private(i, j, k)
1608  !$acc kernels
1609  !$acc loop collapse(3) reduction(+:r0r)
1610  do j = js, je
1611  do i = is, ie
1612  do k = ks, ke
1613  r0r = r0r + r0(k,i,j) * r(k,i,j)
1614  enddo
1615  enddo
1616  enddo
1617  !$acc end kernels
1618 
1619  !$omp parallel do private(i, j, k)
1620  !$acc kernels
1621  do j = js-1, je+1
1622  do i = is-1, ie+1
1623  do k = ks, ke
1624  phi_n(k,i,j) = phi(k,i,j)
1625  end do
1626  end do
1627  end do
1628  !$acc end kernels
1629 
1630  buf(1) = r0r
1631  buf(2) = norm
1632  call mpi_allreduce(mpi_in_place, buf(:), 2, comm_datatype, mpi_sum, comm_world, ierror)
1633  r0r = buf(1)
1634  norm = buf(2)
1635  error2 = norm
1636 
1637  do iter = 1, itmax
1638 
1639  call comm_vars8( p, 1 )
1640 
1641  error = 0.0_rp
1642  !$omp parallel do reduction(+:error) private(i, j, k)
1643  !$acc kernels
1644  !$acc loop collapse(3) reduction(+:error)
1645  do j = js, je
1646  do i = is, ie
1647  do k = ks, ke
1648  error = error + r(k,i,j)**2
1649  enddo
1650  enddo
1651  enddo
1652  !$acc end kernels
1653 
1654  call mpi_allreduce(mpi_in_place, error, 1, comm_datatype, mpi_sum, comm_world, ierror)
1655 
1656  call comm_wait ( p, 1 )
1657 
1658  if ( sqrt(error/norm) < epsilon ) then
1659  log_info("ATMOS_PHY_LT_Efield",'(a,1x,i0,1x,2e15.7)') "Bi-CGSTAB converged:", iter, sqrt(error/norm),norm
1660  exit
1661  endif
1662  error2 = error
1663 
1664  if( flag_preprocessing == 0 ) then !-- No preprocessing
1665 
1666  call mul_matrix( ka, ks, ke, & ! (in)
1667  ia, is, ie, & ! (in)
1668  ja, js, je, & ! (in)
1669  mp, m, p )
1670 
1671  iprod1 = 0.0_rp
1672  !$omp parallel do reduction(+:iprod1) private(i, j, k)
1673  !$acc kernels
1674  !$acc loop collapse(3) reduction(+:iprod1)
1675  do j = js, je
1676  do i = is, ie
1677  do k = ks, ke
1678  iprod1 = iprod1 + r0(k,i,j) * mp(k,i,j)
1679  enddo
1680  enddo
1681  enddo
1682  !$acc end kernels
1683 
1684  else !--- Preprocessing
1685 
1686  if( flag_preprocessing == 1 ) then !--- Gauss-Seidel preprocessing
1687 
1688  call gs( ka, ks, ke, & ! (in)
1689  ia, is, ie, & ! (in)
1690  ja, js, je, & ! (in)
1691  z1, m, p )
1692 
1693  elseif( flag_preprocessing == 2 ) then !--- Synmetric Gauss-Seidel preprocessing (Default)
1694 
1695  call sgs( ka, ks, ke, & ! (in)
1696  ia, is, ie, & ! (in)
1697  ja, js, je, & ! (in)
1698  z1, m, p )
1699 
1700  elseif( flag_preprocessing == 3 ) then !--- Incomplete Cholesky Factorization preprocessing
1701 
1702  !$acc update host(p,diag) ! M is already updated
1703  call solve_ilu( ka, ks, ke, & ! (in)
1704  ia, is, ie, & ! (in)
1705  ja, js, je, & ! (in)
1706  z1, m, p, diag)
1707  !$acc update device(z1)
1708 
1709  endif
1710 
1711  call comm_vars8( z1, 1 )
1712 
1713  call mul_matrix( ka, ks, ke, & ! (in)
1714  ia, is, ie, & ! (in)
1715  ja, js, je, & ! (in)
1716  mp, m, p )
1717 
1718  call comm_wait ( z1, 1 )
1719 
1720  call mul_matrix( ka, ks, ke, & ! (in)
1721  ia, is, ie, & ! (in)
1722  ja, js, je, & ! (in)
1723  mz1, m, z1 )
1724 
1725  iprod1 = 0.0_rp
1726  !$omp parallel do reduction(+:iprod1) private(i, j, k)
1727  !$acc kernels
1728  !$acc loop collapse(3) reduction(+:iprod1)
1729  do j = js, je
1730  do i = is, ie
1731  do k = ks, ke
1732  iprod1 = iprod1 + r0(k,i,j) * mz1(k,i,j)
1733  enddo
1734  enddo
1735  enddo
1736  !$acc end kernels
1737 
1738  endif
1739 
1740  call mpi_allreduce(mpi_in_place, iprod1, 1, comm_datatype, mpi_sum, comm_world, ierror)
1741 
1742  if ( iprod1 == 0.0_rp ) then
1743  log_info("ATMOS_PHY_LT_Efield",'(a,1x,e15.7,1x,i10)') 'Iprod1 is zero(Bi-CGSTAB) skip:', iprod1, iter
1744  exit
1745  endif
1746  al = r0r / iprod1 ! (r0,r) / (r0,Mp)
1747 
1748  if( flag_preprocessing == 0 ) then !-- No preprocessing
1749  !$omp parallel do
1750  !$acc kernels
1751  do j = js, je
1752  do i = is, ie
1753  do k = ks, ke
1754  s(k,i,j) = r(k,i,j) - al*mp(k,i,j)
1755  enddo
1756  enddo
1757  enddo
1758  !$acc end kernels
1759  else ! Preprocessing
1760  !$omp parallel do private(i, j, k)
1761  !$acc kernels
1762  do j = js, je
1763  do i = is, ie
1764  do k = ks, ke
1765  s(k,i,j) = r(k,i,j) - al*mz1(k,i,j)
1766  enddo
1767  enddo
1768  enddo
1769  !$acc end kernels
1770  endif
1771 
1772  call comm_vars8( s, 1 )
1773  call comm_wait ( s, 1 )
1774  if( flag_preprocessing == 0 ) then !--- No Preprocessing
1775 
1776  call mul_matrix( ka, ks, ke, & ! (in)
1777  ia, is, ie, & ! (in)
1778  ja, js, je, & ! (in)
1779  ms, m, s )
1780 
1781  iprod1 = 0.0_rp
1782  iprod2 = 0.0_rp
1783  !$omp parallel do reduction(+:iprod1,iprod2)
1784  !$acc kernels
1785  !$acc loop collapse(3) reduction(+:iprod1,iprod2)
1786  do j = js, je
1787  do i = is, ie
1788  do k = ks, ke
1789  iprod1 = iprod1 + ms(k,i,j) * s(k,i,j)
1790  iprod2 = iprod2 + ms(k,i,j) * ms(k,i,j)
1791  enddo
1792  enddo
1793  enddo
1794  !$acc end kernels
1795 
1796  else !--- Preprocessing
1797 
1798  if( flag_preprocessing == 1 ) then !--- Gauss-Seidel preprocessing
1799 
1800  call gs( ka, ks, ke, & ! (in)
1801  ia, is, ie, & ! (in)
1802  ja, js, je, & ! (in)
1803  z2, m, s )
1804 
1805  elseif( flag_preprocessing == 2 ) then !--- Synmetric Gauss-Seidel preprocessing (Default)
1806 
1807  call sgs( ka, is, ke, & ! (in)
1808  ia, is, ie, & ! (in)
1809  ja, js, je, & ! (in)
1810  z2, m, s )
1811 
1812  elseif( flag_preprocessing == 3 ) then !--- Incomplete Cholesky Factorization preprocessing
1813 
1814  !$acc update host(s,diag) ! M is already updated
1815  call solve_ilu( ka, ks, ke, & ! (in)
1816  ia, is, ie, & ! (in)
1817  ja, js, je, & ! (in)
1818  z2, m, s, diag)
1819  !$acc update device(z2)
1820  endif
1821 
1822  call comm_vars8( z2, 1 )
1823  call comm_wait ( z2, 1 )
1824  call mul_matrix( ka, ks, ke, & ! (in)
1825  ia, is, ie, & ! (in)
1826  ja, js, je, & ! (in)
1827  mz2, m, z2 )
1828 
1829  iprod1 = 0.0_rp
1830  iprod2 = 0.0_rp
1831  !$omp parallel do reduction(+:iprod1,iprod2) private(i, j, k)
1832  !$acc kernels
1833  !$acc loop collapse(3) reduction(+:iprod1,iprod2)
1834  do j = js, je
1835  do i = is, ie
1836  do k = ks, ke
1837  iprod1 = iprod1 + mz2(k,i,j) * s(k,i,j)
1838  iprod2 = iprod2 + mz2(k,i,j) * mz2(k,i,j)
1839  enddo
1840  enddo
1841  enddo
1842  !$acc end kernels
1843 
1844  endif
1845 
1846  buf(1) = iprod1
1847  buf(2) = iprod2
1848  call mpi_allreduce(mpi_in_place, buf(:), 2, comm_datatype, mpi_sum, comm_world, ierror)
1849 
1850  if ( buf(2) == 0.0_rp ) then
1851  log_info("ATMOS_PHY_LT_Efield",'(a,1x,e15.7,1x,i10)') 'Buf(2) is zero(Bi-CGSTAB) skip:', buf(2), iter
1852  exit
1853  endif
1854  w = buf(1) / buf(2) ! (Ms,s) / (Ms,Ms)
1855 
1856  if( flag_preprocessing == 0 ) then !--- No preprocessing
1857 
1858  !$omp parallel do private(i, j, k)
1859  !$acc kernels
1860  do j = js, je
1861  do i = is, ie
1862  do k = ks, ke
1863  phi_n(k,i,j) = phi_n(k,i,j) + al*p(k,i,j) + w*s(k,i,j)
1864  enddo
1865  enddo
1866  enddo
1867  !$acc end kernels
1868 
1869  !$omp parallel do private(i, j, k)
1870  !$acc kernels
1871  do j = js, je
1872  do i = is, ie
1873  do k = ks, ke
1874  rn(k,i,j) = s(k,i,j) - w*ms(k,i,j)
1875  enddo
1876  enddo
1877  enddo
1878  !$acc end kernels
1879 
1880  else
1881 
1882  !$omp parallel do private(i, j, k)
1883  !$acc kernels
1884  do j = js, je
1885  do i = is, ie
1886  do k = ks, ke
1887  phi_n(k,i,j) = phi_n(k,i,j) + al*z1(k,i,j) + w*z2(k,i,j)
1888  enddo
1889  enddo
1890  enddo
1891  !$acc end kernels
1892 
1893  !$omp parallel do private(i, j, k)
1894  !$acc kernels
1895  do j = js, je
1896  do i = is, ie
1897  do k = ks, ke
1898  rn(k,i,j) = s(k,i,j) - w*mz2(k,i,j)
1899  enddo
1900  enddo
1901  enddo
1902  !$acc end kernels
1903 
1904  endif
1905 
1906  iprod1 = 0.0_rp
1907  !$omp parallel do reduction(+:iprod1) private(i, j, k)
1908  !$acc kernels
1909  !$acc loop collapse(3) reduction(+:iprod1)
1910  do j = js, je
1911  do i = is, ie
1912  do k = ks, ke
1913  iprod1 = iprod1 + r0(k,i,j) * rn(k,i,j)
1914  enddo
1915  enddo
1916  enddo
1917  !$acc end kernels
1918 
1919  be = al/w / r0r
1920 
1921  call mpi_allreduce(iprod1, r0r, 1, comm_datatype, mpi_sum, comm_world, ierror)
1922 
1923  be = be * r0r ! al/w * (r0,rn)/(r0,r)
1924 
1925  if( flag_preprocessing == 0 ) then !--- No preprocessing
1926  !$omp parallel do private(i, j, k)
1927  !$acc kernels
1928  do j = js, je
1929  do i = is, ie
1930  do k = ks, ke
1931  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mp(k,i,j) )
1932  enddo
1933  enddo
1934  enddo
1935  !$acc end kernels
1936  else
1937  !$omp parallel do private(i, j, k)
1938  !$acc kernels
1939  do j = js, je
1940  do i = is, ie
1941  do k = ks, ke
1942  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mz1(k,i,j) )
1943  enddo
1944  enddo
1945  enddo
1946  !$acc end kernels
1947  endif
1948 
1949  swap => rn
1950  rn => r
1951  r => swap
1952 
1953  if ( r0r == 0.0_rp ) then
1954  log_info("ATMOS_PHY_LT_Efield",'(a,1x,i0,1x,3e15.7)') "Inner product of r0 and r_itr is zero(Bi-CGSTAB) :", &
1955  iter, r0r, sqrt(error/norm), norm
1956  exit
1957  endif
1958 
1959  enddo
1960 
1961  if ( iter >= itmax ) then
1962  if( prc_ismaster ) then
1963  log_warn("ATMOS_PHY_LT_solve_bicgstab",'(a,1x,2e15.7)') 'Bi-CGSTAB not converged:', error, norm
1964  log_warn_cont('(a,1x,2e15.7)') 'Bi-CGSTAB not converged:', epsilon, sqrt(error/norm)
1965  log_warn_cont('(a,1x,2e15.7)') 'xxx epsilon(set,last)=', epsilon, sqrt(error/norm)
1966  if( error /= error ) then
1967  log_error("ATMOS_PHY_LT_solve_bicgstab",*) 'xxx error or norm is NaN Stop!'
1968  call prc_abort
1969  endif
1970  endif
1971  endif
1972 
1973  !$acc end data
1974  return
1975  end subroutine atmos_phy_lt_solve_bicgstab
1976 
1977  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
1978  subroutine ilu_decomp(KA, KS, KE, &
1979  IA, IS, IE, &
1980  JA, JS, JE, &
1981  M, diag)
1982  implicit none
1983  integer, intent(in) :: KA, KS, KE
1984  integer, intent(in) :: IA, IS, IE
1985  integer, intent(in) :: JA, JS, JE
1986  real(RP), intent(in) :: M(KA,15,IA,JA)
1987  real(RP), intent(out) :: diag(KA,IA,JA)
1988 
1989  integer :: k, i, j
1990 
1991  !$omp parallel do private(i, j, k)
1992  do j=js,je
1993  do i=is,ie
1994  do k=ks,ke
1995  diag(k,i,j)=m(k,1,i,j)
1996  enddo
1997  enddo
1998  enddo
1999 
2000  !$omp parallel do private(i, j, k)
2001  do j=js,je
2002  do i=is,ie
2003  do k=ks,ke
2004  if(k /= ke) then
2005  diag(k+1,i,j) = diag(k+1,i,j) + m(k,3,i,j)*m(k+1,2,i,j)/diag(k,i,j)
2006  if(i /= ie) then
2007  diag(k+1,i+1,j) = diag(k+1,i+1,j) + m(k,13,i,j)*m(k+1,8,i+1,j)/diag(k,i,j)
2008  endif
2009  if(j /= je) then
2010  diag(k+1,i,j+1) = diag(k+1,i,j+1) + m(k,15,i,j)*m(k+1,10,i,j+1)/diag(k,i,j)
2011  endif
2012  endif
2013  if(i /= ie) then
2014  diag(k,i+1,j) = diag(k,i+1,j) +m(k,5,i,j)*m(k,4,i+1,j)/diag(k,i,j)
2015  if(k /= ks) then
2016  diag(k-1,i+1,j) = diag(k-1,i+1,j) + m(k,9,i,j)*m(k-1,12,i+1,j)/diag(k,i,j)
2017  endif
2018  endif
2019  if(j /= je) then
2020  diag(k,i,j+1)=diag(k,i,j+1) + m(k,7,i,j)*m(k,6,i,j+1)/diag(k,i,j)
2021  if(k /= ks) then
2022  diag(k-1,i,j+1) = diag(k-1,i,j+1) + m(k,11,i,j)*m(k-1,14,i,j+1)/diag(k,i,j)
2023  endif
2024  endif
2025  enddo
2026  enddo
2027  enddo
2028 
2029  return
2030  end subroutine ilu_decomp
2031 
2032  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2033  subroutine solve_ilu( KA,KS,KE, &
2034  IA,IS,IE, &
2035  JA,JS,JE, &
2036  Z, M, V, diag)
2037  use scale_prc, only: &
2038  prc_abort, &
2039  prc_ismaster
2040  implicit none
2041  integer, intent(in) :: KA, KS, KE
2042  integer, intent(in) :: IA, IS, IE
2043  integer, intent(in) :: JA, JS, JE
2044  real(RP), intent(in) :: M(KA,15,IA,JA)
2045  real(RP), intent(in) :: diag(KA,IA,JA)
2046  real(RP), intent(in) :: V(KA,IA,JA)
2047  real(RP), intent(out) :: Z(KA,IA,JA)
2048 
2049  real(RP):: Y(KA,IA,JA)
2050  integer :: k,i,j
2051 
2052  call gs_ilu(ka, ks, ke, &
2053  ia, is, ie, &
2054  ja, js, je, &
2055  y, m, v, diag)
2056 
2057  !$omp parallel do
2058  do j=js,je
2059  do i=is,ie
2060  do k=ks,ke
2061  y(k,i,j)=diag(k,i,j)*y(k,i,j)
2062  enddo
2063  enddo
2064  enddo
2065 
2066  call back_sub_ilu(ka, ks, ke, &
2067  ia, is, ie, &
2068  ja, js, je, &
2069  z, m, y, diag)
2070 
2071  return
2072 
2073  end subroutine solve_ilu
2074 
2075  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2076  subroutine back_sub_ilu( KA, KS, KE, &
2077  IA, IS, IE, &
2078  JA, JS, JE, &
2079  Z, M, V,diag)
2080  use scale_prc, only: &
2081  prc_abort, &
2082  prc_ismaster
2083  implicit none
2084  integer, intent(in) :: KA, KS, KE
2085  integer, intent(in) :: IA, IS, IE
2086  integer, intent(in) :: JA, JS, JE
2087  real(RP), intent(in) :: V(KA,IA,JA)
2088  real(RP), intent(in) :: M(KA,15,IA,JA)
2089  real(RP), intent(in) :: diag(KA,IA,JA)
2090  real(RP), intent(out) :: Z(KA,IA,JA)
2091 
2092  integer :: k, i, j
2093 
2094  z(:,:,:)=0.0_rp
2095 
2096  do j = je, js,-1
2097  do i = ie, is,-1
2098 
2099  z(ke,i,j) = ( v(ke,i,j) &
2100  -( m(ke,2,i,j) * z(ke-1,i ,j ) &
2101  + m(ke,4,i,j) * z(ke ,i-1,j ) &
2102  + m(ke,5,i,j) * z(ke ,i+1,j ) &
2103  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2104  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2105  + m(ke,8,i,j) * z(ke-1,i-1,j ) &
2106  + m(ke,9,i,j) * z(ke-1,i+1,j ) &
2107  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2108  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/diag(ke,i,j)
2109 
2110  do k = ke-1, ks+1,-1
2111  z(k,i,j) = (v(k,i,j) &
2112  -( m(k,2,i,j) * z(k-1,i ,j ) &
2113  + m(k,3,i,j) * z(k+1,i ,j ) &
2114  + m(k,4,i,j) * z(k ,i-1,j ) &
2115  + m(k,5,i,j) * z(k ,i+1,j ) &
2116  + m(k,6,i,j) * z(k ,i ,j-1) &
2117  + m(k,7,i,j) * z(k ,i ,j+1) &
2118  + m(k,8,i,j) * z(k-1,i-1,j ) &
2119  + m(k,9,i,j) * z(k-1,i+1,j ) &
2120  + m(k,10,i,j)* z(k-1,i ,j-1) &
2121  + m(k,11,i,j)* z(k-1,i ,j+1) &
2122  + m(k,12,i,j)* z(k+1,i-1,j ) &
2123  + m(k,13,i,j)* z(k+1,i+1,j ) &
2124  + m(k,14,i,j)* z(k+1,i ,j-1) &
2125  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/diag(k,i,j)
2126  enddo
2127 
2128  z(ks,i,j) = (v(ks,i,j) &
2129  -( m(ks,3,i,j) * z(ks+1,i ,j ) &
2130  + m(ks,4,i,j) * z(ks ,i-1,j ) &
2131  + m(ks,5,i,j) * z(ks ,i+1,j ) &
2132  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2133  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2134  + m(ks,12,i,j)* z(ks+1,i-1,j ) &
2135  + m(ks,13,i,j)* z(ks+1,i+1,j ) &
2136  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2137  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /diag(ks,i,j)
2138 
2139  enddo
2140  enddo
2141 
2142  return
2143  end subroutine back_sub_ilu
2144 
2145  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2146  subroutine gs_ilu( KA, KS, KE, &
2147  IA, IS, IE, &
2148  JA, JS, JE, &
2149  Z, M, V , diag) !M*Z=V -> Z=M^{-1} V
2150 
2151  implicit none
2152  integer, intent(in) :: KA, KS, KE
2153  integer, intent(in) :: IA, IS, IE
2154  integer, intent(in) :: JA, JS, JE
2155  real(RP), intent(in) :: V(KA,IA,JA)
2156  real(RP), intent(in) :: M(KA,15,IA,JA)
2157  real(RP), intent(in) :: diag(KA,IA,JA)
2158  real(RP), intent(out) :: Z(KA,IA,JA)
2159 
2160  integer :: k, i, j
2161 
2162  z(:,:,:)=0.0_rp
2163 
2164  do j = js, je
2165  do i = is, ie
2166 
2167  z(ks,i,j) = (v(ks,i,j) &
2168  -( m(ks,3,i,j) * z(ks+1,i ,j ) &
2169  + m(ks,4,i,j) * z(ks ,i-1,j ) &
2170  + m(ks,5,i,j) * z(ks ,i+1,j ) &
2171  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2172  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2173  + m(ks,12,i,j)* z(ks+1,i-1,j ) &
2174  + m(ks,13,i,j)* z(ks+1,i+1,j ) &
2175  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2176  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /diag(ks,i,j)
2177 
2178  do k = ks+1, ke-1
2179  z(k,i,j) = (v(k,i,j) &
2180  -( m(k,2,i,j) * z(k-1,i ,j ) &
2181  + m(k,3,i,j) * z(k+1,i ,j ) &
2182  + m(k,4,i,j) * z(k ,i-1,j ) &
2183  + m(k,5,i,j) * z(k ,i+1,j ) &
2184  + m(k,6,i,j) * z(k ,i ,j-1) &
2185  + m(k,7,i,j) * z(k ,i ,j+1) &
2186  + m(k,8,i,j) * z(k-1,i-1,j ) &
2187  + m(k,9,i,j) * z(k-1,i+1,j ) &
2188  + m(k,10,i,j)* z(k-1,i ,j-1) &
2189  + m(k,11,i,j)* z(k-1,i ,j+1) &
2190  + m(k,12,i,j)* z(k+1,i-1,j ) &
2191  + m(k,13,i,j)* z(k+1,i+1,j ) &
2192  + m(k,14,i,j)* z(k+1,i ,j-1) &
2193  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/diag(k,i,j)
2194  enddo
2195 
2196  z(ke,i,j) = ( v(ke,i,j) &
2197  -( m(ke,2,i,j) * z(ke-1,i ,j ) &
2198  + m(ke,4,i,j) * z(ke ,i-1,j ) &
2199  + m(ke,5,i,j) * z(ke ,i+1,j ) &
2200  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2201  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2202  + m(ke,8,i,j) * z(ke-1,i-1,j ) &
2203  + m(ke,9,i,j) * z(ke-1,i+1,j ) &
2204  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2205  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/diag(ke,i,j)
2206 
2207  enddo
2208  enddo
2209 
2210  return
2211  end subroutine gs_ilu
2212 
2213  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2214  subroutine sgs( KA, KS, KE, &
2215  IA, IS, IE, &
2216  JA, JS, JE, &
2217  Z, M, V)
2218  implicit none
2219  integer, intent(in) :: KA, KS, KE
2220  integer, intent(in) :: IA, IS, IE
2221  integer, intent(in) :: JA, JS, JE
2222  real(RP), intent(in) :: M(KA,15,IA,JA)
2223  real(RP), intent(in) :: V(KA,IA,JA)
2224  real(RP), intent(out) :: Z(KA,IA,JA)
2225  real(RP):: Y(KA,IA,JA)
2226  integer :: k,i,j
2227 
2228  !$acc data copyin(M,V) copyout(Z) create(Y)
2229 
2230  call gs( ka, ks, ke, &
2231  ia, is, ie, &
2232  ja, js, je, &
2233  y, m, v )
2234 
2235  !$omp parallel do private(i, j, k)
2236  !$acc kernels
2237  do j = js, je
2238  do i = is, ie
2239  do k = ks, ke
2240  y(k,i,j)=m(k,1,i,j)*y(k,i,j)
2241  enddo
2242  enddo
2243  enddo
2244  !$acc end kernels
2245 
2246  call back_sub(ka, ks, ke, &
2247  ia, is, ie, &
2248  ja, js, je, &
2249  z, m, y )
2250 
2251  !$acc end data
2252 
2253  return
2254  end subroutine sgs
2255 
2256  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2257  subroutine back_sub ( KA, KS, KE, &
2258  IA, IS, IE, &
2259  JA, JS, JE, &
2260  Z, M, V )
2261  implicit none
2262  integer, intent(in) :: KA, KS, KE
2263  integer, intent(in) :: IA, IS, IE
2264  integer, intent(in) :: JA, JS, JE
2265  real(RP), intent(in) :: V(KA,IA,JA)
2266  real(RP), intent(in) :: M(KA,15,IA,JA)
2267  real(RP), intent(out) :: Z(KA,IA,JA)
2268  integer :: k, i, j
2269 
2270  !$acc data copyin(V,M) copyout(Z)
2271 
2272 !$acc kernels
2273 z(:,:,:)=1d30
2274 !$acc end kernels
2275 
2276  !$omp parallel
2277  !$omp do
2278  !$acc parallel vector_length(32)
2279  do i = is, ie
2280  do k = ks, ke
2281  z(k,i,js-1) = 0.0_rp
2282  z(k,i,je+1) = 0.0_rp
2283  end do
2284  end do
2285  !$acc end parallel
2286  !$omp do
2287  !$acc parallel vector_length(32)
2288  do j = js, je
2289  do k = ks, ke
2290  z(k,is-1,j) = 0.0_rp
2291  z(k,ie+1,j) = 0.0_rp
2292  end do
2293  end do
2294  !$acc end parallel
2295  !$omp end parallel
2296 
2297 
2298 #if COLORING == 2
2299  !$omp parallel do collapse(2)
2300 ! !$acc parallel vector_length(32)
2301  !$acc parallel vector_length(1)
2302  !$acc loop collapse(2)
2303  do j = je, js, -1
2304  do i = ie, is, -1
2305  if ( mod((ie-i)+(je-j),2)==0 ) then
2306 
2307  z(ke,i,j) = ( v(ke,i,j) )/m(ke,1,i,j)
2308 
2309  do k = ke-1, ks, -1
2310  z(k,i,j) = (v(k,i,j) &
2311  -( m(k,3,i,j) * z(k+1,i ,j ) ) )/m(k,1,i,j)
2312  enddo
2313 
2314  end if
2315 
2316 #elif COLORING == 1
2317  !$omp parallel do
2318  !$acc parallel vector_length(1)
2319  do j = je, js, -2
2320  !$acc loop seq
2321  do i = ie, is, -1
2322 
2323  z(ke,i,j) = ( v(ke,i,j) &
2324  -( m(ke,5,i,j) * z(ke ,i+1,j ) ) )/m(ke,1,i,j)
2325 
2326  do k = ke-1, ks, -1
2327  z(k,i,j) = (v(k,i,j) &
2328  -( m(k,3,i,j) * z(k+1,i ,j ) &
2329  + m(k,5,i,j) * z(k ,i+1,j ) &
2330  + m(k,13,i,j)* z(k+1,i+1,j ) ) )/m(k,1,i,j)
2331  enddo
2332 #else
2333  !$acc parallel
2334  !$acc loop seq
2335  do j = je, js, -1
2336  !$acc loop seq
2337  do i = ie, is, -1
2338 
2339  z(ke,i,j) = ( v(ke,i,j) &
2340  -( m(ke,5,i,j) * z(ke ,i+1,j ) &
2341  + m(ke,7,i,j) * z(ke ,i ,j+1) ) )/m(ke,1,i,j)
2342 
2343  do k = ke-1, ks, -1
2344  z(k,i,j) = (v(k,i,j) &
2345  -( m(k,3,i,j) * z(k+1,i ,j ) &
2346  + m(k,5,i,j) * z(k ,i+1,j ) &
2347  + m(k,7,i,j) * z(k ,i ,j+1) &
2348  + m(k,13,i,j)* z(k+1,i+1,j ) &
2349  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/m(k,1,i,j)
2350  enddo
2351 #endif
2352  end do
2353  end do
2354  !$acc end parallel
2355 
2356 #if COLORING > 0
2357 #if COLORING == 2
2358  !$omp parallel do collapse(2)
2359 ! !$acc parallel vector_length(32)
2360  !$acc parallel vector_length(1)
2361  !$acc loop collapse(2)
2362  do j = je, js, -1
2363  do i = ie, is, -1
2364  if ( mod((ie-i)+(je-j),2)==1 ) then
2365 
2366  z(ke,i,j) = ( v(ke,i,j) &
2367  -( m(ke,4,i,j) * z(ke ,i-1,j ) &
2368  + m(ke,5,i,j) * z(ke ,i+1,j ) &
2369  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2370  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2371  + m(ke,8,i,j) * z(ke-1,i-1,j ) &
2372  + m(ke,9,i,j) * z(ke-1,i+1,j ) &
2373  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2374  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/m(ke,1,i,j)
2375 
2376  do k = ke-1, ks+1,-1
2377  z(k,i,j) = (v(k,i,j) &
2378  -( m(k,3,i,j) * z(k+1,i ,j ) &
2379  + m(k,4,i,j) * z(k ,i-1,j ) &
2380  + m(k,5,i,j) * z(k ,i+1,j ) &
2381  + m(k,6,i,j) * z(k ,i ,j-1) &
2382  + m(k,7,i,j) * z(k ,i ,j+1) &
2383  + m(k,8,i,j) * z(k-1,i-1,j ) &
2384  + m(k,9,i,j) * z(k-1,i+1,j ) &
2385  + m(k,10,i,j)* z(k-1,i ,j-1) &
2386  + m(k,11,i,j)* z(k-1,i ,j+1) &
2387  + m(k,12,i,j)* z(k+1,i-1,j ) &
2388  + m(k,13,i,j)* z(k+1,i+1,j ) &
2389  + m(k,14,i,j)* z(k+1,i ,j-1) &
2390  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/m(k,1,i,j)
2391  enddo
2392 
2393  z(ks,i,j) = (v(ks,i,j) &
2394  -( m(ks,3,i,j) * z(ks+1,i ,j ) &
2395  + m(ks,4,i,j) * z(ks ,i-1,j ) &
2396  + m(ks,5,i,j) * z(ks ,i+1,j ) &
2397  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2398  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2399  + m(ks,12,i,j)* z(ks+1,i-1,j ) &
2400  + m(ks,13,i,j)* z(ks+1,i+1,j ) &
2401  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2402  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /m(ks,1,i,j)
2403 
2404  end if
2405 #elif COLORING == 1
2406  !$omp parallel do
2407  !$acc parallel vector_length(1)
2408  do j = je-1, js, -2
2409  !$acc loop seq
2410  do i = ie, is, -1
2411  z(ke,i,j) = ( v(ke,i,j) &
2412  -( m(ke,5,i,j) * z(ke ,i+1,j ) &
2413  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2414  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2415  + m(ke,9,i,j) * z(ke-1,i+1,j ) &
2416  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2417  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/m(ke,1,i,j)
2418 
2419  do k = ke-1, ks+1,-1
2420  z(k,i,j) = (v(k,i,j) &
2421  -( m(k,3,i,j) * z(k+1,i ,j ) &
2422  + m(k,5,i,j) * z(k ,i+1,j ) &
2423  + m(k,6,i,j) * z(k ,i ,j-1) &
2424  + m(k,7,i,j) * z(k ,i ,j+1) &
2425  + m(k,9,i,j) * z(k-1,i+1,j ) &
2426  + m(k,10,i,j)* z(k-1,i ,j-1) &
2427  + m(k,11,i,j)* z(k-1,i ,j+1) &
2428  + m(k,13,i,j)* z(k+1,i+1,j ) &
2429  + m(k,14,i,j)* z(k+1,i ,j-1) &
2430  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/m(k,1,i,j)
2431  enddo
2432 
2433  z(ks,i,j) = (v(ks,i,j) &
2434  -( m(ks,3,i,j) * z(ks+1,i ,j ) &
2435  + m(ks,5,i,j) * z(ks ,i+1,j ) &
2436  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2437  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2438  + m(ks,13,i,j)* z(ks+1,i+1,j ) &
2439  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2440  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /m(ks,1,i,j)
2441 #endif
2442  end do
2443  end do
2444  !$acc end parallel
2445 #endif
2446 
2447  !$acc end data
2448 
2449  return
2450  end subroutine back_sub
2451 
2452  !----- N. Yamashita and T. Iwashita of Hokkaido Univ. created (2021/2/22)-----
2453  subroutine gs ( KA, KS, KE, &
2454  IA, IS, IE, &
2455  JA, JS, JE, &
2456  Z, M, V )
2457  implicit none
2458  integer, intent(in) :: KA, KS, KE
2459  integer, intent(in) :: IA, IS, IE
2460  integer, intent(in) :: JA, JS, JE
2461  real(RP), intent(in) :: V(KA,IA,JA)
2462  real(RP), intent(in) :: M(KA,15,IA,JA)
2463  real(RP), intent(out) :: Z(KA,IA,JA)
2464 
2465  integer :: k, i, j, n
2466 
2467  !$acc data copyin(V,M) copyout(Z)
2468 
2469  !$omp parallel
2470  !$omp do
2471  !$acc parallel vector_length(32)
2472  do i = is, ie
2473  do k = ks, ke
2474  z(k,i,js-1) = 0.0_rp
2475  z(k,i,je+1) = 0.0_rp
2476  end do
2477  end do
2478  !$acc end parallel
2479  !$omp do
2480  !$acc parallel vector_length(32)
2481  do j = js, je
2482  do k = ks, ke
2483  z(k,is-1,j) = 0.0_rp
2484  z(k,ie+1,j) = 0.0_rp
2485  end do
2486  end do
2487  !$acc end parallel
2488  !$omp end parallel
2489 
2490 
2491 #if COLORING == 2
2492  !$omp parallel do collapse(2)
2493 ! !$acc parallel vector_length(32)
2494  !$acc parallel vector_length(1)
2495  !$acc loop collapse(2)
2496  do j = js, je
2497  do i = is, ie
2498  if ( mod((i-is)+(j-js),2)==0 ) then
2499 
2500  z(ks,i,j) = (v(ks,i,j) ) /m(ks,1,i,j)
2501 
2502  do k = ks+1, ke
2503  z(k,i,j) = (v(k,i,j) &
2504  -( m(k,2,i,j) * z(k-1,i ,j ) ) )/m(k,1,i,j)
2505  enddo
2506 
2507  end if
2508 
2509 #elif COLORING == 1
2510  !$omp parallel do
2511  !$acc parallel vector_length(1)
2512  do j = js, je, 2
2513  !$acc loop seq
2514  do i = is, ie
2515 
2516  z(ks,i,j) = (v(ks,i,j) &
2517  -( m(ks,4,i,j) * z(ks ,i-1,j ) ) ) /m(ks,1,i,j)
2518 
2519  do k = ks+1, ke
2520  z(k,i,j) = (v(k,i,j) &
2521  -( m(k,2,i,j) * z(k-1,i ,j ) &
2522  + m(k,4,i,j) * z(k ,i-1,j ) &
2523  + m(k,8,i,j) * z(k-1,i-1,j ) ) )/m(k,1,i,j)
2524  enddo
2525 #else
2526  !$acc parallel
2527  !$acc loop seq
2528  do j = js, je
2529  !$acc loop seq
2530  do i = is, ie
2531 
2532  z(ks,i,j) = (v(ks,i,j) &
2533  -( m(ks,4,i,j) * z(ks ,i-1,j ) &
2534  + m(ks,6,i,j) * z(ks ,i ,j-1) ) ) /m(ks,1,i,j)
2535 
2536  do k = ks+1, ke
2537  z(k,i,j) = (v(k,i,j) &
2538  -( m(k,2,i,j) * z(k-1,i ,j ) &
2539  + m(k,4,i,j) * z(k ,i-1,j ) &
2540  + m(k,6,i,j) * z(k ,i ,j-1) &
2541  + m(k,8,i,j) * z(k-1,i-1,j ) &
2542  + m(k,10,i,j)* z(k-1,i ,j-1) ) )/m(k,1,i,j)
2543  enddo
2544 #endif
2545  end do
2546  end do
2547  !$acc end parallel
2548 
2549 #if COLORING > 0
2550 #if COLORING == 2
2551  !$omp parallel do collapse(2)
2552 ! !$acc parallel vector_length(32)
2553  !$acc parallel vector_length(1)
2554  !$acc loop collapse(2)
2555  do j = js, je
2556  do i = is, ie
2557  if ( mod((i-is)+(j-js),2)==1 ) then
2558 
2559  z(ks,i,j) = (v(ks,i,j) &
2560  -( m(ks,4,i,j) * z(ks ,i-1,j ) &
2561  + m(ks,5,i,j) * z(ks ,i+1,j ) &
2562  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2563  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2564  + m(ks,12,i,j)* z(ks+1,i-1,j ) &
2565  + m(ks,13,i,j)* z(ks+1,i+1,j ) &
2566  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2567  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /m(ks,1,i,j)
2568 
2569  do k = ks+1, ke-1
2570  z(k,i,j) = (v(k,i,j) &
2571  -( m(k,2,i,j) * z(k-1,i ,j ) &
2572  + m(k,4,i,j) * z(k ,i-1,j ) &
2573  + m(k,5,i,j) * z(k ,i+1,j ) &
2574  + m(k,6,i,j) * z(k ,i ,j-1) &
2575  + m(k,7,i,j) * z(k ,i ,j+1) &
2576  + m(k,8,i,j) * z(k-1,i-1,j ) &
2577  + m(k,9,i,j) * z(k-1,i+1,j ) &
2578  + m(k,10,i,j)* z(k-1,i ,j-1) &
2579  + m(k,11,i,j)* z(k-1,i ,j+1) &
2580  + m(k,12,i,j)* z(k+1,i-1,j ) &
2581  + m(k,13,i,j)* z(k+1,i+1,j ) &
2582  + m(k,14,i,j)* z(k+1,i ,j-1) &
2583  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/m(k,1,i,j)
2584  enddo
2585 
2586  z(ke,i,j) = ( v(ke,i,j) &
2587  -( m(ke,2,i,j) * z(ke-1,i ,j ) &
2588  + m(ke,4,i,j) * z(ke ,i-1,j ) &
2589  + m(ke,5,i,j) * z(ke ,i+1,j ) &
2590  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2591  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2592  + m(ke,8,i,j) * z(ke-1,i-1,j ) &
2593  + m(ke,9,i,j) * z(ke-1,i+1,j ) &
2594  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2595  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/m(ke,1,i,j)
2596 
2597  end if
2598 #elif COLORING == 1
2599  !$omp parallel do
2600  !$acc parallel vector_length(1)
2601  do j = js+1, je, 2
2602  !$acc loop seq
2603  do i = is, ie
2604  z(ks,i,j) = (v(ks,i,j) &
2605  -( m(ks,4,i,j) * z(ks ,i-1,j ) &
2606  + m(ks,6,i,j) * z(ks ,i ,j-1) &
2607  + m(ks,7,i,j) * z(ks ,i ,j+1) &
2608  + m(ks,12,i,j)* z(ks+1,i-1,j ) &
2609  + m(ks,14,i,j)* z(ks+1,i ,j-1) &
2610  + m(ks,15,i,j)* z(ks+1,i ,j+1) ) ) /m(ks,1,i,j)
2611 
2612  do k = ks+1, ke-1
2613  z(k,i,j) = (v(k,i,j) &
2614  -( m(k,2,i,j) * z(k-1,i ,j ) &
2615  + m(k,4,i,j) * z(k ,i-1,j ) &
2616  + m(k,6,i,j) * z(k ,i ,j-1) &
2617  + m(k,7,i,j) * z(k ,i ,j+1) &
2618  + m(k,8,i,j) * z(k-1,i-1,j ) &
2619  + m(k,10,i,j)* z(k-1,i ,j-1) &
2620  + m(k,11,i,j)* z(k-1,i ,j+1) &
2621  + m(k,12,i,j)* z(k+1,i-1,j ) &
2622  + m(k,14,i,j)* z(k+1,i ,j-1) &
2623  + m(k,15,i,j)* z(k+1,i ,j+1) ) )/m(k,1,i,j)
2624  enddo
2625 
2626  z(ke,i,j) = ( v(ke,i,j) &
2627  -( m(ke,2,i,j) * z(ke-1,i ,j ) &
2628  + m(ke,4,i,j) * z(ke ,i-1,j ) &
2629  + m(ke,6,i,j) * z(ke ,i ,j-1) &
2630  + m(ke,7,i,j) * z(ke ,i ,j+1) &
2631  + m(ke,8,i,j) * z(ke-1,i-1,j ) &
2632  + m(ke,10,i,j)* z(ke-1,i ,j-1) &
2633  + m(ke,11,i,j)* z(ke-1,i ,j+1) ) )/m(ke,1,i,j)
2634 #endif
2635  end do
2636  end do
2637  !$acc end parallel
2638 #endif
2639 
2640  !$acc end data
2641 
2642  return
2643  end subroutine gs
2644 
2645  subroutine mul_matrix( KA,KS,KE, &
2646  IA,IS,IE, &
2647  JA,JS,JE, &
2648  V, M, C)
2649  implicit none
2650  integer, intent(in) :: KA, KS, KE
2651  integer, intent(in) :: IA, IS, IE
2652  integer, intent(in) :: JA, JS, JE
2653  real(RP), intent(out) :: V(KA,IA,JA)
2654  real(RP), intent(in) :: M(KA,15,IA,JA)
2655  real(RP), intent(in) :: C(KA,IA,JA)
2656 
2657  integer :: k, i, j
2658 
2659  !$omp parallel do private(i,j,k)
2660  !$acc kernels copyin(M,C) copyout(V)
2661  do j = js, je
2662  do i = is, ie
2663  v(ks,i,j) = m(ks,1,i,j) * c(ks ,i ,j ) &
2664  + m(ks,3,i,j) * c(ks+1,i ,j ) &
2665  + m(ks,4,i,j) * c(ks ,i-1,j ) &
2666  + m(ks,5,i,j) * c(ks ,i+1,j ) &
2667  + m(ks,6,i,j) * c(ks ,i ,j-1) &
2668  + m(ks,7,i,j) * c(ks ,i ,j+1) &
2669  + m(ks,12,i,j)* c(ks+1,i-1,j ) &
2670  + m(ks,13,i,j)* c(ks+1,i+1,j ) &
2671  + m(ks,14,i,j)* c(ks+1,i ,j-1) &
2672  + m(ks,15,i,j)* c(ks+1,i ,j+1)
2673  do k = ks+1, ke-1
2674  v(k,i,j) = m(k,1,i,j) * c(k ,i ,j ) &
2675  + m(k,2,i,j) * c(k-1,i ,j ) &
2676  + m(k,3,i,j) * c(k+1,i ,j ) &
2677  + m(k,4,i,j) * c(k ,i-1,j ) &
2678  + m(k,5,i,j) * c(k ,i+1,j ) &
2679  + m(k,6,i,j) * c(k ,i ,j-1) &
2680  + m(k,7,i,j) * c(k ,i ,j+1) &
2681  + m(k,8,i,j) * c(k-1,i-1,j ) &
2682  + m(k,9,i,j) * c(k-1,i+1,j ) &
2683  + m(k,10,i,j)* c(k-1,i ,j-1) &
2684  + m(k,11,i,j)* c(k-1,i ,j+1) &
2685  + m(k,12,i,j)* c(k+1,i-1,j ) &
2686  + m(k,13,i,j)* c(k+1,i+1,j ) &
2687  + m(k,14,i,j)* c(k+1,i ,j-1) &
2688  + m(k,15,i,j)* c(k+1,i ,j+1)
2689  enddo
2690  v(ke,i,j) = m(ke,1,i,j) * c(ke ,i ,j ) &
2691  + m(ke,2,i,j) * c(ke-1,i ,j ) &
2692  + m(ke,4,i,j) * c(ke ,i-1,j ) &
2693  + m(ke,5,i,j) * c(ke ,i+1,j ) &
2694  + m(ke,6,i,j) * c(ke ,i ,j-1) &
2695  + m(ke,7,i,j) * c(ke ,i ,j+1) &
2696  + m(ke,8,i,j) * c(ke-1,i-1,j ) &
2697  + m(ke,9,i,j) * c(ke-1,i+1,j ) &
2698  + m(ke,10,i,j)* c(ke-1,i ,j-1) &
2699  + m(ke,11,i,j)* c(ke-1,i ,j+1)
2700  enddo
2701  enddo
2702  !$acc end kernels
2703 
2704  return
2705  end subroutine mul_matrix
2706 
2707  !-----------------------------------------------------------------------------
2708  function f2h( k,i,j,p )
2709 
2710  use scale_atmos_grid_cartesc, only: &
2711  cdz => atmos_grid_cartesc_cdz
2712  use scale_atmos_grid_cartesc_metric, only: &
2713  gsqrt => atmos_grid_cartesc_metric_gsqrt
2714  use scale_atmos_grid_cartesc_index, only: &
2715  i_xyz
2716 
2717  implicit none
2718 
2719  real(RP) :: f2h
2720  integer, intent(in) :: k, i, j, p
2721 
2722  f2h = (cdz(k+p-1)*gsqrt(k+p-1,i,j,i_xyz) &
2723  / (cdz(k)*gsqrt(k,i,j,i_xyz)+cdz(k+1)*gsqrt(k+1,i,j,i_xyz)))
2724 
2725  end function
2726  !-----------------------------------------------------------------------------
2728  !-----------------------------------------------------------------------------
2729  subroutine atmos_phy_lt_neutralization_mg2001( &
2730  KA, KS, KE, & ! [IN]
2731  IA, IS, IE, & ! [IN]
2732  JA, JS, JE, & ! [IN]
2733  KIJMAX, & ! [IN]
2734  IMAX, & ! [IN]
2735  JMAX, & ! [IN]
2736  Efield, & ! [IN]
2737  E_pot, & ! [IN]
2738  DENS, & ! [IN]
2739  QCRG, & ! [IN]
2740  QHYD, & ! [IN]
2741  NUM_end, & ! [INOUT]
2742  LT_path, & ! [INOUT]
2743  fls_int_p, & ! [OUT]
2744  d_QCRG ) ! [OUT]
2745  use scale_const, only: &
2746  eps => const_eps, &
2747  epsvac => const_epsvac, &
2748  epsair => const_epsair, &
2749  large_num => const_huge
2750  use scale_atmos_grid_cartesc, only: &
2751  cz => atmos_grid_cartesc_cz, &
2752  cx => atmos_grid_cartesc_cx, &
2753  cy => atmos_grid_cartesc_cy, &
2754  cdz => atmos_grid_cartesc_cdz, &
2755  cdx => atmos_grid_cartesc_cdx, &
2756  cdy => atmos_grid_cartesc_cdy, &
2757  rfdz => atmos_grid_cartesc_rfdz, &
2758  rfdx => atmos_grid_cartesc_rfdx, &
2759  rfdy => atmos_grid_cartesc_rfdy
2760  use scale_prc_cartesc, only: &
2761  prc_has_w, &
2762  prc_has_e, &
2763  prc_has_n, &
2764  prc_has_s, &
2765  prc_w, &
2766  prc_n, &
2767  prc_e, &
2768  prc_s, &
2769  prc_nw, &
2770  prc_ne, &
2771  prc_sw, &
2772  prc_se, &
2773  prc_next, &
2774  prc_num_x, &
2775  prc_num_y
2776  use scale_prc, only: &
2777  prc_abort, &
2778  prc_ismaster, &
2779  prc_nprocs, &
2780  prc_masterrank, &
2781  prc_mpibarrier, &
2782  prc_myrank
2783  use scale_comm_cartesc, only: &
2784  comm_datatype, &
2785  comm_world, &
2786  comm_vars8, &
2787  comm_wait, &
2788  comm_bcast
2789  use scale_random, only: &
2790  random_uniform
2791  implicit none
2792 
2793  integer, intent(in) :: KA, KS, KE
2794  integer, intent(in) :: IA, IS, IE
2795  integer, intent(in) :: JA, JS, JE
2796  integer, intent(in) :: KIJMAX, IMAX, JMAX
2797  real(RP), intent(in) :: Efield (KA,IA,JA,4) !-- Electric field [V/m] (4->|E|)
2798  real(RP), intent(in) :: E_pot (KA,IA,JA) !-- Electric potential [V]
2799  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
2800  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
2801  real(RP), intent(in) :: QHYD (KA,IA,JA) !-- Total hydrometeor mixing ratio
2802  real(RP), intent(inout) :: NUM_end (KA,IA,JA,3) !-- Number of Lightning stop grid
2803  real(RP), intent(inout) :: LT_path (KA,IA,JA) !-- Number of path
2804  real(RP), intent(out) :: fls_int_p(KA,IA,JA) !-- Flash initiation point (0->no flash, 1->flash start at the point)
2805  real(RP), intent(out) :: d_QCRG (KA,IA,JA) !-- Charge density [nC/m3]
2806 
2807  real(RP) :: E_det, E_x, E_y, E_z, E_max
2808  real(RP) :: dqrho_pls(KA,IA,JA), dqrho_mns(KA,IA,JA)
2809  real(RP) :: L_path(KA,IA,JA)
2810  real(RP) :: sumdqrho_pls, sumdqrho_mns
2811  real(RP) :: Ncrit, dqrho_cor
2812  real(RP) :: randnum(1)
2813  real(RP) :: rprod1, rprod2, rbuf, rbuf2, phi_end(2)
2814  integer :: Eovid(4,KIJMAX)
2815  integer :: Npls, Nmns, Ntot
2816  integer :: count1(PRC_nprocs)
2817  integer :: countindx(PRC_nprocs+1)
2818  integer :: own_prc_total !--- number of grid whose |E| > Eint-dEint for each process
2819  integer :: iprod(2), ibuf(6), status(MPI_STATUS_SIZE)
2820  integer :: init_point, rank_initpoint, grid_initpoint
2821  integer :: current_prc !--- 1->lightning flash path is include, 0-> nopath in the node
2822  integer :: comm_flag, comm_target, stop_flag, corr_flag(2)
2823  integer :: end_grid(4), wild_grid(6)
2824  real(RP) :: pm
2825  integer :: k, i, j, ipp, iq, direction, ierr
2826  integer :: k1, i1, j1, k2, i2, j2, sign_flash
2827  integer :: wild_flag, count_path
2828  integer :: flg_num_end(KA,IA,JA,3)
2829  real(RP) :: Eint_hgt(KA,IA,JA)
2830 
2831  call prof_rapstart('LT_neut_MG2001', 1)
2832 
2833  !$omp parallel do
2834  do j = js, je
2835  do i = is, ie
2836  do k = ks, ke
2837  fls_int_p(k,i,j) = 0.0_rp
2838  end do
2839  end do
2840  end do
2841 
2842  do count_path = 1, nutr_itmax
2843 
2844 
2845  !--- search grid whose Efield is over threshold of flash inititation
2846  own_prc_total = 0
2847  do j = js, je
2848  do i = is, ie
2849  do k = ks, ke
2850  eint_hgt(k,i,j) = min( eint*dens(k,i,j)/rho0,eint )*flg_eint_hgt &
2851  + ( eint-deleint )*( 1.0_rp-flg_eint_hgt )
2852  if( flg_eint_hgt == 1.0_rp .and. eint_hgt(k,i,j) < estop ) then
2853  eint_hgt(k,i,j) = large_num !--- ignore when Eint < Estop
2854  endif
2855  e_det = efield(k,i,j,i_lt_abs) - eint_hgt(k,i,j)
2856  if( e_det > 0.0_rp ) then
2857  own_prc_total = own_prc_total + 1
2858  eovid(1,own_prc_total) = i
2859  eovid(2,own_prc_total) = j
2860  eovid(3,own_prc_total) = k
2861  eovid(4,own_prc_total) = prc_myrank
2862  endif
2863  enddo
2864  enddo
2865  enddo
2866 
2867  call mpi_allgather( own_prc_total, 1, mpi_integer, &
2868  count1, 1, mpi_integer, &
2869  comm_world, ierr )
2870  countindx(1) = 0
2871 
2872  !**** countindx(PROC_nprocs+1) -> total number of grid with |E|>E_threthold
2873  do ipp = 1, prc_nprocs
2874  countindx(ipp+1) = countindx(ipp) + count1(ipp)
2875  enddo
2876 
2877  !---- select initial point of flash by random select
2878  !**** rank_initpoint -> rank number including initpoint
2879  !**** grid_initpoint -> grid number of init point in rank (rank_initpoint)
2880  if( prc_ismaster ) then
2881  call random_uniform(randnum)
2882  init_point = int( randnum(1) * countindx(prc_nprocs+1) ) + 1
2883  do ipp = 1, prc_nprocs
2884  if ( countindx(ipp+1) >= init_point ) then
2885  rank_initpoint = ipp-1
2886  exit
2887  end if
2888  end do
2889  grid_initpoint = init_point - countindx(rank_initpoint+1)
2890  ibuf(1) = rank_initpoint
2891  ibuf(2) = grid_initpoint
2892  endif
2893 
2894  call comm_bcast( 2, ibuf )
2895 
2896  rank_initpoint = ibuf(1)
2897  grid_initpoint = ibuf(2)
2898 
2899 
2900  !--- Propagate lightning
2901 
2902  !$omp parallel do
2903  do j = js, je
2904  do i = is, ie
2905  do k = ks, ke
2906  l_path(k,i,j) = 0.0_rp !--- +-1 -> path already passed, +-2 -> path calculate current step
2907  flg_num_end(k,i,j,:) = 0
2908  end do
2909  end do
2910  end do
2911  wild_flag = 0
2912  end_grid(:) = 0
2913  comm_target = -999
2914  do iq = 1, 2 !--- loop for direction (1-> parallel, 2-> anti-parallel)
2915  if( iq == 1 ) then
2916  pm = 1.0_rp
2917  elseif( iq == 2 ) then
2918  pm = - 1.0_rp
2919  endif
2920  stop_flag = 0
2921  current_prc = rank_initpoint
2922 
2923  !---- determine initiation point
2924  if( rank_initpoint == prc_myrank ) then
2925  i = eovid(1,grid_initpoint)
2926  j = eovid(2,grid_initpoint)
2927  k = eovid(3,grid_initpoint)
2928  l_path(k,i,j) = pm
2929  if( iq == 1 ) fls_int_p(k,i,j) = fls_int_p(k,i,j) + 1.0_rp
2930  sign_flash = 2 + int( sign( 1.0_rp, qcrg(k,i,j) ) )
2931  else
2932  i = 0
2933  j = 0
2934  k = 0
2935  endif
2936 
2937  loop_path : do ipp = 1, kijmaxg !--- loop for path
2938  comm_flag = 0
2939  !---- calculate path of lightning
2940  if( current_prc == prc_myrank ) then
2941 
2942  !--- determine direction of path
2943  e_x = abs( efield(k,i,j,i_lt_x) )/cdx(i)
2944  e_y = abs( efield(k,i,j,i_lt_y) )/cdy(j)
2945  e_z = abs( efield(k,i,j,i_lt_z) )/cdz(k)
2946  e_det = max( e_x,e_y,e_z )
2947 
2948  i1 = i
2949  j1 = j
2950  k1 = k
2951  if( e_det == e_x ) then
2952  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_x) ) * pm )
2953  i1 = i+direction
2954  elseif( e_det == e_y ) then
2955  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_y) ) * pm )
2956  j1 = j+direction
2957  elseif( e_det == e_z ) then
2958  direction = int( sign( 1.0_rp,efield(k,i,j,i_lt_z) ) * pm )
2959  k1 = k+direction
2960  endif
2961 
2962 
2963  if( efield(k1,i1,j1,i_lt_abs) <= estop ) then
2964  !--- stop if |E| < Estop
2965  phi_end(iq) = qcrg(k,i,j)
2966  wild_flag = 1
2967  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
2968  flg_num_end(k,i,j,sign_flash) = 1
2969  l_path(k,i,j) = pm
2970  stop_flag = 1
2971  end_grid(1) = i
2972  end_grid(2) = j
2973  end_grid(3) = k
2974  end_grid(4) = prc_myrank
2975  if( qhyd(k,i,j) < eps ) then
2976  corr_flag(iq) = 0
2977  else
2978  corr_flag(iq) = 1
2979  endif
2980  elseif( efield(k1,i1,j1,i_lt_abs) > estop ) then
2981  !--- propagate lightning path
2982  l_path(k, i ,j ) = pm
2983  l_path(k1,i1,j1) = pm
2984  endif
2985 
2986  !--- check whether lightning path reach top or bottom layer
2987  if( k1 == ke ) then
2988  !--- reach at top
2989  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
2990  flg_num_end(k,i,j,sign_flash) = 1
2991  l_path(k,i,j) = pm
2992  stop_flag = 1
2993  end_grid(1) = i
2994  end_grid(2) = j
2995  end_grid(3) = k
2996  end_grid(4) = prc_myrank
2997  if( qhyd(k,i,j) < eps ) then
2998  corr_flag(iq) = 0
2999  else
3000  corr_flag(iq) = 1
3001  endif
3002  elseif( cz(k1) < zcg ) then
3003  !--- reach at groud
3004  num_end(k,i,j,2) = num_end(k,i,j,2) + 1.0_rp
3005  flg_num_end(k,i,j,2) = 1
3006  l_path(k,i,j) = pm
3007  stop_flag = 1
3008  end_grid(1) = i
3009  end_grid(2) = j
3010  end_grid(3) = k
3011  end_grid(4) = prc_myrank
3012  if( qhyd(k,i,j) < eps ) then
3013  corr_flag(iq) = 0
3014  else
3015  corr_flag(iq) = 1
3016  endif
3017  endif
3018 
3019  if( stop_flag /= 1 ) then
3020  !--- check whether lightning path reachs boundary in i-direction
3021  if( i1 == is-1 .and. .not. prc_has_w ) then
3022  !--- reach west boundary of global domain(stop)
3023  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
3024  flg_num_end(k,i,j,sign_flash) = 1
3025  l_path(k,i,j) = pm
3026  l_path(k1,i1,j1) = pm
3027  stop_flag = 1
3028  end_grid(1) = i
3029  end_grid(2) = j
3030  end_grid(3) = k
3031  end_grid(4) = prc_myrank
3032  if( qhyd(k,i,j) < eps ) then
3033  corr_flag(iq) = 0
3034  else
3035  corr_flag(iq) = 1
3036  endif
3037  elseif( i1 == is-1 .and. prc_has_w ) then
3038  !--- reach west boundary of local domain(propagate)
3039  l_path(k,i,j) = pm
3040  l_path(k1,i1,j1) = pm
3041  comm_flag = 1
3042  comm_target = prc_next(prc_w)
3043  i1 = ie
3044  elseif( i1 == ie+1 .and. .not. prc_has_e ) then
3045  !--- reach east boundary of global domain(stop)
3046  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
3047  flg_num_end(k,i,j,sign_flash) = 1
3048  l_path(k,i,j) = pm
3049  l_path(k1,i1,j1) = pm
3050  stop_flag = 1
3051  end_grid(1) = i
3052  end_grid(2) = j
3053  end_grid(3) = k
3054  end_grid(4) = prc_myrank
3055  if( qhyd(k,i,j) < eps ) then
3056  corr_flag(iq) = 0
3057  else
3058  corr_flag(iq) = 1
3059  endif
3060  elseif( i1 == ie+1 .and. prc_has_e ) then
3061  !--- reach east boundary of local domain(propagate)
3062  l_path(k,i,j) = pm
3063  l_path(k1,i1,j1) = pm
3064  comm_flag = 1
3065  comm_target = prc_next(prc_e)
3066  i1 = is
3067  endif
3068 
3069  !--- check whether lightning path reachs boundary in i-direction
3070  if( j1 == js-1 .and. .not. prc_has_s ) then
3071  !--- reach south boundary of global domain(stop)
3072  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
3073  flg_num_end(k,i,j,sign_flash) = 1
3074  l_path(k,i,j) = pm
3075  l_path(k1,i1,j1) = pm
3076  stop_flag = 1
3077  end_grid(1) = i
3078  end_grid(2) = j
3079  end_grid(3) = k
3080  end_grid(4) = prc_myrank
3081  if( qhyd(k,i,j) < eps ) then
3082  corr_flag(iq) = 0
3083  else
3084  corr_flag(iq) = 1
3085  endif
3086  elseif( j1 == js-1 .and. prc_has_s ) then
3087  !--- reach south boundary of local domain(propagate)
3088  l_path(k,i,j) = pm
3089  l_path(k1,i1,j1) = pm
3090  comm_flag = 1
3091  comm_target = prc_next(prc_s)
3092  j1 = je
3093  elseif( j1 == je+1 .and. .not. prc_has_n ) then
3094  !--- reach north boundary of global domain(stop)
3095  num_end(k,i,j,sign_flash) = num_end(k,i,j,sign_flash) + 1.0_rp
3096  flg_num_end(k,i,j,sign_flash) = 1
3097  l_path(k,i,j) = pm
3098  l_path(k1,i1,j1) = pm
3099  comm_flag = 1
3100  stop_flag = 1
3101  end_grid(1) = i
3102  end_grid(2) = j
3103  end_grid(3) = k
3104  end_grid(4) = prc_myrank
3105  if( qhyd(k,i,j) < eps ) then
3106  corr_flag(iq) = 0
3107  else
3108  corr_flag(iq) = 1
3109  endif
3110  elseif( j1 == je+1 .and. prc_has_n ) then
3111  !--- reach north boundary of local domain(propagate)
3112  l_path(k,i,j) = pm
3113  l_path(k1,i1,j1) = pm
3114  comm_flag = 1
3115  comm_target = prc_next(prc_n)
3116  j1 = js
3117  endif
3118  endif
3119 
3120  endif
3121 
3122  !---- determine stop or not
3123  call mpi_bcast(stop_flag, 1, mpi_integer, current_prc, comm_world, ierr)
3124 
3125  if( stop_flag == 1 ) then
3126  !---- send flag wildfire
3127  call mpi_bcast(wild_flag, 1, mpi_integer, current_prc, comm_world, ierr)
3128 
3129  if ( prc_myrank == current_prc ) then
3130  ibuf(1:4) = end_grid(:)
3131  ibuf(5:6) = corr_flag(:)
3132  end if
3133  call mpi_bcast(ibuf, 6, mpi_integer, current_prc, comm_world, ierr)
3134  end_grid(:) = ibuf(1:4)
3135  corr_flag(:) = ibuf(5:6)
3136 
3137  stop_flag = 0
3138  !---- If lightning path reaches end_grid stop
3139  exit loop_path
3140  endif
3141 
3142  !---- determine wether process change occurs or not
3143  call mpi_bcast(comm_flag, 1, mpi_integer, current_prc, comm_world, ierr)
3144 
3145  !--- process change occurs
3146  if( comm_flag == 1 ) then
3147  call mpi_bcast(comm_target, 1, mpi_integer, current_prc, comm_world, ierr)
3148 
3149  if ( prc_myrank == comm_target ) then
3150  call mpi_recv(ibuf, 6, mpi_integer, current_prc, 1, comm_world, status, ierr)
3151  k1 = ibuf(1)
3152  i1 = ibuf(2)
3153  j1 = ibuf(3)
3154  corr_flag(:) = ibuf(4:5)
3155  sign_flash = ibuf(6)
3156  end if
3157 
3158  if ( prc_myrank == current_prc ) then
3159  ibuf(1) = k1
3160  ibuf(2) = i1
3161  ibuf(3) = j1
3162  ibuf(4:5) = corr_flag(:)
3163  ibuf(6) = sign_flash
3164  call mpi_send(ibuf, 6, mpi_integer, comm_target, 1, comm_world, ierr)
3165  end if
3166 
3167  !--- change current proc
3168  current_prc = comm_target
3169  endif
3170 
3171  if( current_prc == prc_myrank ) then
3172  k = k1
3173  j = j1
3174  i = i1
3175  endif
3176 
3177  enddo loop_path
3178 
3179 ! call PRC_MPIbarrier
3180 
3181  !---- Wildfire method
3182  if( wild_flag == 1 .and. end_grid(4) == prc_myrank ) then
3183  i1 = end_grid(1)
3184  j1 = end_grid(2)
3185  k1 = end_grid(3)
3186  i2 = end_grid(1)+1
3187  j2 = end_grid(2)
3188  k2 = end_grid(3)
3189  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3190  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3191  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3192  endif
3193  i1 = end_grid(1)
3194  j1 = end_grid(2)
3195  k1 = end_grid(3)
3196  i2 = end_grid(1)-1
3197  j2 = end_grid(2)
3198  k2 = end_grid(3)
3199  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3200  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3201  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3202  endif
3203  i1 = end_grid(1)
3204  j1 = end_grid(2)
3205  k1 = end_grid(3)
3206  i2 = end_grid(1)
3207  j2 = end_grid(2)+1
3208  k2 = end_grid(3)
3209  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3210  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3211  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3212  endif
3213  i1 = end_grid(1)
3214  j1 = end_grid(2)
3215  k1 = end_grid(3)
3216  i2 = end_grid(1)
3217  j2 = end_grid(2)-1
3218  k2 = end_grid(3)
3219  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3220  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3221  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3222  endif
3223  i1 = end_grid(1)
3224  j1 = end_grid(2)
3225  k1 = end_grid(3)
3226  i2 = end_grid(1)
3227  j2 = end_grid(2)
3228  k2 = end_grid(3)+1
3229  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3230  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3231  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3232  endif
3233  i1 = end_grid(1)
3234  j1 = end_grid(2)
3235  k1 = end_grid(3)
3236  i2 = end_grid(1)
3237  j2 = end_grid(2)
3238  k2 = end_grid(3)-1
3239  if( abs( qcrg(k2,i2,j2) ) > qrho_chan/dens(k2,i2,j2) .and. &
3240  abs( e_pot(k2,i2,j2) ) > abs( phi_end(iq) ) ) then
3241  l_path( k2,i2,j2 ) = l_path( k1,i1,j1 )
3242  endif
3243  endif
3244 
3245 ! call PRC_MPIbarrier
3246 
3247  enddo !--- loop for direction
3248 
3249 ! call COMM_vars8( L_path,1 )
3250 ! call COMM_wait ( L_path,1 )
3251 
3252  !---- Neutralization
3253  npls = 0
3254  nmns = 0
3255  sumdqrho_pls = 0.0_rp
3256  sumdqrho_mns = 0.0_rp
3257  !$omp parallel do reduction(+:Npls,Nmns,sumdqrho_pls,sumdqrho_mns)
3258  do j = js, je
3259  do i = is, ie
3260  do k = ks, ke
3261 
3262  !---- whether the grid is on lightning path or not
3263  if( l_path(k,i,j) /= 0.0_rp .and. abs( qcrg(k,i,j) ) > qrho_chan ) then
3264  if ( qcrg(k,i,j) >= 0.0_rp ) then
3265  npls = npls + 1
3266  dqrho_pls(k,i,j) = fp * ( abs( qcrg(k,i,j) )-qrho_neut )
3267  sumdqrho_pls = sumdqrho_pls &
3268  + fp * ( abs( qcrg(k,i,j) )-qrho_neut )
3269  dqrho_mns(k,i,j) = 0.0_rp
3270  else
3271  nmns = nmns + 1
3272  dqrho_mns(k,i,j) = fp * ( abs( qcrg(k,i,j) )-qrho_neut )
3273  sumdqrho_mns = sumdqrho_mns &
3274  + fp * ( abs( qcrg(k,i,j) )-qrho_neut )
3275  dqrho_pls(k,i,j) = 0.0_rp
3276  endif
3277  else
3278  dqrho_pls(k,i,j) = 0.0_rp
3279  dqrho_mns(k,i,j) = 0.0_rp
3280  endif
3281 
3282  enddo
3283  enddo
3284  enddo
3285  iprod(1) = npls
3286  iprod(2) = nmns
3287  call mpi_allreduce(iprod, ibuf, 2, mpi_integer, mpi_sum, comm_world, ierr)
3288  npls = ibuf(1)
3289  nmns = ibuf(2)
3290 
3291  ntot = npls + nmns
3292 
3293  if ( ntot > 0 ) exit
3294 
3295  end do
3296 
3297  if ( count_path > nutr_itmax ) then
3298  log_info("ATMOS_PHY_LT_Efield",*) "Reach limit iteration for searching flash path", count_path, npls, nmns, current_prc
3299  endif
3300 
3301  if( corr_flag(1) == 1 .and. corr_flag(2) == 1 ) then
3302  dqrho_cor = sumdqrho_pls - sumdqrho_mns
3303  rprod2 = dqrho_cor
3304  call mpi_allreduce(rprod2, rbuf, 1, comm_datatype, mpi_sum, comm_world, ierr)
3305  dqrho_cor = rbuf
3306 
3307  dqrho_cor = dqrho_cor / ntot !---- Eq.(2) of MacGorman et al. (2001)
3308  else
3309  dqrho_cor = 0.0_rp !---- No correlation if flash to ground and out of cloud
3310  endif
3311 
3312  !$omp parallel do
3313  do j = js, je
3314  do i = is, ie
3315  do k = ks, ke
3316  d_qcrg(k,i,j) = 0.0_rp
3317  end do
3318  end do
3319  end do
3320  rprod1 = 0.0_rp
3321  rprod2 = 0.0_rp
3322  !--- pseud-2D experiment for x-direction
3323  if( prc_num_x == 1 .and. imax == 2 ) then
3324 
3325  !$omp parallel do
3326  do j = js, je
3327  do k = ks, ke
3328  if( l_path(k,is,j) /= 0.0_rp ) then
3329  l_path(k,ie,j) = l_path(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
3330  if( dqrho_pls(k,is,j) /= 0.0_rp ) then
3331  d_qcrg(k,is,j) = - dqrho_pls(k,is,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3332  d_qcrg(k,ie,j) = d_qcrg(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
3333  elseif( dqrho_mns(k,is,j) /= 0.0_rp ) then
3334  d_qcrg(k,is,j) = dqrho_mns(k,is,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3335  d_qcrg(k,ie,j) = d_qcrg(k,is,j) !--- copy IS value to IE (which is regarded as same grid in psued-2D)
3336  endif
3337  elseif( l_path(k,ie,j) /= 0.0_rp ) then
3338  l_path(k,is,j) = l_path(k,ie,j) !--- copy IE value to IS (which is regarded as same grid in psued-2D)
3339  if( dqrho_pls(k,ie,j) /= 0.0_rp ) then
3340  d_qcrg(k,ie,j) = - dqrho_pls(k,ie,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3341  d_qcrg(k,is,j) = d_qcrg(k,ie,j) !--- copy IS value to IS (which is regarded as same grid in psued-2D)
3342  elseif( dqrho_mns(k,ie,j) /= 0.0_rp ) then
3343  d_qcrg(k,ie,j) = dqrho_mns(k,ie,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3344  d_qcrg(k,is,j) = d_qcrg(k,ie,j) !--- copy IS value to IS (which is regarded as same grid in psued-2D)
3345  endif
3346  endif
3347 
3348  do iq = 1, 3
3349  if( flg_num_end(k,is,j,iq) == 1 ) then
3350  num_end(k,ie,j,iq) = num_end(k,ie,j,iq) + 1.0_rp
3351  endif
3352  if( flg_num_end(k,ie,j,iq) == 1 ) then
3353  num_end(k,is,j,iq) = num_end(k,is,j,iq) + 1.0_rp
3354  endif
3355  enddo
3356 
3357  enddo
3358  enddo
3359 
3360  !--- pseud-2D experiment for y-direction
3361  elseif( prc_num_y == 1 .and. jmax == 2 ) then
3362 
3363  !$omp parallel do
3364  do i = is, ie
3365  do k = ks, ke
3366  if( l_path(k,i,js) /= 0.0_rp ) then
3367  l_path(k,i,je) = l_path(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
3368  if( dqrho_pls(k,i,js) /= 0.0_rp ) then
3369  d_qcrg(k,i,js) = - dqrho_pls(k,i,js) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3370  d_qcrg(k,i,je) = d_qcrg(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
3371  elseif( dqrho_mns(k,i,js) /= 0.0_rp ) then
3372  d_qcrg(k,i,js) = dqrho_mns(k,i,js) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3373  d_qcrg(k,i,je) = d_qcrg(k,i,js) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
3374  endif
3375  elseif( l_path(k,i,je) /= 0.0_rp ) then
3376  l_path(k,i,js) = l_path(k,i,je) !--- copy JE value to JS (which is regarded as same grid in psued-2D)
3377  if( dqrho_pls(k,i,je) /= 0.0_rp ) then
3378  d_qcrg(k,i,je) = - dqrho_pls(k,i,je) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3379  d_qcrg(k,i,js) = d_qcrg(k,i,je) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
3380  elseif( dqrho_mns(k,i,je) /= 0.0_rp ) then
3381  d_qcrg(k,i,je) = dqrho_mns(k,i,je) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3382  d_qcrg(k,i,js) = d_qcrg(k,i,je) !--- copy JS value to JE (which is regarded as same grid in psued-2D)
3383  endif
3384  endif
3385 
3386  do iq = 1, 3
3387  if( flg_num_end(k,i,js,iq) == 1 ) then
3388  num_end(k,i,je,iq) = num_end(k,i,je,iq) + 1.0_rp
3389  endif
3390  if( flg_num_end(k,i,je,iq) == 1 ) then
3391  num_end(k,i,js,iq) = num_end(k,i,js,iq) + 1.0_rp
3392  endif
3393  enddo
3394 
3395  enddo
3396  enddo
3397 
3398  else
3399 
3400  !$omp parallel do
3401  do j = js, je
3402  do i = is, ie
3403  do k = ks, ke
3404  if( l_path(k,i,j) /= 0.0_rp .and. dqrho_pls(k,i,j) /= 0.0_rp ) then
3405  d_qcrg(k,i,j) = - dqrho_pls(k,i,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3406  elseif( l_path(k,i,j) /= 0.0_rp .and. dqrho_mns(k,i,j) /= 0.0_rp ) then
3407  d_qcrg(k,i,j) = dqrho_mns(k,i,j) - dqrho_cor !--- Eq.(3) MacGorman et al. (2001)
3408  endif
3409  enddo
3410  enddo
3411  enddo
3412 
3413  endif
3414 
3415  !$omp parallel do
3416  do j = js, je
3417  do i = is, ie
3418  do k = ks, ke
3419  lt_path(k,i,j) = lt_path(k,i,j) + l_path(k,i,j)
3420  enddo
3421  enddo
3422  enddo
3423 
3424  call prof_rapend('LT_neut_MG2001', 1)
3425 
3426  return
3427 
3428  end subroutine atmos_phy_lt_neutralization_mg2001
3429  !-----------------------------------------------------------------------------
3431  !-----------------------------------------------------------------------------
3432  subroutine atmos_phy_lt_neutralization_f2013( &
3433  KA, KS, KE, & ! [IN]
3434  IA, IS, IE, & ! [IN]
3435  JA, JS, JE, & ! [IN]
3436  KIJMAX, & ! [IN]
3437  Efield, & ! [IN]
3438  E_pot, & ! [IN]
3439  DENS, & ! [IN]
3440  QCRG, & ! [IN]
3441  QHYD, & ! [IN]
3442  NUM_end, & ! [INOUT]
3443  LT_path, & ! [INOUT]
3444  fls_int_p, & ! [OUT]
3445  d_QCRG, & ! [OUT]
3446  B_OUT ) ! [OUT]
3447  use scale_const, only: &
3448  eps => const_eps, &
3449  epsvac => const_epsvac, &
3450  epsair => const_epsair
3451  use scale_atmos_grid_cartesc, only: &
3452  cz => atmos_grid_cartesc_cz, &
3453  cx => atmos_grid_cartesc_cx, &
3454  cy => atmos_grid_cartesc_cy, &
3455  cdz => atmos_grid_cartesc_cdz, &
3456  cdx => atmos_grid_cartesc_cdx, &
3457  cdy => atmos_grid_cartesc_cdy, &
3458  rfdz => atmos_grid_cartesc_rfdz, &
3459  rfdx => atmos_grid_cartesc_rfdx, &
3460  rfdy => atmos_grid_cartesc_rfdy
3461  use scale_prc_cartesc, only: &
3462  prc_has_w, &
3463  prc_has_e, &
3464  prc_has_n, &
3465  prc_has_s, &
3466  prc_w, &
3467  prc_n, &
3468  prc_e, &
3469  prc_s, &
3470  prc_nw, &
3471  prc_ne, &
3472  prc_sw, &
3473  prc_se, &
3474  prc_next, &
3475  prc_num_x,&
3476  prc_num_y
3477  use scale_prc, only: &
3478  prc_abort, &
3479  prc_ismaster, &
3480  prc_nprocs, &
3481  prc_masterrank, &
3482  prc_mpibarrier, &
3483  prc_myrank
3484  use scale_comm_cartesc, only: &
3485  comm_datatype, &
3486  comm_world
3487 ! use scale_random, only: &
3488 ! RANDOM_reset, &
3489 ! RANDOM_uniform
3490  implicit none
3491 
3492  integer, intent(in) :: KA, KS, KE
3493  integer, intent(in) :: IA, IS, IE
3494  integer, intent(in) :: JA, JS, JE
3495  integer, intent(in) :: KIJMAX
3496  real(RP), intent(in) :: Efield (KA,IA,JA,4) !-- Electric field [V/m] (4->|E|)
3497  real(RP), intent(in) :: E_pot (KA,IA,JA) !-- Electric potential [V]
3498  real(RP), intent(in) :: QCRG (KA,IA,JA) !-- Charge density [nC/m3]
3499  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
3500  real(RP), intent(in) :: QHYD (KA,IA,JA) !-- Total hydrometeor mixing ratio
3501  real(RP), intent(inout) :: NUM_end (KA,IA,JA,3) !-- Number of Lightning stop grid
3502  real(RP), intent(inout) :: LT_path (KA,IA,JA) !-- Number of path
3503  real(RP), intent(out) :: fls_int_p(KA,IA,JA) !-- Flash initiation point (0->no flash, 1->flash start at the point)
3504  real(RP), intent(out) :: d_QCRG (KA,IA,JA) !-- Charge density [nC/m3]
3505  real(RP), intent(out) :: B_OUT (IA,JA) !-- B value for output
3506 
3507  real(RP), parameter :: q_thre = 0.1_rp ! threshold of discharge zone (Fierro et al. 2013) [nC/m3]
3508 
3509  real(RP) :: B(IA,JA) !--- B in Fierro et al. 2013
3510  integer :: C(IA,JA) !--- flg for cylinder (in cylinder -> C=1)
3511  real(RP) :: Q_d, Fpls, Fmns
3512  real(RP) :: Spls, Smns
3513  real(RP) :: Spls_g, Smns_g
3514  real(RP) :: distance
3515  real(RP) :: Edif(KS:KE), Edif_max, abs_qcrg_max
3516  real(RP) :: sw
3517 
3518  !integer, allocatable :: proc_num(:), proc_numg(:)
3519  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)
3520  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)
3521  real(RP) :: exce_grid(2,KIJMAX)
3522  integer :: count1(PRC_nprocs)
3523  integer :: countindx(PRC_nprocs+1)
3524  integer :: num_own !--- number of column whose |E| > Eint-dEint for each process
3525  integer :: num_total !--- total number of column whose |E| > Eint-dEint
3526  real(RP) :: rbuf1, rbuf2
3527  integer :: k, i, j, ipp, iq, ierr
3528 
3529  call prof_rapstart('LT_neut_F2013', 1)
3530 
3531  !$acc data copyin(Efield,E_pot,QCRG,DENS,QHYD,CX,CY) &
3532  !$acc copy(NUM_end,LT_path) &
3533  !$acc copyout(fls_int_p,d_QCRG,B_OUT) &
3534  !$acc create(B,C,Edif,exce_grid)
3535 
3536  !$omp parallel do
3537  !$acc kernels
3538  do j = js, je
3539  do i = is, ie
3540  do k = ks, ie
3541  do iq = 1, 3
3542  num_end(k,i,j,iq) = 0.0_rp
3543  end do
3544  end do
3545  end do
3546  end do
3547  !$acc end kernels
3548 
3549  !--- search grid whose Efield is over threshold of flash inititation
3550  num_own = 0
3551  !$acc kernels
3552  !$acc loop collapse(2) independent
3553  do j = js, je
3554  do i = is, ie
3555  edif_max = -1.0_rp
3556  !$acc loop reduction(max:Edif_max)
3557  do k = ks, ke
3558  edif(k) = efield(k,i,j,i_lt_abs) - ( eint-deleint )
3559  edif_max = max( edif_max, edif(k) )
3560  enddo
3561  if ( edif_max > 0.0_rp ) then
3562  !$acc atomic
3563  num_own = num_own + 1
3564  !$acc end atomic
3565  exce_grid(1,num_own) = cx(i)
3566  exce_grid(2,num_own) = cy(j)
3567  do k = ks, ke
3568  fls_int_p(k,i,j) = 0.5_rp + sign( 0.5_rp, edif(k)-eps )
3569  enddo
3570  else
3571  do k = ks, ke
3572  fls_int_p(k,i,j) = 0.0_rp
3573  enddo
3574  endif
3575  enddo
3576  enddo
3577  !$acc end kernels
3578 
3579  !############## calculation by CPU ################
3580  !$acc update host(exce_grid)
3581  !**** proc_num(0~) -> process number of each grid with |E|> E_threthold (local)
3582 ! allocate(proc_num(own_prc_total))
3583  allocate(e_exce_x(num_own))
3584  allocate(e_exce_y(num_own))
3585 
3586 ! proc_num(:) = PRC_myrank
3587  e_exce_x(1:num_own) = exce_grid(1,1:num_own)
3588  e_exce_y(1:num_own) = exce_grid(2,1:num_own)
3589 
3590  call mpi_allgather( num_own, 1, mpi_integer, &
3591  count1(:), 1, mpi_integer, &
3592  comm_world, ierr )
3593 
3594  countindx(1) = 0
3595  do ipp = 1, prc_nprocs
3596  countindx(ipp+1) = countindx(ipp) + count1(ipp)
3597  enddo
3598 
3599  !---- Create global version of proc_num(proc_numg)
3600  !**** countindx(PROC_nprocs) -> total number of grid with |E|>E_threthold
3601  !**** proc_numg(0~) -> process number of each grid with |E|> E_threthold (global)
3602  num_total = countindx(prc_nprocs+1)
3603  allocate(e_exce_x_g(num_total))
3604  allocate(e_exce_y_g(num_total))
3605 
3606  call mpi_allgatherv( e_exce_x, num_own, comm_datatype, &
3607  e_exce_x_g, count1, countindx, comm_datatype, &
3608  comm_world, ierr )
3609 
3610  call mpi_allgatherv( e_exce_y, num_own, comm_datatype, &
3611  e_exce_y_g, count1, countindx, comm_datatype, &
3612  comm_world, ierr )
3613  !$acc data copyin(E_exce_x_g,E_exce_y_g)
3614  !############## calculation by CPU ################
3615 
3616  !$acc kernels
3617  c(:,:) = 0
3618  do ipp = 1, num_total
3619  do j = js, je
3620  do i = is, ie
3621  distance = sqrt( ( cx(i)-e_exce_x_g(ipp) )**2 &
3622  + ( cy(j)-e_exce_y_g(ipp) )**2 )
3623  if( distance <= r_neut ) then
3624  c(i,j) = 1
3625  endif
3626  enddo
3627  enddo
3628  enddo
3629  !$acc end kernels
3630 
3631  !$acc end data
3632 
3633  spls = 0.0_rp
3634  smns = 0.0_rp
3635  !$omp parallel do reduction(+:Spls,Smns) &
3636  !$omp private(abs_qcrg_max,sw)
3637  !$acc kernels
3638  !$acc loop collapse(2) reduction(+:Spls,Smns)
3639  do j = js, je
3640  do i = is, ie
3641 
3642  b(i,j) = 0.0_rp
3643  if ( c(i,j) == 1 ) then
3644  abs_qcrg_max = 0.0_rp
3645  !$acc loop reduction(max:abs_qcrg_max)
3646  do k = ks, ke
3647  abs_qcrg_max = max( abs_qcrg_max, abs( qcrg(k,i,j) ) )
3648  end do
3649  if ( abs_qcrg_max >= q_thre ) then
3650  b(i,j) = 1.0_rp
3651  do k = ks, ke
3652  sw = 0.5_rp + sign( 0.5_rp, qcrg(k,i,j) )
3653  spls = spls + qcrg(k,i,j) * sw
3654  smns = smns - qcrg(k,i,j) * ( 1.0_rp - sw )
3655  enddo
3656  end if
3657  end if
3658 
3659  enddo
3660  enddo
3661  !$acc end kernels
3662 
3663  rbuf1 = spls
3664  call mpi_allreduce( rbuf1, rbuf2, 1, comm_datatype, &
3665  mpi_sum, comm_world, ierr )
3666  spls_g = rbuf2
3667 
3668  rbuf1 = smns
3669  call mpi_allreduce( rbuf1, rbuf2, 1, comm_datatype, &
3670  mpi_sum, comm_world, ierr )
3671  smns_g = rbuf2
3672 
3673  if( max( spls_g,smns_g )*0.3_rp < min( spls_g,smns_g ) ) then
3674  q_d = 0.3_rp * max( spls_g,smns_g )
3675  else
3676  q_d = min( spls_g,smns_g )
3677  endif
3678 
3679  if( spls_g /= 0.0_rp ) then
3680  fpls = q_d / spls_g
3681  else
3682  fpls = 0.0_rp
3683  endif
3684 
3685  if( smns_g /= 0.0_rp ) then
3686  fmns = q_d / smns_g
3687  else
3688  fmns = 0.0_rp
3689  endif
3690 
3691  !---- select initial point of flash by random select
3692  !$omp parallel do
3693  !$acc kernels
3694  do j = js, je
3695  do i = is, ie
3696  do k = ks, ke
3697  if( qcrg(k,i,j) > q_thre ) then
3698  d_qcrg(k,i,j) = - fpls * ( qcrg(k,i,j) - q_thre ) * b(i,j)
3699  lt_path(k,i,j) = lt_path(k,i,j) + b(i,j)
3700  elseif( qcrg(k,i,j) < -q_thre ) then
3701  d_qcrg(k,i,j) = + fmns * ( - qcrg(k,i,j) - q_thre ) * b(i,j)
3702  lt_path(k,i,j) = lt_path(k,i,j) - b(i,j)
3703  else
3704  d_qcrg(k,i,j) = 0.0_rp
3705  endif
3706  enddo
3707  enddo
3708  enddo
3709  !$acc end kernels
3710 
3711  !$acc kernels
3712  do j = js, je
3713  do i = is, ie
3714  b_out(i,j) = b(i,j)
3715  enddo
3716  enddo
3717  !$acc end kernels
3718 
3719  deallocate(e_exce_x)
3720  deallocate(e_exce_y)
3721  deallocate(e_exce_x_g)
3722  deallocate(e_exce_y_g)
3723 
3724  !$acc end data
3725 
3726  call prof_rapend('LT_neut_F2013', 1)
3727 
3728  return
3729  end subroutine atmos_phy_lt_neutralization_f2013
3730  !-----------------------------------------------------------------------------
3732  !-----------------------------------------------------------------------------
3733  subroutine atmos_phy_lt_judge_abse( &
3734  KA, KS, KE, & ! [IN]
3735  IA, IS, IE, & ! [IN]
3736  JA, JS, JE, & ! [IN]
3737  DENS, & ! [IN]
3738  Efield, & ! [IN]
3739  Emax, & ! [OUT]
3740  flg_lt_neut ) ! [OUT]
3741  use scale_comm_cartesc, only: &
3742  comm_datatype, &
3743  comm_world
3744  use scale_const, only: &
3745  large_num => const_huge
3746  implicit none
3747 
3748  integer, intent(in) :: KA, KS, KE
3749  integer, intent(in) :: IA, IS, IE
3750  integer, intent(in) :: JA, JS, JE
3751  real(RP), intent(in) :: DENS (KA,IA,JA) !-- Total density [kg/m3]
3752  real(RP), intent(in) :: Efield (KA,IA,JA) !-- Absolute value of Electric field [V/m]
3753  real(RP), intent(out) :: Emax !-- Maximum value of Electric field [V/m]
3754  integer, intent(out) :: flg_lt_neut !-- flg 1-> neutralization, 0->no neutralization
3755 
3756  real(RP) :: Eint_hgt(KA,IA,JA)
3757  real(RP) :: E_det, rbuf, rprod1
3758  integer :: own_prc_total, iprod1, buf
3759  integer :: k, i, j, ierr
3760 
3761  !$acc data &
3762  !$acc copyin(DENS,Efield) &
3763  !$acc create(Eint_hgt)
3764 
3765  !--- search grid whose Efield is over threshold of flash inititation
3766  own_prc_total = 0
3767  if( flg_eint_hgt == 1.0_rp ) then
3768  !$omp parallel do &
3769  !$omp reduction(+:own_prc_total) &
3770  !$omp private(E_det)
3771  !$acc kernels
3772  !$acc loop collapse(3) reduction(+:own_prc_total)
3773  do j = js, je
3774  do i = is, ie
3775  do k = ks, ke
3776  eint_hgt(k,i,j) = min( eint*dens(k,i,j)/rho0,eint )
3777  if( eint_hgt(k,i,j) < estop ) then
3778  eint_hgt(k,i,j) = large_num !--- ignore when Eint < Estop
3779  endif
3780  e_det = efield(k,i,j) - eint_hgt(k,i,j)
3781  if( e_det > 0.0_rp ) then
3782  own_prc_total = own_prc_total + 1
3783  endif
3784  enddo
3785  enddo
3786  enddo
3787  !$acc end kernels
3788  else
3789  !$omp parallel do &
3790  !$omp reduction(+:own_prc_total) &
3791  !$omp private(E_det)
3792  !$acc kernels
3793  !$acc loop collapse(3) reduction(+:own_prc_total)
3794  do j = js, je
3795  do i = is, ie
3796  do k = ks, ke
3797  e_det = efield(k,i,j) - ( eint-deleint )
3798  if( e_det > 0.0_rp ) then
3799  own_prc_total = own_prc_total + 1
3800  endif
3801  enddo
3802  enddo
3803  enddo
3804  !$acc end kernels
3805  endif
3806 
3807  !--- Add number of grids, whose |E| are over threshold of flash initiaion, for all process
3808  iprod1 = own_prc_total
3809  call mpi_allreduce(iprod1, buf, 1, mpi_integer, mpi_sum, comm_world, ierr)
3810  rprod1 = 0.0_rp
3811  !$acc kernels
3812  !$acc loop collapse(3) reduction(max:rprod1)
3813  do j = js, je
3814  do i = is, ie
3815  do k = ks, ke
3816  rprod1 = max(rprod1,efield(k,i,j))
3817  enddo
3818  enddo
3819  enddo
3820  !$acc end kernels
3821  call mpi_allreduce(rprod1, rbuf, 1, comm_datatype, mpi_max, comm_world, ierr)
3822  !--- exit when no grid point with |E| over threshold of flash initiation exist
3823  emax = rbuf
3824  if( buf == 0 ) then
3825  flg_lt_neut = 0
3826  else
3827  flg_lt_neut = 1
3828  if( emax < eint .and. emax >= eint-deleint .and. flg_eint_hgt == 0.0_rp ) then
3829  flg_lt_neut = 2
3830  endif
3831  endif
3832 
3833  !$acc end data
3834 
3835  end subroutine atmos_phy_lt_judge_abse
3836  !-----------------------------------------------------------------------------
3838  !-----------------------------------------------------------------------------
3839  subroutine atmos_phy_lt_sato2019_select_dqcrg_from_lut( &
3840  KA, KS, KE, & ! [IN]
3841  IA, IS, IE, & ! [IN]
3842  JA, JS, JE, & ! [IN]
3843  NLIQ, & ! [IN]
3844  TEMP, & ! [IN]
3845  DENS, & ! [IN]
3846  QLIQ, & ! [IN]
3847  dqcrg, & ! [OUT]
3848  beta_crg ) ! [OUT]
3849  use scale_const, only: &
3850  t00 => const_tem00
3851  implicit none
3852 
3853  integer, intent(in) :: ka, ks, ke
3854  integer, intent(in) :: ia, is, ie
3855  integer, intent(in) :: ja, js, je
3856  integer, intent(in) :: nliq
3857  real(rp), intent(in) :: temp(ka,ia,ja)
3858  real(rp), intent(in) :: dens(ka,ia,ja)
3859  real(rp), intent(in) :: qliq(ka,ia,ja,nliq)
3860  real(rp), intent(out) :: dqcrg(ka,ia,ja)
3861  real(rp), intent(out) :: beta_crg(ka,ia,ja)
3862 
3863  integer :: i, j, k, pp, qq, iq
3864  integer :: grid1, grid2
3865  real(rp) :: cwc
3866  real(rp) :: diffx, diffy, tmp
3867 
3868  !$acc data copyin(TEMP,DENS,QLIQ) copyout(dqcrg,beta_crg)
3869 
3870  !$omp parallel do &
3871  !$omp private(cwc,diffx,diffy,tmp,grid1,grid2)
3872  !$acc kernels
3873  do j = js, je
3874  do i = is, ie
3875  do k = ks, ke
3876  if( temp(k,i,j) <= t00 .and. temp(k,i,j) >= tcrglimit ) then
3877  cwc = 0.0_rp
3878  !$acc loop seq
3879  do iq = 1, nliq
3880  cwc = cwc + qliq(k,i,j,iq) * dens(k,i,j) * 1.0e+3_rp ![g/m3]
3881  enddo
3882  diffx = abs( grid_lut_t(1)-temp(k,i,j) )
3883  grid1 = 1
3884  !$acc loop seq
3885  do pp = 2, nxlut_lt
3886  tmp = abs( grid_lut_t(pp)-temp(k,i,j) )
3887  if ( tmp < diffx ) then
3888  diffx = tmp
3889  grid1 = pp
3890  end if
3891  enddo
3892  diffy = abs( grid_lut_l(1)-cwc )
3893  grid2 = 1
3894  !$acc loop seq
3895  do qq = 2, nylut_lt
3896  tmp = abs( grid_lut_l(qq)-cwc )
3897  if ( tmp < diffy ) then
3898  diffy = tmp
3899  grid2 = qq
3900  end if
3901  enddo
3902  dqcrg(k,i,j) = dq_chrg( grid1, grid2 ) &
3903  *( 0.5_rp + sign( 0.5_rp,cwc-1.0e-2_rp ) ) !--- no charge separation when cwc < 0.01 [g/m3]
3904  else
3905  dqcrg(k,i,j) = 0.0_rp
3906  endif
3907  if( temp(k,i,j) >= -30.0_rp+t00 ) then
3908  beta_crg(k,i,j) = 1.0_rp
3909  elseif( temp(k,i,j) < -30.0_rp+t00 .and. temp(k,i,j) >= -43.0_rp+t00 ) then
3910  beta_crg(k,i,j) = 1.0_rp - ( ( temp(k,i,j)-t00+30.0_rp )/13.0_rp )**2
3911 ! elseif( TEMP(k,i,j) < -43.0_RP+T00 ) then
3912  else
3913  beta_crg(k,i,j) = 0.0_rp
3914  endif
3915  enddo
3916  enddo
3917  enddo
3918  !$acc end kernels
3919 
3920  !$acc end data
3921 
3922  return
3923  end subroutine atmos_phy_lt_sato2019_select_dqcrg_from_lut
3924  !-----------------------------------------------------------------------------
3925 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:100
scale_comm_cartesc::comm_datatype
integer, public comm_datatype
datatype of variable
Definition: scale_comm_cartesC.F90:105
scale_atmos_grid_cartesc_index::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:101
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_prc_cartesc::prc_n
integer, parameter, public prc_n
[node direction] north
Definition: scale_prc_cartesC.F90:34
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:181
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:43
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:1326
scale_prc_cartesc::prc_sw
integer, parameter, public prc_sw
[node direction] southwest
Definition: scale_prc_cartesC.F90:39
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:68
scale_prc_cartesc::prc_next
integer, dimension(8), public prc_next
node ID of 8 neighbour process
Definition: scale_prc_cartesC.F90:46
scale_prc_cartesc::prc_has_s
logical, public prc_has_s
Definition: scale_prc_cartesC.F90:51
scale_const::const_epsvac
real(rp), parameter, public const_epsvac
parts par million
Definition: scale_const.F90:103
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:520
scale_atmos_phy_lt_sato2019::atmos_phy_lt_sato2019_finalize
subroutine, public atmos_phy_lt_sato2019_finalize
finalize
Definition: scale_atmos_phy_lt_sato2019.F90:490
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:3741
scale_atmos_phy_lt_sato2019::ilu_decomp
subroutine ilu_decomp(KA, KS, KE, IA, IS, IE, JA, JS, JE, M, diag)
Definition: scale_atmos_phy_lt_sato2019.F90:1982
scale_atmos_grid_cartesc_index::i_uv
integer, public i_uv
Definition: scale_atmos_grid_cartesC_index.F90:102
scale_atmos_phy_lt_sato2019::back_sub_ilu
subroutine back_sub_ilu(KA, KS, KE, IA, IS, IE, JA, JS, JE, Z, M, V, diag)
Definition: scale_atmos_phy_lt_sato2019.F90:2080
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:66
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:174
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:62
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:91
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:373
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:45
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:40
scale_prc_cartesc::prc_has_n
logical, public prc_has_n
Definition: scale_prc_cartesC.F90:49
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:91
scale_prc_cartesc::prc_has_e
logical, public prc_has_e
Definition: scale_prc_cartesC.F90:50
scale_prc
module PROCESS
Definition: scale_prc.F90:11
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, B_OUT)
calculate charge neutralization Fierro et al. (2013)
Definition: scale_atmos_phy_lt_sato2019.F90:3447
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:35
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:99
scale_prc::prc_mpibarrier
subroutine, public prc_mpibarrier
Barrier MPI.
Definition: scale_prc.F90:828
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:95
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:67
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_prc::prc_masterrank
integer, parameter, public prc_masterrank
master process in each communicator
Definition: scale_prc.F90:67
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:38
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_prc_cartesc::prc_w
integer, parameter, public prc_w
[node direction] west
Definition: scale_prc_cartesC.F90:33
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:37
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:93
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:40
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
scale_comm_cartesc::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm_cartesC.F90:106
scale_prc_cartesc::prc_s
integer, parameter, public prc_s
[node direction] south
Definition: scale_prc_cartesC.F90:36
scale_prc::prc_nprocs
integer, public prc_nprocs
myrank in local communicator
Definition: scale_prc.F90:90
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:99
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:37
scale_prc_cartesc::prc_num_y
integer, public prc_num_y
y length of 2D processor topology
Definition: scale_prc_cartesC.F90:43
scale_prc_cartesc::prc_e
integer, parameter, public prc_e
[node direction] east
Definition: scale_prc_cartesC.F90:35
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:57
scale_const::const_epsair
real(rp), public const_epsair
parts par million
Definition: scale_const.F90:104
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:46
scale_atmos_phy_lt_sato2019
module atmosphere / physics / lightninh / SATO2019
Definition: scale_atmos_phy_lt_sato2019.F90:27
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:685
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:246
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:41
scale_prc_cartesc::prc_num_x
integer, public prc_num_x
x length of 2D processor topology
Definition: scale_prc_cartesC.F90:42
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:61
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:39
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:3849
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:41
scale_atmos_grid_cartesc_index::i_uvz
integer, public i_uvz
Definition: scale_atmos_grid_cartesC_index.F90:97
scale_atmos_phy_lt_sato2019::gs_ilu
subroutine gs_ilu(KA, KS, KE, IA, IS, IE, JA, JS, JE, Z, M, V, diag)
Definition: scale_atmos_phy_lt_sato2019.F90:2150
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:2745
scale_prc_cartesc::prc_ne
integer, parameter, public prc_ne
[node direction] northeast
Definition: scale_prc_cartesC.F90:38
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:69
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:1500
scale_prc_cartesc::prc_has_w
logical, public prc_has_w
Definition: scale_prc_cartesC.F90:48
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
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:56
scale_prc::prc_ismaster
logical, public prc_ismaster
master process in local communicator?
Definition: scale_prc.F90:92
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:92
Generated by   doxygen 1.8.17