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  do j = js, je
136  do i = is, ie
137 
138  diffq(:) = 0.0_rp
139 
140  do k = ks, ke
141  eng(k) = cvtot(k,i,j) * temp(k,i,j)
142  end do
143 
144  if ( rhoh_out ) then
145  do k = ks, ke
146  eng0(k) = eng(k)
147  end do
148  end if
149 
150  if ( qla > 0 ) then
151  call atmos_hydrometeor_lhv( ka, ks, ke, temp(:,i,j), lhv(:) )
152 
153  do iq = 1, qla
154  do k = ks, ke
155  work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to liq)
156  if ( work > 0.0_rp ) then
157  eng(k) = eng(k) + work * lhv(k)
158  diffq(k) = diffq(k) - work
159  cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_water - cv_vapor )
160  cptot(k,i,j) = cptot(k,i,j) + work * ( cp_water - cp_vapor )
161  ! remove negative value of hydrometeors (mass)
162  qtrc(k,i,j,iq) = 0.0_rp
163  end if
164  enddo
165  enddo
166  end if
167 
168  if ( qia > 0 ) then
169  call atmos_hydrometeor_lhs( ka, ks, ke, temp(:,i,j), lhs(:) )
170 
171  do iq = qla+1, qla+qia
172  do k = ks, ke
173  work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to ice)
174  if ( work > 0.0_rp ) then
175  eng(k) = eng(k) + work * lhs(k)
176  diffq(k) = diffq(k) - work
177  cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_ice - cv_vapor )
178  cptot(k,i,j) = cptot(k,i,j) + work * ( cp_ice - cp_vapor )
179  ! remove negative value of hydrometeors (mass)
180  qtrc(k,i,j,iq) = 0.0_rp
181  end if
182  enddo
183  enddo
184  end if
185 
186 
187  if ( abs(limit_negative) > 0.0_rp ) then
188  do k = ks, ke
189  if ( diffq(k) < - abs(limit_negative) ) then
190  diffq_min = min( diffq_min, diffq(k) )
191  log_error("ATMOS_PHY_MP_negative_fixer",*) 'large negative is found'
192  log_error_cont(*) 'value = ', diffq(k), ' at (', k, ',', i, ',', j, ')'
193  end if
194  end do
195  end if
196 
197  do k = ks, ke
198  if ( diffq(k) < 0.0_rp ) then
199  ! Compensate for the lack of hydrometeors by the water vapor
200  qv(k,i,j) = qv(k,i,j) + diffq(k)
201  end if
202  end do
203 
204  if ( rhoh_out ) then
205  do k = ks, ke
206  rhoh(k,i,j) = ( eng(k) - eng0(k) ) * dens(k,i,j)
207  end do
208  end if
209 
210 
211  if ( dens_out ) then
212  do k = ks, ke
213  dens0(k) = dens(k,i,j)
214  end do
215  end if
216  if ( engi_out ) then
217  do k = ks, ke
218  eng0(k) = eng(k) * dens(k,i,j)
219  end do
220  end if
221 
222  ! fix negative QV
223  do k = ks, ke
224  work = - min( qv(k,i,j), 0.0_rp )
225  if ( work > 0.0_rp ) then
226  qv(k,i,j) = 0.0_rp
227  dens(k,i,j) = dens(k,i,j) * ( 1.0_rp + work ) ! not to change mass of dry air
228  cvtot(k,i,j) = ( cvtot(k,i,j) + work * cv_vapor ) / ( 1.0_rp + work )
229  cptot(k,i,j) = ( cptot(k,i,j) + work * cp_vapor ) / ( 1.0_rp + work )
230  eng(k) = ( eng(k) + work * cv_vapor * temp(k,i,j) ) / ( 1.0_rp + work )
231  end if
232  enddo
233 
234  do k = ks, ke
235  ! update
236  temp(k,i,j) = eng(k) / cvtot(k,i,j)
237  end do
238 
239  if ( dens_out ) then
240  do k = ks, ke
241  dens_diff(k,i,j) = dens(k,i,j) - dens0(k)
242  end do
243  end if
244  if ( engi_out ) then
245  do k = ks, ke
246  engi_diff(k,i,j) = eng(k) * dens(k,i,j) - eng0(k)
247  end do
248  end if
249 
250  enddo
251  enddo
252 
253 
254  if ( abs(limit_negative) > 0.0_rp &
255  .AND. abs(limit_negative) < abs(diffq_min) ) then
256  log_error_cont(*) 'maximum negative hydrometeor ', diffq_min, ' < ', - abs(limit_negative)
257  call prc_abort
258  endif
259 
260  call prof_rapend('MP_filter', 3)
261 
262  return
263  end subroutine atmos_phy_mp_negative_fixer
264 
265  !-----------------------------------------------------------------------------
267  !-----------------------------------------------------------------------------
268  subroutine atmos_phy_mp_saturation_adjustment_3d( &
269  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
270  DENS, &
271  flag_liquid, &
272  TEMP, &
273  QV, QC, QI, &
274  CPtot, CVtot, &
275  RHOE_d )
276  use scale_prc, only: &
277  prc_abort
278  use scale_atmos_hydrometeor, only: &
279  cp_vapor, &
280  cp_water, &
281  cp_ice, &
282  cv_vapor, &
283  cv_water, &
284  cv_ice, &
285  lhv, &
286  lhf
287  use scale_atmos_saturation, only: &
288  atmos_saturation_moist_conversion_dens_all, &
289  atmos_saturation_moist_conversion_dens_liq
290  implicit none
291 
292  integer, intent(in) :: KA, KS, KE
293  integer, intent(in) :: IA, IS, IE
294  integer, intent(in) :: JA, JS, JE
295 
296  real(RP), intent(in) :: DENS(KA,IA,JA)
297  logical , intent(in) :: flag_liquid
298 
299  real(RP), intent(inout) :: TEMP (KA,IA,JA)
300  real(RP), intent(inout) :: QV (KA,IA,JA)
301  real(RP), intent(inout) :: QC (KA,IA,JA)
302  real(RP), intent(inout) :: QI (KA,IA,JA)
303  real(RP), intent(inout) :: CPtot(KA,IA,JA)
304  real(RP), intent(inout) :: CVtot(KA,IA,JA)
305 
306  real(RP), intent(out) :: RHOE_d(KA,IA,JA)
307 
308  ! working
309  real(RP) :: QV1
310  real(RP) :: QC1
311  real(RP) :: QI1
312 
313  real(RP) :: Emoist ! moist internal energy
314 
315  logical :: converged, error
316 
317  integer :: k, i, j
318  !---------------------------------------------------------------------------
319 
320  call prof_rapstart('MP_Saturation_adjustment', 2)
321 
322  error = .false.
323 
324  if ( flag_liquid ) then ! warm rain
325 
326  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
327  !$omp shared (IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID) &
328  !$omp shared (KS,KE,IS,IE,JS,JE, &
329  !$omp CP_VAPOR,CP_WATER,CV_VAPOR,CV_WATER,LHV,LHF, &
330  !$omp DENS,QV,QC,TEMP,CPtot,CVtot,RHOE_d,error) &
331  !$omp private(i,j,k, &
332  !$omp QV1,QC1,Emoist,converged)
333  do j = js, je
334  do i = is, ie
335  do k = ks, ke
336 
337  qv1 = qv(k,i,j)
338  qc1 = qc(k,i,j)
339 
340  emoist = temp(k,i,j) * cvtot(k,i,j) &
341  + qv1 * lhv
342 
343  call atmos_saturation_moist_conversion_dens_liq( dens(k,i,j), emoist, & ! [IN]
344  temp(k,i,j), qv1, qc1, & ! [INOUT]
345  cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
346  converged ) ! [OUT]
347 
348  if ( .NOT. converged ) then
349  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j
350  error = .true.
351  exit
352  endif
353 
354  rhoe_d(k,i,j) = - lhv * ( qv1 - qv(k,i,j) ) * dens(k,i,j)
355 
356  qv(k,i,j) = qv1
357  qc(k,i,j) = qc1
358 
359  end do
360  end do
361  end do
362 
363  if ( error ) call prc_abort
364 
365  else ! cold rain
366 
367  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
368  !$omp shared (IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID) &
369  !$omp shared (KS,KE,IS,IE,JS,JE, &
370  !$omp CP_VAPOR,CP_WATER,CP_ICE,CV_VAPOR,CV_WATER,CV_ICE,LHV,LHF, &
371  !$omp DENS,QV,QC,QI,TEMP,CPtot,CVtot,RHOE_d,error) &
372  !$omp private(i,j,k, &
373  !$omp QV1,QC1,QI1,Emoist,converged)
374  do j = js, je
375  do i = is, ie
376  do k = ks, ke
377  qv1 = qv(k,i,j)
378  qc1 = qc(k,i,j)
379  qi1 = qi(k,i,j)
380 
381  emoist = temp(k,i,j) * cvtot(k,i,j) &
382  + qv1 * lhv &
383  - qi1 * lhf
384 
385  call atmos_saturation_moist_conversion_dens_all( dens(k,i,j), emoist, & ! [IN]
386  temp(k,i,j), qv1, qc1, qi1, & ! [INOUT]
387  cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
388  converged ) ! [OUT]
389 
390  if ( .NOT. converged ) then
391  log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j
392  error = .true.
393  exit
394  endif
395 
396  rhoe_d(k,i,j) = ( - lhv * ( qv1 - qv(k,i,j) ) &
397  + lhf * ( qi1 - qi(k,i,j) ) ) * dens(k,i,j)
398 
399  qv(k,i,j) = qv1
400  qc(k,i,j) = qc1
401  qi(k,i,j) = qi1
402 
403  end do
404  end do
405  end do
406 
407  if ( error ) call prc_abort
408  endif
409 
410  call prof_rapend ('MP_Saturation_adjustment', 2)
411 
412  return
413  end subroutine atmos_phy_mp_saturation_adjustment_3d
414 
415  !-----------------------------------------------------------------------------
416 !OCL SERIAL
417  subroutine atmos_phy_mp_precipitation_upwind( &
418  KA, KS, KE, QHA, QLA, QIA, &
419  TEMP, vterm, FDZ, RCDZ, dt, &
420  i, j, &
421  DENS, RHOQ, CPtot, CVtot, RHOE, &
422  mflx, sflx, esflx )
423  use scale_const, only: &
424  grav => const_grav
425  use scale_atmos_hydrometeor, only: &
426  cp_water, &
427  cp_ice, &
428  cv_water, &
429  cv_ice
430  implicit none
431  integer, intent(in) :: ka, ks, ke
432  integer, intent(in) :: qha, qla, qia
433 
434  real(rp), intent(in) :: temp (ka)
435  real(rp), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
436  real(rp), intent(in) :: fdz (ka)
437  real(rp), intent(in) :: rcdz (ka)
438  real(dp), intent(in) :: dt
439  integer, intent(in) :: i, j ! for debug
440 
441  real(rp), intent(inout) :: dens (ka)
442  real(rp), intent(inout) :: rhoq (ka,qha)
443  real(rp), intent(inout) :: cptot(ka)
444  real(rp), intent(inout) :: cvtot(ka)
445  real(rp), intent(inout) :: rhoe (ka)
446 
447  real(rp), intent(out) :: mflx (ka)
448  real(rp), intent(out) :: sflx (2)
449  real(rp), intent(out) :: esflx
450 
451  real(rp) :: vtermh(ka)
452  real(rp) :: qflx (ka)
453  real(rp) :: eflx (ka)
454  real(rp) :: rhocp (ka)
455  real(rp) :: rhocv (ka)
456  real(rp) :: ddens
457  real(rp) :: cp, cv
458 
459  integer :: k, iq
460  !---------------------------------------------------------------------------
461 
462  ! tracer/energy transport by falldown
463  ! 1st order upwind, forward euler, velocity is always negative
464 
465  mflx(:) = 0.0_rp
466  sflx(:) = 0.0_rp
467  esflx = 0.0_rp
468  qflx(ke) = 0.0_rp
469  eflx(ke) = 0.0_rp
470 
471  do k = ks, ke
472  rhocp(k) = cptot(k) * dens(k)
473  rhocv(k) = cvtot(k) * dens(k)
474  end do
475 
476  do iq = 1, qha
477  do k = ks, ke-1
478  vtermh(k) = 0.5_rp * ( vterm(k+1,iq) + vterm(k,iq) )
479  enddo
480  vtermh(ks-1) = vterm(ks,iq)
481 
482  !--- mass flux for each tracer, upwind with vel < 0
483  do k = ks-1, ke-1
484  qflx(k) = vtermh(k) * rhoq(k+1,iq)
485  enddo
486 
487  !--- update falling tracer
488  do k = ks, ke
489  rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
490  enddo ! falling (water mass & number) tracer
491 
492  ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
493  if ( iq > qla + qia ) cycle
494 
495  do k = ks-1, ke-1
496  mflx(k) = mflx(k) + qflx(k)
497  end do
498 
499  if ( iq > qla ) then ! ice water
500  cp = cp_ice
501  cv = cv_ice
502  sflx(2) = sflx(2) + qflx(ks-1)
503  else ! liquid water
504  cp = cp_water
505  cv = cv_water
506  sflx(1) = sflx(1) + qflx(ks-1)
507  end if
508 
509  !--- update density
510  do k = ks, ke
511  ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
512  rhocp(k) = rhocp(k) + cp * ddens
513  rhocv(k) = rhocv(k) + cv * ddens
514  dens(k) = dens(k) + ddens
515  end do
516 
517  ! internal energy flux
518  do k = ks-1, ke-1
519  eflx(k) = qflx(k) * temp(k+1) * cv
520  end do
521  esflx = esflx + eflx(ks-1)
522 
523  !--- update internal energy
524  do k = ks, ke
525  rhoe(k) = rhoe(k) - ( ( eflx(k) - eflx(k-1) ) & ! contribution with the transport of internal energy
526  + qflx(k) * fdz(k) * grav & ! contribution with the release of potential energy
527  ) * rcdz(k) * dt
528  end do
529 
530  end do
531 
532  do k = ks, ke
533  cptot(k) = rhocp(k) / dens(k)
534  cvtot(k) = rhocv(k) / dens(k)
535  end do
536 
537  return
538  end subroutine atmos_phy_mp_precipitation_upwind
539 
540  !-----------------------------------------------------------------------------
541 !OCL SERIAL
542  subroutine atmos_phy_mp_precipitation_semilag( &
543  KA, KS, KE, QHA, QLA, QIA, &
544  TEMP, vterm, FZ, FDZ, RCDZ, dt, &
545  i, j, &
546  DENS, RHOQ, CPtot, CVtot, RHOE, &
547  mflx, sflx, esflx )
548  use scale_const, only: &
549  grav => const_grav
550  use scale_atmos_hydrometeor, only: &
551  cp_water, &
552  cp_ice, &
553  cv_water, &
554  cv_ice
555  implicit none
556  integer, intent(in) :: ka, ks, ke
557  integer, intent(in) :: qha, qla, qia
558 
559  real(rp), intent(in) :: temp (ka)
560  real(rp), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
561  real(rp), intent(in) :: fz (ka)
562  real(rp), intent(in) :: fdz (ka)
563  real(rp), intent(in) :: rcdz (ka)
564  real(dp), intent(in) :: dt
565  integer, intent(in) :: i, j ! for debug
566 
567  real(rp), intent(inout) :: dens (ka)
568  real(rp), intent(inout) :: rhoq (ka,qha)
569  real(rp), intent(inout) :: cptot(ka)
570  real(rp), intent(inout) :: cvtot(ka)
571  real(rp), intent(inout) :: rhoe (ka)
572 
573  real(rp), intent(out) :: mflx (ka)
574  real(rp), intent(out) :: sflx (2)
575  real(rp), intent(out) :: esflx
576 
577  real(rp) :: qflx(ka)
578  real(rp) :: eflx(ka)
579  real(rp) :: rhocp(ka)
580  real(rp) :: rhocv(ka)
581  real(rp) :: ddens
582  real(rp) :: cp, cv
583 
584  real(rp) :: vtermh(ka)
585  real(rp) :: dvterm(ka)
586  real(rp) :: cdz (ka)
587  real(rp) :: rfdz2 (ka)
588  real(rp) :: dist (ka)
589  real(rp) :: z_src
590  real(rp) :: flx
591  integer :: k_src (ka)
592  integer :: k_dst
593 
594  integer :: k, iq
595  !---------------------------------------------------------------------------
596 
597  ! tracer/energy transport by falldown
598  ! velocity is always negative
599 
600  mflx(:) = 0.0_rp
601  sflx(:) = 0.0_rp
602  eflx(:) = 0.0_rp
603 
604  do k = ks, ke
605  rhocp(k) = cptot(k) * dens(k)
606  rhocv(k) = cvtot(k) * dens(k)
607  end do
608 
609  do k = ks, ke
610  cdz(k) = 1.0_rp / rcdz(k)
611  end do
612  do k = ks, ke-1
613  rfdz2(k) = 1.0_rp / ( cdz(k) + cdz(k+1) )
614  end do
615 
616  do iq = 1, qha
617 
618  qflx(:) = 0.0_rp
619 
620  do k = ks, ke-1
621  vtermh(k) = ( cdz(k) * vterm(k+1,iq) + cdz(k+1) * vterm(k,iq) ) * rfdz2(k)
622  enddo
623  vtermh(ks-1) = vterm(ks,iq)
624 
625  do k = ks, ke
626  dvterm(k) = vtermh(k) - vtermh(k-1)
627  enddo
628 
629  ! Movement distance of the cell wall by the fall
630  ! the midpoint method (second-order Runge-Kutta)
631  ! dz/dt = v(z + v dt/2) ~ v(z) + v dt/2 dv/dz + 1/2 (v dt/2)^2 d^2v/dz^2
632  do k = ks, ke-1
633  dist(k) = - vtermh(k) * dt &
634  + vtermh(k) * dt**2 / 2.0_rp * ( dvterm(k+1)+dvterm(k) ) * rfdz2(k) &
635  - vtermh(k)**2 * dt**3 / 4.0_rp * ( dvterm(k+1)*rcdz(k+1) - dvterm(k)*rcdz(k) ) * rfdz2(k)
636  dist(k) = max( dist(k), 0.0_rp )
637  enddo
638  dist(ks-1) = - vtermh(ks-1) * dt &
639  + vtermh(ks-1) * dt**2 / 2.0_rp * dvterm(ks)*rcdz(ks)
640  dist(ks-1) = max( dist(ks-1), 0.0_rp )
641 
642  ! wall cannot overtake
643  do k = ke-2, ks-1, -1
644  dist(k) = min( dist(k), dist(k+1) + cdz(k+1) )
645  end do
646 
647 ! LOG_INFO_CONT(*) "distance", iq
648 ! do k = KA, 1, -1
649 ! LOG_INFO_CONT('(1x,I5,3F9.3,ES15.5)') k, dist(k), vtermh(k), vterm(k,iq), RHOQ(k,iq)
650 ! enddo
651 
652  ! search number of source cell
653  do k_dst = ks-1, ke-1
654  z_src = fz(k_dst) + dist(k_dst)
655 
656  k_src(k_dst) = k_dst
657  do k = k_dst, ke-1
658  if ( z_src > fz(k ) &
659  .AND. z_src <= fz(k+1) ) then
660  k_src(k_dst) = k
661  exit
662  endif
663  enddo
664  if ( z_src > fz(ke) ) k_src(k_dst) = ke
665  enddo
666 
667 ! LOG_INFO_CONT(*) "seek", iq
668 ! do k = KA, 1, -1
669 ! LOG_INFO_CONT('(1x,2I5,2F9.3)') k, k_src(k), FZ(k), FZ(k)+dist(k)
670 ! enddo
671 
672  if ( iq > qla ) then ! ice water
673  cp = cp_ice
674  cv = cv_ice
675  else ! liquid water
676  cp = cp_water
677  cv = cv_water
678  end if
679 
680  do k_dst = ks-1, ke-1
681  do k = k_dst, k_src(k_dst)-1
682  flx = rhoq(k+1,iq) * cdz(k+1) / dt ! sum column mass rhoq*dz
683  qflx(k_dst) = qflx(k_dst) - flx
684  eflx(k_dst) = eflx(k_dst) - flx * temp(k+1) * cv ! internal energy flux
685  dist(k_dst) = dist(k_dst) - cdz(k+1) ! residual
686  enddo
687  if ( k_src(k_dst) < ke ) then
688  ! residual (simple upwind)
689  flx = rhoq(k_src(k_dst)+1,iq) * dist(k_dst) / dt
690  qflx(k_dst) = qflx(k_dst) - flx ! sum column mass rhoq*dz
691  eflx(k_dst) = eflx(k_dst) -flx * temp(k_src(k_dst)+1) * cv ! internal energy flux
692  end if
693  enddo
694 
695 ! LOG_INFO_CONT(*) "flux", iq
696 ! do k = KA, 1, -1
697 ! LOG_INFO_CONT('(1x,2I5,F9.3,2ES15.5)') k, k_src(k), dist(k), qflx(k), vtermh(k)*RHOQ(k+1,iq)
698 ! enddo
699 
700  !--- update falling tracer
701  do k = ks, ke
702  rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
703  enddo ! falling (water mass & number) tracer
704 
705 ! LOG_INFO_CONT(*) "tendency", iq
706 ! do k = KA, 1, -1
707 ! LOG_INFO_CONT('(1x,I5,ES15.5)') k, - dt * ( qflx(k) - qflx(k-1) ) * RCDZ(k)
708 ! enddo
709 
710  ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
711  if ( iq > qla + qia ) cycle
712 
713  do k = ks-1, ke-1
714  mflx(k) = mflx(k) + qflx(k)
715  end do
716 
717  if ( iq > qla ) then ! ice water
718  sflx(2) = sflx(2) + qflx(ks-1)
719  else ! liquid water
720  sflx(1) = sflx(1) + qflx(ks-1)
721  end if
722 
723  !--- update density
724  do k = ks, ke
725  ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
726  rhocp(k) = rhocp(k) + cp * ddens
727  rhocv(k) = rhocv(k) + cv * ddens
728  dens(k) = dens(k) + ddens
729  end do
730 
731  !--- update internal energy
732  do k = ks, ke
733  rhoe(k) = rhoe(k) - ( ( eflx(k) - eflx(k-1) ) & ! contribution with the transport of internal energy
734  + qflx(k) * fdz(k) * grav & ! contribution with the release of potential energy
735  ) * rcdz(k) * dt
736  end do
737  esflx = eflx(ks-1)
738 
739  end do
740 
741  do k = ks, ke
742  cptot(k) = rhocp(k) / dens(k)
743  cvtot(k) = rhocv(k) / dens(k)
744  end do
745 
746  return
747  end subroutine atmos_phy_mp_precipitation_semilag
748 
749  !-----------------------------------------------------------------------------
750 !OCL SERIAL
751  subroutine atmos_phy_mp_precipitation_momentum( &
752  KA, KS, KE, &
753  DENS, MOMZ, U, V, mflx, &
754  RCDZ, RFDZ, &
755  MOMZ_t, RHOU_t, RHOV_t )
756  implicit none
757 
758  integer, intent(in) :: ka, ks, ke
759  real(rp), intent(in) :: dens (ka)
760  real(rp), intent(in) :: momz (ka)
761  real(rp), intent(in) :: u (ka)
762  real(rp), intent(in) :: v (ka)
763  real(rp), intent(in) :: mflx (ka)
764  real(rp), intent(in) :: rcdz (ka)
765  real(rp), intent(in) :: rfdz (ka)
766  real(rp), intent(out) :: momz_t(ka)
767  real(rp), intent(out) :: rhou_t(ka)
768  real(rp), intent(out) :: rhov_t(ka)
769 
770  real(rp) :: flx(ka)
771 
772  integer :: k
773  !---------------------------------------------------------------------------
774 
775  flx(ke) = 0.0_rp
776 
777  !--- momentum z (half level)
778  do k = ks, ke-1
779  flx(k) = ( mflx(k) + mflx(k-1) ) * momz(k) / ( dens(k+1) + dens(k) )
780  enddo
781  do k = ks, ke-1
782  momz_t(k) = - ( flx(k+1) - flx(k) ) * rfdz(k)
783  enddo
784  momz_t(ke) = 0.0_rp
785 
786  !--- momentum x
787  do k = ks-1, ke-1
788  flx(k) = mflx(k) * u(k+1)
789  enddo
790  do k = ks, ke
791  rhou_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
792  enddo
793 
794  !--- momentum y
795  do k = ks-1, ke-1
796  flx(k) = mflx(k) * v(k+1)
797  enddo
798  do k = ks, ke
799  rhov_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
800  enddo
801 
802  return
803  end subroutine atmos_phy_mp_precipitation_momentum
804 
805 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:46
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:342
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:133
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:159
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:130
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:548
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:44
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:126
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:57
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:756
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:56
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:276
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:423
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_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
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:128
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
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:131
scale_atmos_hydrometeor::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:132
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:134
scale_atmos_hydrometeor::cp_ice
real(rp), public cp_ice
CP for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:135
Generated by   doxygen 1.8.17