SCALE-RM
scale_atmos_phy_mp_common.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
11 !-------------------------------------------------------------------------------
12 #include "scalelib.h"
13 module scale_atmos_phy_mp_common
14  !-----------------------------------------------------------------------------
15  !
16  !++ used modules
17  !
18  use scale_precision
19  use scale_io
20  use scale_prof
21  use scale_atmos_grid_cartesc_index
22  use scale_tracer
23  !-----------------------------------------------------------------------------
24  implicit none
25  private
26  !-----------------------------------------------------------------------------
27  !
28  !++ Public procedure
29  !
30  public :: atmos_phy_mp_negative_fixer
31  public :: atmos_phy_mp_saturation_adjustment
32  public :: atmos_phy_mp_precipitation_upwind
33  public :: atmos_phy_mp_precipitation_semilag
34  public :: atmos_phy_mp_precipitation_momentum
35 
36  interface atmos_phy_mp_saturation_adjustment
37  module procedure atmos_phy_mp_saturation_adjustment_3d
38  end interface atmos_phy_mp_saturation_adjustment
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  !-----------------------------------------------------------------------------
53 contains
54 
55  !-----------------------------------------------------------------------------
59  !-----------------------------------------------------------------------------
60  subroutine atmos_phy_mp_negative_fixer( &
61  KA, KS, KE, IA, IS, IE, JA, JS, JE, QLA, QIA, &
62  limit_negative, &
63  DENS, TEMP, &
64  CVtot, CPtot, &
65  QV, QTRC, &
66  RHOH, &
67  DENS_diff, &
68  ENGI_diff )
69  use scale_prc, only: &
70  prc_abort
71  use scale_const, only: &
72  cvdry => const_cvdry, &
73  cpdry => const_cpdry
74  use scale_atmos_hydrometeor, only: &
75  atmos_hydrometeor_lhv, &
76  atmos_hydrometeor_lhs, &
77  cv_vapor, &
78  cp_vapor, &
79  cv_water, &
80  cp_water, &
81  cv_ice, &
82  cp_ice
83  implicit none
84  integer, intent(in) :: ka, ks, ke
85  integer, intent(in) :: ia, is, ie
86  integer, intent(in) :: ja, js, je
87  integer, intent(in) :: qla, qia
88 
89  real(rp), intent(in) :: limit_negative
90 
91  real(rp), intent(inout) :: dens (ka,ia,ja)
92  real(rp), intent(inout) :: temp (ka,ia,ja)
93  real(rp), intent(inout) :: cvtot(ka,ia,ja)
94  real(rp), intent(inout) :: cptot(ka,ia,ja)
95  real(rp), intent(inout) :: qv (ka,ia,ja)
96  real(rp), intent(inout) :: qtrc (ka,ia,ja,qla+qia)
97 
98  real(rp), intent(out), optional :: rhoh (ka,ia,ja)
99  real(rp), intent(out), optional :: dens_diff(ka,ia,ja)
100  real(rp), intent(out), optional :: engi_diff(ka,ia,ja)
101 
102  real(rp) :: lhv(ka)
103  real(rp) :: lhs(ka)
104  real(rp) :: eng (ka)
105  real(rp) :: eng0 (ka)
106  real(rp) :: dens0(ka)
107  real(rp) :: diffq(ka)
108  real(rp) :: diffq_min
109  real(rp) :: work
110 
111  logical :: rhoh_out
112  logical :: dens_out
113  logical :: engi_out
114 
115  integer :: k, i, j, iq
116  !---------------------------------------------------------------------------
117 
118  call prof_rapstart('MP_filter', 3)
119 
120  rhoh_out = present( rhoh )
121  dens_out = present( dens_diff )
122  engi_out = present( engi_diff )
123 
124  diffq_min = 0.0_rp
125 
126  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
127  !$omp reduction(min:diffq_min) &
128  !$omp private(i,j,k,iq, &
129  !$omp LHV,LHS,eng,eng0,dens0,diffq,work) &
130  !$omp shared(KA,KS,KE,IS,IE,JS,JE,QLA,QIA, &
131  !$omp IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID, &
132  !$omp CVdry,CPdry,CV_VAPOR,CP_VAPOR,CV_WATER,CP_WATER,CV_ICE,CP_ICE, &
133  !$omp DENS,TEMP,CVtot,CPtot,QV,QTRC, &
134  !$omp RHOH,DENS_diff,ENGI_diff,rhoh_out,dens_out,engi_out,limit_negative)
135  !$acc kernels copy(DENS, TEMP, CVtot, CPtot, QV, QTRC) copyout(RHOH, DENS_diff, ENGI_diff)
136  !$acc loop reduction(min:diffq_min)
137  do j = js, je
138  !$acc loop private(LHV, LHS, eng, eng0, dens0, diffq) reduction(min:diffq_min)
139  do i = is, ie
140 
141  diffq(:) = 0.0_rp
142 
143  do k = ks, ke
144  eng(k) = cvtot(k,i,j) * temp(k,i,j)
145  end do
146 
147  if ( rhoh_out ) then
148  do k = ks, ke
149  eng0(k) = eng(k)
150  end do
151  end if
152 
153  if ( qla > 0 ) then
154  call atmos_hydrometeor_lhv( ka, ks, ke, temp(:,i,j), lhv(:) )
155 
156  do iq = 1, qla
157  do k = ks, ke
158  work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to liq)
159  if ( work > 0.0_rp ) then
160  eng(k) = eng(k) + work * lhv(k)
161  diffq(k) = diffq(k) - work
162  cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_water - cv_vapor )
163  cptot(k,i,j) = cptot(k,i,j) + work * ( cp_water - cp_vapor )
164  ! remove negative value of hydrometeors (mass)
165  qtrc(k,i,j,iq) = 0.0_rp
166  end if
167  enddo
168  enddo
169  end if
170 
171  if ( qia > 0 ) then
172  call atmos_hydrometeor_lhs( ka, ks, ke, temp(:,i,j), lhs(:) )
173 
174  do iq = qla+1, qla+qia
175  do k = ks, ke
176  work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to ice)
177  if ( work > 0.0_rp ) then
178  eng(k) = eng(k) + work * lhs(k)
179  diffq(k) = diffq(k) - work
180  cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_ice - cv_vapor )
181  cptot(k,i,j) = cptot(k,i,j) + work * ( cp_ice - cp_vapor )
182  ! remove negative value of hydrometeors (mass)
183  qtrc(k,i,j,iq) = 0.0_rp
184  end if
185  enddo
186  enddo
187  end if
188 
189 
190  if ( abs(limit_negative) > 0.0_rp ) then
191  !$acc loop reduction(min:diffq_min)
192  do k = ks, ke
193  if ( diffq(k) < - abs(limit_negative) ) then
194  diffq_min = min( diffq_min, diffq(k) )
195 #ifndef _OPENACC
196  log_error("ATMOS_PHY_MP_negative_fixer",*) 'large negative is found'
197  log_error_cont(*) 'value = ', diffq(k), ' at (', k, ',', i, ',', j, ')'
198 #endif
199  end if
200  end do
201  end if
202 
203  do k = ks, ke
204  if ( diffq(k) < 0.0_rp ) then
205  ! Compensate for the lack of hydrometeors by the water vapor
206  qv(k,i,j) = qv(k,i,j) + diffq(k)
207  end if
208  end do
209 
210  if ( rhoh_out ) then
211  do k = ks, ke
212  rhoh(k,i,j) = ( eng(k) - eng0(k) ) * dens(k,i,j)
213  end do
214  end if
215 
216 
217  if ( dens_out ) then
218  do k = ks, ke
219  dens0(k) = dens(k,i,j)
220  end do
221  end if
222  if ( engi_out ) then
223  do k = ks, ke
224  eng0(k) = eng(k) * dens(k,i,j)
225  end do
226  end if
227 
228  ! fix negative QV
229  do k = ks, ke
230  work = - min( qv(k,i,j), 0.0_rp )
231  if ( work > 0.0_rp ) then
232  qv(k,i,j) = 0.0_rp
233  dens(k,i,j) = dens(k,i,j) * ( 1.0_rp + work ) ! not to change mass of dry air
234  cvtot(k,i,j) = ( cvtot(k,i,j) + work * cv_vapor ) / ( 1.0_rp + work )
235  cptot(k,i,j) = ( cptot(k,i,j) + work * cp_vapor ) / ( 1.0_rp + work )
236  eng(k) = ( eng(k) + work * cv_vapor * temp(k,i,j) ) / ( 1.0_rp + work )
237  end if
238  enddo
239 
240  do k = ks, ke
241  ! update
242  temp(k,i,j) = eng(k) / cvtot(k,i,j)
243  end do
244 
245  if ( dens_out ) then
246  do k = ks, ke
247  dens_diff(k,i,j) = dens(k,i,j) - dens0(k)
248  end do
249  end if
250  if ( engi_out ) then
251  do k = ks, ke
252  engi_diff(k,i,j) = eng(k) * dens(k,i,j) - eng0(k)
253  end do
254  end if
255 
256  enddo
257  enddo
258  !$acc end kernels
259 
260 
261  if ( abs(limit_negative) > 0.0_rp &
262  .AND. abs(limit_negative) < abs(diffq_min) ) then
263  log_error_cont(*) 'maximum negative hydrometeor ', diffq_min, ' < ', - abs(limit_negative)
264  call prc_abort
265  endif
266 
267  call prof_rapend('MP_filter', 3)
268 
269  return
270  end subroutine atmos_phy_mp_negative_fixer
271 
272  !-----------------------------------------------------------------------------
274  !-----------------------------------------------------------------------------
275  subroutine atmos_phy_mp_saturation_adjustment_3d( &
276  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
277  DENS, &
278  flag_liquid, &
279  TEMP, &
280  QV, QC, QI, &
281  CPtot, CVtot, &
282  RHOE_d )
283  use scale_prc, only: &
284  prc_abort
285  use scale_atmos_hydrometeor, only: &
286  cp_vapor, &
287  cp_water, &
288  cp_ice, &
289  cv_vapor, &
290  cv_water, &
291  cv_ice, &
292  lhv, &
293  lhf
294  use scale_atmos_saturation, only: &
295  atmos_saturation_moist_conversion_dens_all, &
296  atmos_saturation_moist_conversion_dens_liq
297  implicit none
298 
299  integer, intent(in) :: KA, KS, KE
300  integer, intent(in) :: IA, IS, IE
301  integer, intent(in) :: JA, JS, JE
302 
303  real(RP), intent(in) :: DENS(KA,IA,JA)
304  logical , intent(in) :: flag_liquid
305 
306  real(RP), intent(inout) :: TEMP (KA,IA,JA)
307  real(RP), intent(inout) :: QV (KA,IA,JA)
308  real(RP), intent(inout) :: QC (KA,IA,JA)
309  real(RP), intent(inout) :: QI (KA,IA,JA)
310  real(RP), intent(inout) :: CPtot(KA,IA,JA)
311  real(RP), intent(inout) :: CVtot(KA,IA,JA)
312 
313  real(RP), intent(out) :: RHOE_d(KA,IA,JA)
314 
315  ! working
316  real(RP) :: QV1
317  real(RP) :: QC1
318  real(RP) :: QI1
319 
320  real(RP) :: Emoist ! moist internal energy
321 
322  logical :: converged, error
323 
324  integer :: k, i, j
325  !---------------------------------------------------------------------------
326 
327  call prof_rapstart('MP_Saturation_adjustment', 2)
328 
329  error = .false.
330 
331  if ( flag_liquid ) then ! warm rain
332 
333  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
334  !$omp shared (IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID) &
335  !$omp shared (KS,KE,IS,IE,JS,JE, &
336  !$omp CP_VAPOR,CP_WATER,CV_VAPOR,CV_WATER,LHV,LHF, &
337  !$omp DENS,QV,QC,TEMP,CPtot,CVtot,RHOE_d,error) &
338  !$omp private(i,j,k, &
339  !$omp QV1,QC1,Emoist,converged)
340  !$acc kernels copyin(DENS) copyout(RHOE_d) copy(TEMP, QV, QC, CPtot, CVtot)
341  !$acc loop collapse(3) independent &
342  !$acc reduction(.or.:error) private(qv1, qc1, emoist, converged)
343  do j = js, je
344  do i = is, ie
345  do k = ks, ke
346 
347  qv1 = qv(k,i,j)
348  qc1 = qc(k,i,j)
349 
350  emoist = temp(k,i,j) * cvtot(k,i,j) &
351  + qv1 * lhv
352 
353  call atmos_saturation_moist_conversion_dens_liq( dens(k,i,j), emoist, & ! [IN]
354  temp(k,i,j), qv1, qc1, & ! [INOUT]
355  cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
356  converged ) ! [OUT]
357 
358  if ( .NOT. converged ) then
359  error = .true.
360 #ifndef _OPENACC
361  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j, dens(k,i,j), emoist, temp(k,i,j), qv(k,i,j), qc(k,i,j), cptot(k,i,j), cvtot(k,i,j)
362  exit
363 #endif
364  endif
365 
366  rhoe_d(k,i,j) = - lhv * ( qv1 - qv(k,i,j) ) * dens(k,i,j)
367 
368  qv(k,i,j) = qv1
369  qc(k,i,j) = qc1
370 
371  end do
372  end do
373  end do
374  !$acc end kernels
375 
376  if ( error ) then
377 #ifndef _OPENACC
378  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged!'
379 #endif
380  call prc_abort
381  end if
382  else ! cold rain
383 
384  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
385  !$omp shared (IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID) &
386  !$omp shared (KS,KE,IS,IE,JS,JE, &
387  !$omp CP_VAPOR,CP_WATER,CP_ICE,CV_VAPOR,CV_WATER,CV_ICE,LHV,LHF, &
388  !$omp DENS,QV,QC,QI,TEMP,CPtot,CVtot,RHOE_d,error) &
389  !$omp private(i,j,k, &
390  !$omp QV1,QC1,QI1,Emoist,converged)
391  !$acc kernels copyin(DENS) copyout(RHOE_d) copy(TEMP, QV, QC, QI, CPtot, CVtot)
392  !$acc loop collapse(3) independent &
393  !$acc reduction(.or.:error) private(qv1, qc1, qi1, emoist, converged)
394  do j = js, je
395  do i = is, ie
396  do k = ks, ke
397  qv1 = qv(k,i,j)
398  qc1 = qc(k,i,j)
399  qi1 = qi(k,i,j)
400 
401  emoist = temp(k,i,j) * cvtot(k,i,j) &
402  + qv1 * lhv &
403  - qi1 * lhf
404 
405  call atmos_saturation_moist_conversion_dens_all( dens(k,i,j), emoist, & ! [IN]
406  temp(k,i,j), qv1, qc1, qi1, & ! [INOUT]
407  cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
408  converged ) ! [OUT]
409 
410  if ( .NOT. converged ) then
411  error = .true.
412 #ifndef _OPENACC
413  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j, dens(k,i,j), emoist, temp(k,i,j), qv(k,i,j), qc(k,i,j), qi(k,i,j), cptot(k,i,j), cvtot(k,i,j)
414  exit
415 #endif
416  endif
417 
418  rhoe_d(k,i,j) = ( - lhv * ( qv1 - qv(k,i,j) ) &
419  + lhf * ( qi1 - qi(k,i,j) ) ) * dens(k,i,j)
420 
421  qv(k,i,j) = qv1
422  qc(k,i,j) = qc1
423  qi(k,i,j) = qi1
424 
425  end do
426  end do
427  end do
428  !$acc end kernels
429 
430  if ( error ) then
431 #ifdef _OPENACC
432  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged!'
433 #endif
434  call prc_abort
435  end if
436  endif
437 
438  call prof_rapend ('MP_Saturation_adjustment', 2)
439 
440  return
441  end subroutine atmos_phy_mp_saturation_adjustment_3d
442 
443  !-----------------------------------------------------------------------------
444 !OCL SERIAL
445  subroutine atmos_phy_mp_precipitation_upwind( &
446  KA, KS, KE, QHA, QLA, QIA, &
447  TEMP, vterm, FDZ, RCDZ, dt, &
448  i, j, &
449  DENS, RHOQ, CPtot, CVtot, RHOE, &
450  mflx, sflx, esflx )
451  !$acc routine vector
452  use scale_const, only: &
453  grav => const_grav
454  use scale_atmos_hydrometeor, only: &
455  cp_water, &
456  cp_ice, &
457  cv_water, &
458  cv_ice
459  implicit none
460  integer, intent(in), value :: ka, ks, ke
461  integer, intent(in), value :: qha, qla, qia
462 
463  real(rp), intent(in) :: temp (ka)
464  real(rp), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
465  real(rp), intent(in) :: fdz (ka)
466  real(rp), intent(in) :: rcdz (ka)
467  real(rp), intent(in) :: dt
468  integer, intent(in) :: i, j ! for debug
469 
470  real(rp), intent(inout) :: dens (ka)
471  real(rp), intent(inout) :: rhoq (ka,qha)
472  real(rp), intent(inout) :: cptot(ka)
473  real(rp), intent(inout) :: cvtot(ka)
474  real(rp), intent(inout) :: rhoe (ka)
475 
476  real(rp), intent(out) :: mflx (ka)
477  real(rp), intent(out) :: sflx (2)
478  real(rp), intent(out) :: esflx
479 
480 #ifdef _OPENACC
481  real(rp) :: qflx0
482  real(rp) :: qflx1
483  real(rp) :: eflx0
484  real(rp) :: eflx1
485  real(rp) :: rhocp_pu
486  real(rp) :: rhocv_pu
487 #define RHOCP_pu(k) RHOCP_pu
488 #define RHOCV_pu(k) RHOCV_pu
489 #else
490  real(rp) :: qflx(ka)
491  real(rp) :: eflx(ka)
492  real(rp) :: rhocp_pu(ka)
493  real(rp) :: rhocv_pu(ka)
494 #endif
495 
496  real(rp) :: ddens
497  real(rp) :: cp, cv
498 
499  integer :: k, iq
500  !---------------------------------------------------------------------------
501 
502  ! tracer/energy transport by falldown
503  ! 1st order upwind, forward euler, velocity is always negative
504 
505  mflx(:) = 0.0_rp
506  sflx(:) = 0.0_rp
507  esflx = 0.0_rp
508 
509 #ifndef _OPENACC
510  qflx(ke) = 0.0_rp
511  eflx(ke) = 0.0_rp
512 #endif
513 
514  do k = ks, ke
515 
516  rhocp_pu(k) = cptot(k) * dens(k)
517  rhocv_pu(k) = cvtot(k) * dens(k)
518 #ifndef _OPENACC
519  end do
520 #endif
521 
522  !$acc loop seq
523  do iq = 1, qha
524 
525  !--- mass flux for each tracer, upwind with vel < 0
526 #ifdef _OPENACC
527  if ( k == ks ) then
528  qflx0 = vterm(ks,iq) * rhoq(ks,iq)
529  else
530  qflx0 = 0.5_rp * ( vterm(k,iq) + vterm(k-1,iq) ) * rhoq(k,iq)
531  end if
532  if ( k == ke ) then
533  qflx1 = 0.0_rp
534  else
535  qflx1 = 0.5_rp * ( vterm(k+1,iq) + vterm(k,iq) ) * rhoq(k+1,iq)
536  end if
537 #else
538  qflx(ks-1) = vterm(ks,iq) * rhoq(ks,iq)
539  do k = ks, ke-1
540  qflx(k) = 0.5_rp * ( vterm(k+1,iq) + vterm(k,iq) ) * rhoq(k+1,iq)
541  enddo
542 #endif
543 
544  !--- update falling tracer
545 #ifdef _OPENACC
546  rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx1 - qflx0 ) * rcdz(k)
547 #else
548  do k = ks, ke
549  rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
550  enddo ! falling (water mass & number) tracer
551 #endif
552 
553  ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
554  if ( iq > qla + qia ) cycle
555 
556 #ifdef _OPENACC
557  mflx(k-1) = mflx(k-1) + qflx0
558 #else
559  do k = ks-1, ke-1
560  mflx(k) = mflx(k) + qflx(k)
561  end do
562 #endif
563 
564  if ( iq > qla ) then ! ice water
565  cp = cp_ice
566  cv = cv_ice
567 #ifdef _OPENACC
568  if ( k == ks ) sflx(2) = sflx(2) + qflx0
569 #else
570  sflx(2) = sflx(2) + qflx(ks-1)
571 #endif
572  else ! liquid water
573  cp = cp_water
574  cv = cv_water
575 #ifdef _OPENACC
576  if ( k == ks ) sflx(1) = sflx(1) + qflx0
577 #else
578  sflx(1) = sflx(1) + qflx(ks-1)
579 #endif
580  end if
581 
582  !--- update density
583 #ifdef _OPENACC
584  ddens = - ( qflx1 - qflx0 ) * rcdz(k) * dt
585 #else
586  do k = ks, ke
587  ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
588 #endif
589  rhocp_pu(k) = rhocp_pu(k) + cp * ddens
590  rhocv_pu(k) = rhocv_pu(k) + cv * ddens
591  dens(k) = dens(k) + ddens
592 
593  ! internal energy flux
594 #ifdef _OPENACC
595  eflx0 = qflx0 * temp(k ) * cv
596  eflx1 = qflx1 * temp(k+1) * cv
597  if ( k == ks ) esflx = esflx + eflx0
598 #else
599  end do
600 
601  do k = ks-1, ke-1
602  eflx(k) = qflx(k) * temp(k+1) * cv
603  end do
604  esflx = esflx + eflx(ks-1)
605 #endif
606 
607  !--- update internal energy
608 #ifdef _OPENACC
609  rhoe(k) = rhoe(k) - ( ( eflx1 - eflx0 ) & ! contribution with the transport of internal energy
610  + qflx1 * fdz(k) * grav & ! contribution with the release of potential energy
611  ) * rcdz(k) * dt
612 #else
613  do k = ks, ke
614  rhoe(k) = rhoe(k) - ( ( eflx(k) - eflx(k-1) ) & ! contribution with the transport of internal energy
615  + qflx(k) * fdz(k) * grav & ! contribution with the release of potential energy
616  ) * rcdz(k) * dt
617  end do
618 #endif
619 
620  end do
621 
622 #ifndef _OPENACC
623  do k = ks, ke
624 #endif
625  cptot(k) = rhocp_pu(k) / dens(k)
626  cvtot(k) = rhocv_pu(k) / dens(k)
627 
628  end do
629 
630  return
631  end subroutine atmos_phy_mp_precipitation_upwind
632 
633  !-----------------------------------------------------------------------------
634 !OCL SERIAL
635  subroutine atmos_phy_mp_precipitation_semilag( &
636  KA, KS, KE, QHA, QLA, QIA, &
637  TEMP, vterm, FZ, FDZ, RCDZ, dt, &
638  i, j, &
639  DENS, RHOQ, CPtot, CVtot, RHOE, &
640  mflx, sflx, esflx )
641  !$acc routine vector
642  use scale_const, only: &
643  grav => const_grav
644  use scale_atmos_hydrometeor, only: &
645  cp_water, &
646  cp_ice, &
647  cv_water, &
648  cv_ice
649  implicit none
650  integer, intent(in), value :: ka, ks, ke
651  integer, intent(in), value :: qha, qla, qia
652 
653  real(rp), intent(in) :: temp (ka)
654  real(rp), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
655  real(rp), intent(in) :: fz (ka)
656  real(rp), intent(in) :: fdz (ka)
657  real(rp), intent(in) :: rcdz (ka)
658  real(rp), intent(in) :: dt
659  integer, intent(in) :: i, j ! for debug
660 
661  real(rp), intent(inout) :: dens (ka)
662  real(rp), intent(inout) :: rhoq (ka,qha)
663  real(rp), intent(inout) :: cptot(ka)
664  real(rp), intent(inout) :: cvtot(ka)
665  real(rp), intent(inout) :: rhoe (ka)
666 
667  real(rp), intent(out) :: mflx (ka)
668  real(rp), intent(out) :: sflx (2)
669  real(rp), intent(out) :: esflx
670 
671  real(rp) :: qflx(ka)
672  real(rp) :: eflx(ka)
673  real(rp) :: rhocp(ka)
674  real(rp) :: rhocv(ka)
675  real(rp) :: ddens
676  real(rp) :: cp, cv
677 
678  real(rp) :: vtermh(ka)
679  real(rp) :: dvterm(ka)
680  real(rp) :: cdz (ka)
681  real(rp) :: rfdz2 (ka)
682  real(rp) :: dist (ka)
683  real(rp) :: z_src
684  real(rp) :: flx
685  integer :: k_src (ka)
686  integer :: k_dst
687 
688  integer :: k, iq
689  !---------------------------------------------------------------------------
690 
691  ! tracer/energy transport by falldown
692  ! velocity is always negative
693 
694  mflx(:) = 0.0_rp
695  sflx(:) = 0.0_rp
696 
697  do k = ks, ke
698  rhocp(k) = cptot(k) * dens(k)
699  rhocv(k) = cvtot(k) * dens(k)
700  end do
701 
702  do k = ks, ke
703  cdz(k) = 1.0_rp / rcdz(k)
704  end do
705  do k = ks, ke-1
706  rfdz2(k) = 1.0_rp / ( cdz(k) + cdz(k+1) )
707  end do
708 
709  !$acc loop seq
710  do iq = 1, qha
711 
712  qflx(:) = 0.0_rp
713  eflx(:) = 0.0_rp
714 
715  do k = ks, ke-1
716  vtermh(k) = ( cdz(k) * vterm(k+1,iq) + cdz(k+1) * vterm(k,iq) ) * rfdz2(k)
717  enddo
718  vtermh(ks-1) = vterm(ks,iq)
719 
720  do k = ks, ke
721  dvterm(k) = vtermh(k) - vtermh(k-1)
722  enddo
723 
724  ! Movement distance of the cell wall by the fall
725  ! the midpoint method (second-order Runge-Kutta)
726  ! dz/dt = v(z + v dt/2) ~ v(z) + v dt/2 dv/dz + 1/2 (v dt/2)^2 d^2v/dz^2
727  do k = ks, ke-1
728  dist(k) = - vtermh(k) * dt &
729  + vtermh(k) * dt**2 / 2.0_rp * ( dvterm(k+1)+dvterm(k) ) * rfdz2(k) &
730  - vtermh(k)**2 * dt**3 / 4.0_rp * ( dvterm(k+1)*rcdz(k+1) - dvterm(k)*rcdz(k) ) * rfdz2(k)
731  dist(k) = max( dist(k), 0.0_rp )
732  enddo
733  dist(ks-1) = - vtermh(ks-1) * dt &
734  + vtermh(ks-1) * dt**2 / 2.0_rp * dvterm(ks)*rcdz(ks)
735  dist(ks-1) = max( dist(ks-1), 0.0_rp )
736 
737  ! wall cannot overtake
738  !$acc loop seq
739  do k = ke-2, ks-1, -1
740  dist(k) = min( dist(k), dist(k+1) + cdz(k+1) )
741  end do
742 
743 ! LOG_INFO_CONT(*) "distance", iq
744 ! do k = KA, 1, -1
745 ! LOG_INFO_CONT('(1x,I5,3F9.3,ES15.5)') k, dist(k), vtermh(k), vterm(k,iq), RHOQ(k,iq)
746 ! enddo
747 
748  ! search number of source cell
749  do k_dst = ks-1, ke-1
750  z_src = fz(k_dst) + dist(k_dst)
751 
752  k_src(k_dst) = k_dst
753  do k = k_dst, ke-1
754  if ( z_src > fz(k ) &
755  .AND. z_src <= fz(k+1) ) then
756  k_src(k_dst) = k
757  exit
758  endif
759  enddo
760  if ( z_src > fz(ke) ) k_src(k_dst) = ke
761  enddo
762 
763 ! LOG_INFO_CONT(*) "seek", iq
764 ! do k = KA, 1, -1
765 ! LOG_INFO_CONT('(1x,2I5,2F9.3)') k, k_src(k), FZ(k), FZ(k)+dist(k)
766 ! enddo
767 
768  if ( iq > qla ) then ! ice water
769  cp = cp_ice
770  cv = cv_ice
771  else ! liquid water
772  cp = cp_water
773  cv = cv_water
774  end if
775 
776  do k_dst = ks-1, ke-1
777  do k = k_dst, k_src(k_dst)-1
778  flx = rhoq(k+1,iq) * cdz(k+1) / dt ! sum column mass rhoq*dz
779  qflx(k_dst) = qflx(k_dst) - flx
780  eflx(k_dst) = eflx(k_dst) - flx * temp(k+1) * cv ! internal energy flux
781  dist(k_dst) = dist(k_dst) - cdz(k+1) ! residual
782  enddo
783  if ( k_src(k_dst) < ke ) then
784  ! residual (simple upwind)
785  flx = rhoq(k_src(k_dst)+1,iq) * dist(k_dst) / dt
786  qflx(k_dst) = qflx(k_dst) - flx ! sum column mass rhoq*dz
787  eflx(k_dst) = eflx(k_dst) -flx * temp(k_src(k_dst)+1) * cv ! internal energy flux
788  end if
789  enddo
790 
791 ! LOG_INFO_CONT(*) "flux", iq
792 ! do k = KA, 1, -1
793 ! LOG_INFO_CONT('(1x,2I5,F9.3,2ES15.5)') k, k_src(k), dist(k), qflx(k), vtermh(k)*RHOQ(k+1,iq)
794 ! enddo
795 
796  !--- update falling tracer
797  do k = ks, ke
798  rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
799  enddo ! falling (water mass & number) tracer
800 
801 ! LOG_INFO_CONT(*) "tendency", iq
802 ! do k = KA, 1, -1
803 ! LOG_INFO_CONT('(1x,I5,ES15.5)') k, - dt * ( qflx(k) - qflx(k-1) ) * RCDZ(k)
804 ! enddo
805 
806  ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
807  if ( iq > qla + qia ) cycle
808 
809  do k = ks-1, ke-1
810  mflx(k) = mflx(k) + qflx(k)
811  end do
812 
813  if ( iq > qla ) then ! ice water
814  sflx(2) = sflx(2) + qflx(ks-1)
815  else ! liquid water
816  sflx(1) = sflx(1) + qflx(ks-1)
817  end if
818 
819  !--- update density
820  do k = ks, ke
821  ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
822  rhocp(k) = rhocp(k) + cp * ddens
823  rhocv(k) = rhocv(k) + cv * ddens
824  dens(k) = dens(k) + ddens
825  end do
826 
827  !--- update internal energy
828  do k = ks, ke
829  rhoe(k) = rhoe(k) - ( ( eflx(k) - eflx(k-1) ) & ! contribution with the transport of internal energy
830  + qflx(k) * fdz(k) * grav & ! contribution with the release of potential energy
831  ) * rcdz(k) * dt
832  end do
833  esflx = esflx + eflx(ks-1)
834 
835  end do
836 
837  do k = ks, ke
838  cptot(k) = rhocp(k) / dens(k)
839  cvtot(k) = rhocv(k) / dens(k)
840  end do
841 
842  return
843  end subroutine atmos_phy_mp_precipitation_semilag
844 
845  !-----------------------------------------------------------------------------
846 !OCL SERIAL
847  subroutine atmos_phy_mp_precipitation_momentum( &
848  KA, KS, KE, &
849  DENS, MOMZ, U, V, mflx, &
850  RCDZ, RFDZ, &
851  MOMZ_t, RHOU_t, RHOV_t )
852  !$acc routine vector
853  implicit none
854 
855  integer, intent(in), value :: ka, ks, ke
856  real(rp), intent(in) :: dens (ka)
857  real(rp), intent(in) :: momz (ka)
858  real(rp), intent(in) :: u (ka)
859  real(rp), intent(in) :: v (ka)
860  real(rp), intent(in) :: mflx (ka)
861  real(rp), intent(in) :: rcdz (ka)
862  real(rp), intent(in) :: rfdz (ka)
863  real(rp), intent(out) :: momz_t(ka)
864  real(rp), intent(out) :: rhou_t(ka)
865  real(rp), intent(out) :: rhov_t(ka)
866 
867 #ifdef _OPENACC
868  real(rp) :: flx0, flx1
869 #else
870  real(rp) :: flx(ka)
871 #endif
872 
873  integer :: k
874  !---------------------------------------------------------------------------
875 
876 #ifdef _OPENACC
877  do k = ks, ke
878 #else
879  flx(ke) = 0.0_rp
880 #endif
881 
882  !--- momentum z (half level)
883 #ifdef _OPENACC
884  if ( k < ke ) then
885  flx0 = ( mflx(k ) + mflx(k-1) ) * momz(k ) / ( dens(k+1) + dens(k ) )
886  flx1 = ( mflx(k+1) + mflx(k ) ) * momz(k+1) / ( dens(k+2) + dens(k+1) )
887  momz_t(k) = - ( flx1 - flx0 ) * rfdz(k)
888  else ! k = KE
889  momz_t(k) = 0.0_rp
890  end if
891 #else
892  do k = ks, ke-1
893  flx(k) = ( mflx(k) + mflx(k-1) ) * momz(k) / ( dens(k+1) + dens(k) )
894  enddo
895  do k = ks, ke-1
896  momz_t(k) = - ( flx(k+1) - flx(k) ) * rfdz(k)
897  enddo
898  momz_t(ke) = 0.0_rp
899 #endif
900 
901  !--- momentum x
902 #ifdef _OPENACC
903  flx0 = mflx(k-1) * u(k)
904  if ( k < ke ) then
905  flx1 = mflx(k) * u(k+1)
906  else
907  flx1 = 0.0_rp
908  end if
909  rhou_t(k) = - ( flx1 - flx0 ) * rcdz(k)
910 #else
911  do k = ks-1, ke-1
912  flx(k) = mflx(k) * u(k+1)
913  enddo
914  do k = ks, ke
915  rhou_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
916  enddo
917 #endif
918 
919  !--- momentum y
920 #ifdef _OPENACC
921  flx0 = mflx(k-1) * v(k)
922  if ( k < ke ) then
923  flx1 = mflx(k) * v(k+1)
924  end if
925  rhov_t(k) = - ( flx1 - flx0 ) * rcdz(k)
926 
927  end do
928 #else
929  do k = ks-1, ke-1
930  flx(k) = mflx(k) * v(k+1)
931  enddo
932  do k = ks, ke
933  rhov_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
934  enddo
935 #endif
936 
937  return
938  end subroutine atmos_phy_mp_precipitation_momentum
939 
940 end module scale_atmos_phy_mp_common
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:152
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:174
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:149
scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer
subroutine, public atmos_phy_mp_negative_fixer(KA, KS, KE, IA, IS, IE, JA, JS, JE, QLA, QIA, limit_negative, DENS, TEMP, CVtot, CPtot, QV, QTRC, RHOH, DENS_diff, ENGI_diff)
ATMOS_PHY_MP_negative_fixer negative fixer.
Definition: scale_atmos_phy_mp_common.F90:69
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_semilag
subroutine, public atmos_phy_mp_precipitation_semilag(KA, KS, KE, QHA, QLA, QIA, TEMP, vterm, FZ, FDZ, RCDZ, dt, i, j, DENS, RHOQ, CPtot, CVtot, RHOE, mflx, sflx, esflx)
Definition: scale_atmos_phy_mp_common.F90:641
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_hydrometeor::lhv
real(rp), public lhv
latent heat of vaporization for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:144
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:61
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_momentum
subroutine, public atmos_phy_mp_precipitation_momentum(KA, KS, KE, DENS, MOMZ, U, V, mflx, RCDZ, RFDZ, MOMZ_t, RHOU_t, RHOV_t)
Definition: scale_atmos_phy_mp_common.F90:852
scale_atmos_phy_mp_common
module ATMOSPHERE / Physics Cloud Microphysics - Common
Definition: scale_atmos_phy_mp_common.F90:13
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:60
scale_atmos_phy_mp_common::atmos_phy_mp_saturation_adjustment_3d
subroutine atmos_phy_mp_saturation_adjustment_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, flag_liquid, TEMP, QV, QC, QI, CPtot, CVtot, RHOE_d)
Saturation adjustment.
Definition: scale_atmos_phy_mp_common.F90:283
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_upwind
subroutine, public atmos_phy_mp_precipitation_upwind(KA, KS, KE, QHA, QLA, QIA, TEMP, vterm, FDZ, RCDZ, dt, i, j, DENS, RHOQ, CPtot, CVtot, RHOE, mflx, sflx, esflx)
Definition: scale_atmos_phy_mp_common.F90:451
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_hydrometeor::lhf
real(rp), public lhf
latent heat of fusion for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:146
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:246
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_hydrometeor::cp_vapor
real(rp), public cp_vapor
CP for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:150
scale_atmos_hydrometeor::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:151
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:153
scale_atmos_hydrometeor::cp_ice
real(rp), public cp_ice
CP for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:154
Generated by   doxygen 1.8.17