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