SCALE-RM
scale_atmos_phy_mp_common.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
14 !-------------------------------------------------------------------------------
15 #include "inc_openmp.h"
16 module scale_atmos_phy_mp_common
17  !-----------------------------------------------------------------------------
18  !
19  !++ used modules
20  !
21  use scale_precision
22  use scale_stdio
23  use scale_prof
24  use scale_grid_index
25  use scale_tracer
26  use scale_const, only: &
27  undef => const_undef
28  !-----------------------------------------------------------------------------
29  implicit none
30  private
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public procedure
34  !
35  public :: atmos_phy_mp_negative_fixer
36  public :: atmos_phy_mp_saturation_adjustment
37  public :: atmos_phy_mp_precipitation
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Public parameters & variables
42  !
43  !-----------------------------------------------------------------------------
44  !
45  !++ Private procedure
46  !
47  private :: moist_conversion_liq
48  private :: moist_conversion_all
49 
50  !-----------------------------------------------------------------------------
51  !
52  !++ Private parameters & variables
53  !
54  !-----------------------------------------------------------------------------
55 contains
56 
57  !-----------------------------------------------------------------------------
59  !-----------------------------------------------------------------------------
60  subroutine atmos_phy_mp_negative_fixer( &
61  DENS, &
62  RHOT, &
63  QTRC )
64  implicit none
65 
66  real(RP), intent(inout) :: DENS(ka,ia,ja)
67  real(RP), intent(inout) :: RHOT(ka,ia,ja)
68  real(RP), intent(inout) :: QTRC(ka,ia,ja,qa)
69 
70  real(RP) :: diffq
71 
72  integer :: k, i, j, iq
73  !---------------------------------------------------------------------------
74 
75  call prof_rapstart('MP_filter', 3)
76 
77  !$omp parallel do private(i,j,diffq) OMP_SCHEDULE_ collapse(2)
78  do j = 1, ja
79  do i = 1, ia
80  do k = ks, ke
81 
82  diffq = 0.0_rp
83  do iq = qqs+1, qqe
84  ! total hydrometeor (before correction)
85  diffq = diffq + qtrc(k,i,j,iq)
86  ! remove negative value of hydrometeors (mass)
87  qtrc(k,i,j,iq) = max( qtrc(k,i,j,iq), 0.0_rp )
88  enddo
89 
90  do iq = qqs+1, qqe
91  ! difference between before and after correction
92  diffq = diffq - qtrc(k,i,j,iq)
93  enddo
94 
95  ! Compensate for the lack of hydrometeors by the water vapor
96  qtrc(k,i,j,i_qv) = qtrc(k,i,j,i_qv) + diffq
97 
98  ! TODO: We have to consider energy conservation (but very small)
99 
100  ! remove negative value of water vapor (mass)
101  diffq = qtrc(k,i,j,i_qv)
102  qtrc(k,i,j,i_qv) = max( qtrc(k,i,j,i_qv), 0.0_rp )
103  diffq = diffq - qtrc(k,i,j,i_qv)
104 
105  ! Apply correction to total density
106  ! TODO: We have to consider energy conservation (but very small)
107  dens(k,i,j) = dens(k,i,j) * ( 1.0_rp - diffq ) ! diffq is negative
108  rhot(k,i,j) = rhot(k,i,j) * ( 1.0_rp - diffq )
109 
110  enddo
111  enddo
112  enddo
113 
114  call prof_rapend('MP_filter', 3)
115 
116  return
117  end subroutine atmos_phy_mp_negative_fixer
118 
119  !-----------------------------------------------------------------------------
121  !-----------------------------------------------------------------------------
122  subroutine atmos_phy_mp_saturation_adjustment( &
123  RHOE_t, &
124  QTRC_t, &
125  RHOE0, &
126  QTRC0, &
127  DENS0, &
128  flag_liquid )
129  use scale_const, only: &
130  lhv => const_lhv, &
131  lhf => const_lhf
132  use scale_time, only: &
133  dt => time_dtsec_atmos_phy_mp
134  use scale_atmos_thermodyn, only: &
135  thermodyn_qd => atmos_thermodyn_qd, &
136  thermodyn_cv => atmos_thermodyn_cv, &
137  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
138  use scale_atmos_saturation, only: &
139  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
140  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
141  implicit none
142 
143  real(RP), intent(inout) :: RHOE_t(ka,ia,ja) ! tendency rhoe [J/m3/s]
144  real(RP), intent(inout) :: QTRC_t(ka,ia,ja,qa) ! tendency tracer [kg/kg/s]
145  real(RP), intent(inout) :: RHOE0 (ka,ia,ja) ! density * internal energy [J/m3]
146  real(RP), intent(inout) :: QTRC0 (ka,ia,ja,qa) ! mass concentration [kg/kg]
147  real(RP), intent(in) :: DENS0 (ka,ia,ja) ! density [kg/m3]
148  logical, intent(in) :: flag_liquid ! use scheme only for the liquid water?
149 
150  ! working
151  real(RP) :: TEMP0 (ka,ia,ja)
152  real(RP) :: PRES0 (ka,ia,ja)
153  real(RP) :: QDRY0 (ka,ia,ja)
154  real(RP) :: CVtot (ka,ia,ja)
155 
156  real(RP) :: Emoist(ka,ia,ja) ! moist internal energy
157  real(RP) :: QSUM1 (ka,ia,ja) ! QV+QC+QI
158  real(RP) :: TEMP1 (ka,ia,ja)
159 
160  real(RP) :: RHOE1 (ka,ia,ja)
161  real(RP) :: QTRC1 (ka,ia,ja,qa)
162  real(RP) :: rdt
163 
164  integer :: k, i, j, iq
165  !---------------------------------------------------------------------------
166 
167 #ifndef DRY
168 
169  call prof_rapstart('MP_Saturation_adjustment', 2)
170 
171  rdt = 1.0_rp / dt
172 
173  !$omp parallel do private(i,j,k,iq) OMP_SCHEDULE_ collapse(4)
174  do iq = qqs, qqe
175  do j = jsb, jeb
176  do i = isb, ieb
177  do k = ks, ke
178  qtrc1(k,i,j,iq) = qtrc0(k,i,j,iq)
179  enddo
180  enddo
181  enddo
182  enddo
183 
184  call thermodyn_temp_pres_e( temp0(:,:,:), & ! [OUT]
185  pres0(:,:,:), & ! [OUT]
186  dens0(:,:,:), & ! [IN]
187  rhoe0(:,:,:), & ! [IN]
188  qtrc0(:,:,:,:) ) ! [IN]
189 
190  ! qdry dont change through the process
191  call thermodyn_qd( qdry0(:,:,:), & ! [OUT]
192  qtrc0(:,:,:,:) ) ! [IN]
193 
194  call thermodyn_cv( cvtot(:,:,:), & ! [OUT]
195  qtrc0(:,:,:,:), & ! [IN]
196  qdry0(:,:,:) ) ! [IN]
197 
198  if ( i_qi <= 0 .OR. flag_liquid ) then ! warm rain
199 
200  ! Turn QC into QV with consistency of moist internal energy
201  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
202  do j = jsb, jeb
203  do i = isb, ieb
204  do k = ks, ke
205  emoist(k,i,j) = temp0(k,i,j) * cvtot(k,i,j) &
206  + qtrc1(k,i,j,i_qv) * lhv
207 
208  qsum1(k,i,j) = qtrc1(k,i,j,i_qv) &
209  + qtrc1(k,i,j,i_qc)
210 
211  qtrc1(k,i,j,i_qv) = qsum1(k,i,j)
212  qtrc1(k,i,j,i_qc) = 0.0_rp
213  enddo
214  enddo
215  enddo
216 
217  call thermodyn_cv( cvtot(:,:,:), & ! [OUT]
218  qtrc1(:,:,:,:), & ! [IN]
219  qdry0(:,:,:) ) ! [IN]
220 
221  ! new temperature (after QC evaporation)
222  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
223  do j = jsb, jeb
224  do i = isb, ieb
225  do k = ks, ke
226  temp1(k,i,j) = ( emoist(k,i,j) - qtrc1(k,i,j,i_qv) * lhv ) / cvtot(k,i,j)
227  enddo
228  enddo
229  enddo
230 
231  call moist_conversion_liq( temp1(:,:,:), & ! [INOUT]
232  qtrc1(:,:,:,:), & ! [INOUT]
233  dens0(:,:,:), & ! [IN]
234  qsum1(:,:,:), & ! [IN]
235  qdry0(:,:,:), & ! [IN]
236  emoist(:,:,:) ) ! [IN]
237 
238  else ! cold rain
239 
240  ! Turn QC & QI into QV with consistency of moist internal energy
241  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
242  do j = jsb, jeb
243  do i = isb, ieb
244  do k = ks, ke
245  emoist(k,i,j) = temp0(k,i,j) * cvtot(k,i,j) &
246  + qtrc1(k,i,j,i_qv) * lhv &
247  - qtrc1(k,i,j,i_qi) * lhf
248 
249  qsum1(k,i,j) = qtrc1(k,i,j,i_qv) &
250  + qtrc1(k,i,j,i_qc) &
251  + qtrc1(k,i,j,i_qi)
252 
253  qtrc1(k,i,j,i_qv) = qsum1(k,i,j)
254  qtrc1(k,i,j,i_qc) = 0.0_rp
255  qtrc1(k,i,j,i_qi) = 0.0_rp
256  enddo
257  enddo
258  enddo
259 
260  call thermodyn_cv( cvtot(:,:,:), & ! [OUT]
261  qtrc1(:,:,:,:), & ! [IN]
262  qdry0(:,:,:) ) ! [IN]
263 
264  ! new temperature (after QC & QI evaporation)
265  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
266  do j = jsb, jeb
267  do i = isb, ieb
268  do k = ks, ke
269  temp1(k,i,j) = ( emoist(k,i,j) - qtrc1(k,i,j,i_qv) * lhv ) / cvtot(k,i,j)
270  enddo
271  enddo
272  enddo
273 
274  call moist_conversion_all( temp1(:,:,:), & ! [INOUT]
275  qtrc1(:,:,:,:), & ! [INOUT]
276  dens0(:,:,:), & ! [IN]
277  qsum1(:,:,:), & ! [IN]
278  qdry0(:,:,:), & ! [IN]
279  emoist(:,:,:) ) ! [IN]
280 
281  endif
282 
283  call thermodyn_cv( cvtot(:,:,:), & ! [OUT]
284  qtrc1(:,:,:,:), & ! [IN]
285  qdry0(:,:,:) ) ! [IN]
286 
287  ! mass & energy update
288  !$omp parallel do private(i,j,k,iq) OMP_SCHEDULE_ collapse(4)
289  do iq = qqs, qqe
290  do j = js, je
291  do i = is, ie
292  do k = ks, ke
293  qtrc_t(k,i,j,iq) = qtrc_t(k,i,j,iq) + ( qtrc1(k,i,j,iq) - qtrc0(k,i,j,iq) ) * rdt
294 
295  qtrc0(k,i,j,iq) = qtrc1(k,i,j,iq)
296  enddo
297  enddo
298  enddo
299  enddo
300 
301  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
302  do j = js, je
303  do i = is, ie
304  do k = ks, ke
305  rhoe1(k,i,j) = dens0(k,i,j) * temp1(k,i,j) * cvtot(k,i,j)
306 
307  rhoe_t(k,i,j) = rhoe_t(k,i,j) + ( rhoe1(k,i,j) - rhoe0(k,i,j) ) * rdt
308 
309  rhoe0(k,i,j) = rhoe1(k,i,j)
310  enddo
311  enddo
312  enddo
313 
314  call prof_rapend ('MP_Saturation_adjustment', 2)
315 
316 #else
317  rhoe_t = undef
318  qtrc_t = undef
319  rhoe0 = undef
320  qtrc0 = undef
321 #endif
322  return
323  end subroutine atmos_phy_mp_saturation_adjustment
324 
325  !-----------------------------------------------------------------------------
327  !-----------------------------------------------------------------------------
328  subroutine moist_conversion_liq( &
329  TEMP1, &
330  QTRC1, &
331  DENS0, &
332  QSUM1, &
333  QDRY0, &
334  Emoist )
335  use scale_const, only: &
336  lhv => const_lhv
337  use scale_process, only: &
338  prc_mpistop
339  use scale_atmos_thermodyn, only: &
340  thermodyn_cv => atmos_thermodyn_cv, &
341  cvw => aq_cv
342  use scale_atmos_saturation, only: &
343  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
344  cvovr_liq, &
345  lovr_liq
346  implicit none
347 
348  real(RP), intent(inout) :: TEMP1 (ka,ia,ja)
349  real(RP), intent(inout) :: QTRC1 (ka,ia,ja,qa)
350  real(RP), intent(in) :: DENS0 (ka,ia,ja)
351  real(RP), intent(in) :: QSUM1 (ka,ia,ja)
352  real(RP), intent(in) :: QDRY0 (ka,ia,ja)
353  real(RP), intent(in) :: Emoist(ka,ia,ja)
354 
355  real(RP) :: QSAT(ka,ia,ja) ! saturated water vapor
356 
357  ! working
358  real(RP) :: temp
359  real(RP) :: q(qa)
360  real(RP) :: CVtot
361  real(RP) :: qsatl_new
362  real(RP) :: Emoist_new ! moist internal energy
363 
364  ! d(X)/dT
365  real(RP) :: dqsatl_dT
366  real(RP) :: dqc_dT
367  real(RP) :: dCVtot_dT
368  real(RP) :: dEmoist_dT
369  real(RP) :: dtemp
370 
371  integer :: ijk_sat
372  integer :: index_sat(ka*ia*ja,3) ! list vector
373 
374  integer, parameter :: itelim = 100
375  real(RP) :: dtemp_criteria
376  logical :: converged
377  integer :: k, i, j, ijk, iq, ite
378  !---------------------------------------------------------------------------
379 
380 #ifndef DRY
381 
382  dtemp_criteria = 0.1_rp**(2+rp/2)
383 
384  call saturation_dens2qsat_liq( qsat(:,:,:), & ! [OUT]
385  temp1(:,:,:), & ! [IN]
386  dens0(:,:,:) ) ! [IN]
387 
388  ijk_sat = 0
389  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
390  do j = js, je
391  do i = is, ie
392  do k = ks, ke
393  if ( qsum1(k,i,j) > qsat(k,i,j) ) then
394  !$omp critical
395  ijk_sat = ijk_sat + 1
396  index_sat(ijk_sat,1) = k
397  index_sat(ijk_sat,2) = i
398  index_sat(ijk_sat,3) = j
399  !$omp end critical
400  endif
401  enddo
402  enddo
403  enddo
404 
405  do ijk = 1, ijk_sat
406  k = index_sat(ijk,1)
407  i = index_sat(ijk,2)
408  j = index_sat(ijk,3)
409 
410  ! store to work
411  temp = temp1(k,i,j)
412  do iq = qqs, qqe
413  q(iq) = qtrc1(k,i,j,iq)
414  enddo
415 
416  converged = .false.
417  do ite = 1, itelim
418 
419  call saturation_dens2qsat_liq( qsatl_new, & ! [OUT]
420  temp, & ! [IN]
421  dens0(k,i,j) ) ! [IN]
422 
423  ! Separation
424  q(i_qv) = qsatl_new
425  q(i_qc) = qsum1(k,i,j) - qsatl_new
426 
427  call thermodyn_cv( cvtot, & ! [OUT]
428  q(:), & ! [IN]
429  qdry0(k,i,j) ) ! [IN]
430 
431  emoist_new = temp * cvtot + qsatl_new * lhv
432 
433  ! dX/dT
434  dqsatl_dt = ( lovr_liq / ( temp*temp ) + cvovr_liq / temp ) * qsatl_new
435 
436  dqc_dt = - dqsatl_dt
437 
438  dcvtot_dt = dqsatl_dt * cvw(i_qv) &
439  + dqc_dt * cvw(i_qc)
440 
441  demoist_dt = qsatl_new * dcvtot_dt + cvtot + dqsatl_dt * lhv
442 
443  dtemp = ( emoist_new - emoist(k,i,j) ) / demoist_dt
444  temp = temp - dtemp
445 
446  if ( abs(dtemp) < dtemp_criteria ) then
447  converged = .true.
448  exit
449  endif
450 
451  if( temp*0.0_rp /= 0.0_rp) exit
452  enddo
453 
454  if ( .NOT. converged ) then
455  write(*,*) 'xxx [moist_conversion] not converged! dtemp=', dtemp,k,i,j,ite
456  call prc_mpistop
457  endif
458 
459  temp1(k,i,j) = temp
460  do iq = qqs, qqe
461  qtrc1(k,i,j,iq) = q(iq)
462  enddo
463  enddo
464 
465 #else
466  temp1 = undef
467  qtrc1 = undef
468 #endif
469 
470  return
471  end subroutine moist_conversion_liq
472 
473  !-----------------------------------------------------------------------------
475  !-----------------------------------------------------------------------------
476  subroutine moist_conversion_all( &
477  TEMP1, &
478  QTRC1, &
479  DENS0, &
480  QSUM1, &
481  QDRY0, &
482  Emoist )
483  use scale_const, only: &
484  lhv => const_lhv, &
485  lhf => const_lhf
486  use scale_process, only: &
487  prc_mpistop
488  use scale_atmos_thermodyn, only: &
489  thermodyn_cv => atmos_thermodyn_cv, &
490  cvw => aq_cv
491  use scale_atmos_saturation, only: &
492  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all, &
493  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
494  saturation_dens2qsat_ice => atmos_saturation_dens2qsat_ice, &
495  saturation_alpha => atmos_saturation_alpha, &
496  saturation_dalphadt => atmos_saturation_dalphadt, &
497  cvovr_liq, &
498  cvovr_ice, &
499  lovr_liq, &
500  lovr_ice
501  implicit none
502 
503  real(RP), intent(inout) :: TEMP1 (ka,ia,ja)
504  real(RP), intent(inout) :: QTRC1 (ka,ia,ja,qa)
505  real(RP), intent(in) :: DENS0 (ka,ia,ja)
506  real(RP), intent(in) :: QSUM1 (ka,ia,ja)
507  real(RP), intent(in) :: QDRY0 (ka,ia,ja)
508  real(RP), intent(in) :: Emoist(ka,ia,ja)
509 
510  real(RP) :: QSAT(ka,ia,ja) ! saturated water vapor
511 
512  ! working
513  real(RP) :: temp
514  real(RP) :: q(qa)
515  real(RP) :: CVtot
516  real(RP) :: alpha
517  real(RP) :: qsat_new, qsatl_new, qsati_new
518  real(RP) :: Emoist_new ! moist internal energy
519 
520  ! d(X)/dT
521  real(RP) :: dalpha_dT
522  real(RP) :: dqsat_dT, dqsatl_dT, dqsati_dT
523  real(RP) :: dqc_dT, dqi_dT
524  real(RP) :: dCVtot_dT
525  real(RP) :: dEmoist_dT
526  real(RP) :: dtemp
527 
528  integer :: ijk_sat
529  integer :: index_sat(ka*ia*ja,3) ! list vector
530 
531  integer, parameter :: itelim = 100
532  real(RP) :: dtemp_criteria
533  logical :: converged
534  integer :: k, i, j, ijk, iq, ite
535  !---------------------------------------------------------------------------
536 
537 #ifndef DRY
538  dtemp_criteria = 0.1_rp**(2+rp/2)
539 
540  call saturation_dens2qsat_all( qsat(:,:,:), & ! [OUT]
541  temp1(:,:,:), & ! [IN]
542  dens0(:,:,:) ) ! [IN]
543 
544  ijk_sat = 0
545  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
546  do j = js, je
547  do i = is, ie
548  do k = ks, ke
549  if ( qsum1(k,i,j) > qsat(k,i,j) ) then
550  !$omp critical
551  ijk_sat = ijk_sat + 1
552  index_sat(ijk_sat,1) = k
553  index_sat(ijk_sat,2) = i
554  index_sat(ijk_sat,3) = j
555  !$omp end critical
556  endif
557  enddo
558  enddo
559  enddo
560 
561  do ijk = 1, ijk_sat
562  k = index_sat(ijk,1)
563  i = index_sat(ijk,2)
564  j = index_sat(ijk,3)
565 
566  ! store to work
567  temp = temp1(k,i,j)
568  do iq = qqs, qqe
569  q(iq) = qtrc1(k,i,j,iq)
570  enddo
571 
572  converged = .false.
573  do ite = 1, itelim
574 
575  ! liquid/ice separation factor
576  call saturation_alpha( alpha, temp )
577  ! Saturation
578  call saturation_dens2qsat_all( qsat_new, temp, dens0(k,i,j) )
579  call saturation_dens2qsat_liq( qsatl_new, temp, dens0(k,i,j) )
580  call saturation_dens2qsat_ice( qsati_new, temp, dens0(k,i,j) )
581 
582  ! Separation
583  q(i_qv) = qsat_new
584  q(i_qc) = ( qsum1(k,i,j)-qsat_new ) * ( alpha )
585  q(i_qi) = ( qsum1(k,i,j)-qsat_new ) * ( 1.0_rp-alpha )
586 
587  call thermodyn_cv( cvtot, & ! [OUT]
588  q(:), & ! [IN]
589  qdry0(k,i,j) ) ! [IN]
590 
591  emoist_new = temp * cvtot + qsat_new * lhv - q(i_qi) * lhf
592 
593  ! dX/dT
594  call saturation_dalphadt( dalpha_dt, temp )
595 
596  dqsatl_dt = ( lovr_liq / ( temp*temp ) + cvovr_liq / temp ) * qsatl_new
597  dqsati_dt = ( lovr_ice / ( temp*temp ) + cvovr_ice / temp ) * qsati_new
598 
599  dqsat_dt = qsatl_new * dalpha_dt + dqsatl_dt * ( alpha ) &
600  - qsati_new * dalpha_dt + dqsati_dt * ( 1.0_rp-alpha )
601 
602  dqc_dt = ( qsum1(k,i,j)-qsat_new ) * dalpha_dt - dqsat_dt * ( alpha )
603  dqi_dt = -( qsum1(k,i,j)-qsat_new ) * dalpha_dt - dqsat_dt * ( 1.0_rp-alpha )
604 
605  dcvtot_dt = dqsat_dt * cvw(i_qv) &
606  + dqc_dt * cvw(i_qc) &
607  + dqi_dt * cvw(i_qi)
608 
609  demoist_dt = temp * dcvtot_dt + cvtot + dqsat_dt * lhv - dqi_dt * lhf
610 
611  dtemp = ( emoist_new - emoist(k,i,j) ) / demoist_dt
612  temp = temp - dtemp
613 
614  if ( abs(dtemp) < dtemp_criteria ) then
615  converged = .true.
616  exit
617  endif
618 
619  if( temp*0.0_rp /= 0.0_rp) exit
620  enddo
621 
622  if ( .NOT. converged ) then
623  write(*,*) temp1(k,i,j), dens0(k,i,j),q,qtrc1(k,i,j,qqs:qqe)
624  write(*,*) 'xxx [moist_conversion] not converged! dtemp=', dtemp, k,i,j,ite
625  call prc_mpistop
626  endif
627 
628  temp1(k,i,j) = temp
629  do iq = qqs, qqe
630  qtrc1(k,i,j,iq) = q(iq)
631  enddo
632  enddo
633 
634 #else
635  temp1 = undef
636  qtrc1 = undef
637 #endif
638 
639  return
640  end subroutine moist_conversion_all
641 
642  !-----------------------------------------------------------------------------
644  !-----------------------------------------------------------------------------
645  subroutine atmos_phy_mp_precipitation( &
646  flux_rain, &
647  flux_snow, &
648  DENS, &
649  MOMZ, &
650  MOMX, &
651  MOMY, &
652  RHOE, &
653  QTRC, &
654  vterm, &
655  temp, &
656  dt )
657  use scale_const, only: &
658  grav => const_grav
659  use scale_grid_real, only: &
660  real_cz, &
661  real_fz
662  use scale_gridtrans, only: &
663  j33g => gtrans_j33g
664  use scale_comm, only: &
665  comm_vars8, &
666  comm_wait
667  use scale_atmos_thermodyn, only: &
668  cvw => aq_cv
669  implicit none
670 
671  real(RP), intent(out) :: flux_rain(ka,ia,ja)
672  real(RP), intent(out) :: flux_snow(ka,ia,ja)
673  real(RP), intent(inout) :: DENS (ka,ia,ja)
674  real(RP), intent(inout) :: MOMZ (ka,ia,ja)
675  real(RP), intent(inout) :: MOMX (ka,ia,ja)
676  real(RP), intent(inout) :: MOMY (ka,ia,ja)
677  real(RP), intent(inout) :: RHOE (ka,ia,ja)
678  real(RP), intent(inout) :: QTRC (ka,ia,ja,qa)
679  real(RP), intent(inout) :: vterm (ka,ia,ja,qa) ! terminal velocity of cloud mass
680  real(RP), intent(in) :: temp (ka,ia,ja)
681  real(DP), intent(in) :: dt
682 
683  real(RP) :: rhoq(ka,ia,ja,qa) ! rho * q before precipitation
684  real(RP) :: qflx(ka,ia,ja,qa)
685  real(RP) :: eflx(ka,ia,ja)
686 
687  real(RP) :: rcdz (ka,ia,ja)
688  real(RP) :: rcdz_u(ka,ia,ja)
689  real(RP) :: rcdz_v(ka,ia,ja)
690  real(RP) :: rfdz (ka,ia,ja)
691 
692  integer :: k, i, j, iq
693  !---------------------------------------------------------------------------
694 
695  call prof_rapstart('MP_Precipitation', 2)
696 
697  if ( tracer_type /= 'SUZUKI10' ) then
698  do iq = i_qc, qqe
699  call comm_vars8( vterm(:,:,:,iq), iq )
700  enddo
701  endif
702  do iq = i_qc, qqe
703  call comm_vars8( qtrc(:,:,:,iq), qqe+iq )
704  enddo
705 
706 !OCL XFILL
707  flux_rain(:,:,:) = 0.0_rp
708 !OCL XFILL
709  flux_snow(:,:,:) = 0.0_rp
710 
711  ! tracer/energy transport by falldown
712  ! 1st order upwind, forward euler, velocity is always negative
713 
714  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
715  do j = js, je
716  do i = is, ie
717  rfdz(ks-1,i,j) = 1.0_rp / ( real_cz(ks,i,j) - real_fz(ks-1,i,j) )
718  do k = ks, ke
719  rcdz(k,i,j) = 1.0_rp / ( real_fz(k,i,j) - real_fz(k-1,i,j) )
720  rcdz_u(k,i,j) = 2.0_rp / ( ( real_fz(k,i+1,j) - real_fz(k-1,i+1,j) ) &
721  + ( real_fz(k,i ,j) - real_fz(k-1,i ,j) ) )
722  rcdz_v(k,i,j) = 2.0_rp / ( ( real_fz(k,i,j+1) - real_fz(k-1,i,j+1) ) &
723  + ( real_fz(k,i,j ) - real_fz(k-1,i,j ) ) )
724  rfdz(k,i,j) = 1.0_rp / ( real_cz(k+1,i,j) - real_cz(k,i,j) )
725  enddo
726  enddo
727  enddo
728 
729  do iq = i_qc, qqe
730 
731  if ( tracer_type /= 'SUZUKI10' ) then
732  call comm_wait( vterm(:,:,:,iq), iq )
733  endif
734  call comm_wait( qtrc(:,:,:,iq), qqe+iq )
735 
736  do j = js, je
737  do i = is, ie
738  !--- mass flux for each mass tracer, upwind with vel < 0
739  do k = ks-1, ke-1
740  qflx(k,i,j,iq) = vterm(k+1,i,j,iq) * dens(k+1,i,j) * qtrc(k+1,i,j,iq) * j33g
741  enddo
742  qflx(ke,i,j,iq) = 0.0_rp
743 
744  !--- internal energy
745  eflx(ks-1,i,j) = qflx(ks-1,i,j,iq) * temp(ks,i,j) * cvw(iq)
746  do k = ks, ke-1
747  eflx(k,i,j) = qflx(k,i,j,iq) * temp(k+1,i,j) * cvw(iq)
748  rhoe(k,i,j) = rhoe(k,i,j) - dt * ( eflx(k,i,j) - eflx(k-1,i,j) ) * rcdz(k,i,j)
749  enddo
750  eflx(ke,i,j) = 0.0_rp
751  rhoe(ke,i,j) = rhoe(ke,i,j) - dt * ( - eflx(ke-1,i,j) ) * rcdz(ke,i,j)
752 
753  !--- potential energy
754  eflx(ks-1,i,j) = qflx(ks-1,i,j,iq) * grav / rfdz(ks-1,i,j)
755  do k = ks, ke-1
756  eflx(k,i,j) = qflx(k,i,j,iq) * grav / rfdz(k,i,j)
757  rhoe(k,i,j) = rhoe(k,i,j) - dt * ( eflx(k,i,j) - eflx(k-1,i,j) ) * rcdz(k,i,j)
758  enddo
759  rhoe(ke,i,j) = rhoe(ke,i,j) - dt * ( - eflx(ke-1,i,j) ) * rcdz(ke,i,j)
760 
761  !--- momentum z (half level)
762  do k = ks-1, ke-1
763  eflx(k,i,j) = 0.25_rp * ( vterm(k+1,i,j,iq) + vterm(k,i,j,iq) ) &
764  * ( qtrc(k+1,i,j,iq) + qtrc(k,i,j,iq) ) &
765  * momz(k,i,j)
766  enddo
767  do k = ks, ke-1
768  momz(k,i,j) = momz(k,i,j) - dt * ( eflx(k+1,i,j) - eflx(k,i,j) ) * rfdz(k,i,j)
769  enddo
770 
771  !--- momentum x
772  eflx(ks-1,i,j) = 0.25_rp * ( vterm(ks,i,j,iq) + vterm(ks,i+1,j,iq) ) &
773  * ( qtrc(ks,i,j,iq) + qtrc(ks,i+1,j,iq) ) &
774  * momx(ks,i,j)
775  do k = ks, ke-1
776  eflx(k,i,j) = 0.25_rp * ( vterm(k+1,i,j,iq) + vterm(k+1,i+1,j,iq) ) &
777  * ( qtrc(k+1,i,j,iq) + qtrc(k+1,i+1,j,iq) ) &
778  * momx(k+1,i,j)
779  momx(k,i,j) = momx(k,i,j) - dt * ( eflx(k,i,j) - eflx(k-1,i,j) ) * rcdz_u(k,i,j)
780  enddo
781  momx(ke,i,j) = momx(ke,i,j) - dt * ( - eflx(ke-1,i,j) ) * rcdz_u(ke,i,j)
782 
783  !--- momentum y
784  eflx(ks-1,i,j) = 0.25_rp * ( vterm(ks,i,j,iq) + vterm(ks,i,j+1,iq) ) &
785  * ( qtrc(ks,i,j,iq) + qtrc(ks,i,j+1,iq) ) &
786  * momy(ks,i,j)
787  do k = ks, ke-1
788  eflx(k,i,j) = 0.25_rp * ( vterm(k+1,i,j,iq) + vterm(k+1,i,j+1,iq) ) &
789  * ( qtrc(k+1,i,j,iq) + qtrc(k+1,i,j+1,iq) ) &
790  * momy(k+1,i,j)
791  momy(k,i,j) = momy(k,i,j) - dt * ( eflx(k,i,j) - eflx(k-1,i,j) ) * rcdz_v(k,i,j)
792  enddo
793  momy(ke,i,j) = momy(ke,i,j) - dt * ( - eflx(ke-1,i,j) ) * rcdz_v(ke,i,j)
794  enddo
795  enddo
796 
797  enddo ! mass tracer loop
798 
799  ! save previous value
800  do iq = 1, qa
801  do j = js, je
802  do i = is, ie
803  rhoq(ks-1,i,j,iq) = dens(ks,i,j) * qtrc(ks,i,j,iq)
804  do k = ks, ke
805  rhoq(k,i,j,iq) = dens(k,i,j) * qtrc(k,i,j,iq)
806  enddo
807  enddo
808  enddo
809  enddo
810 
811  !--- mass flux for each tracer, upwind with vel < 0
812  do iq = i_qc, qa
813 
814  do j = js, je
815  do i = is, ie
816  do k = ks-1, ke-1
817  qflx(k,i,j,iq) = vterm(k+1,i,j,iq) * rhoq(k+1,i,j,iq)
818  enddo
819  qflx(ke,i,j,iq) = 0.0_rp
820  enddo
821  enddo
822  enddo
823 
824  !--- update total density
825  do iq = i_qc, qqe
826  do j = js, je
827  do i = is, ie
828  do k = ks, ke
829  dens(k,i,j) = dens(k,i,j) - dt * ( qflx(k,i,j,iq) - qflx(k-1,i,j,iq) ) * rcdz(k,i,j)
830  enddo
831  enddo
832  enddo
833  enddo ! mass tracer loop
834 
835  !--- update falling tracer (use updated total density)
836  do iq = i_qc, qa
837  do j = js, je
838  do i = is, ie
839  do k = ks, ke
840  qtrc(k,i,j,iq) = ( rhoq(k,i,j,iq) - dt * ( qflx(k,i,j,iq) - qflx(k-1,i,j,iq) ) * rcdz(k,i,j) ) / dens(k,i,j)
841  enddo
842  enddo
843  enddo
844  enddo
845 
846  !--- update no-falling tracer (use updated total density)
847  do j = js, je
848  do i = is, ie
849  do k = ks, ke
850  qtrc(k,i,j,i_qv) = rhoq(k,i,j,i_qv) / dens(k,i,j)
851  enddo
852  enddo
853  enddo
854 
855  !--- lowermost flux is saved for land process
856  if ( qws > 0 ) then
857  do j = js, je
858  do i = is, ie
859  do k = ks-1, ke
860  do iq = qws, qwe
861  flux_rain(k,i,j) = flux_rain(k,i,j) - qflx(k,i,j,iq)
862  enddo
863  enddo
864  enddo
865  enddo
866  endif
867  if ( qis > 0 ) then
868  do j = js, je
869  do i = is, ie
870  do k = ks-1, ke
871  do iq = qis, qie
872  flux_snow(k,i,j) = flux_snow(k,i,j) - qflx(k,i,j,iq)
873  enddo
874  enddo
875  enddo
876  enddo
877  endif
878 
879 
880  call prof_rapend ('MP_Precipitation', 2)
881 
882  return
883  end subroutine atmos_phy_mp_precipitation
884 
885 end module scale_atmos_phy_mp_common
scale_tracer::qie
integer, public qie
Definition: scale_tracer.F90:66
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_grid_index::jeb
integer, public jeb
Definition: scale_grid_index.F90:74
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation
subroutine, public atmos_phy_mp_precipitation(flux_rain, flux_snow, DENS, MOMZ, MOMX, MOMY, RHOE, QTRC, vterm, temp, dt)
precipitation transport
Definition: scale_atmos_phy_mp_common.F90:657
scale_atmos_saturation::cvovr_liq
real(rp), public cvovr_liq
Definition: scale_atmos_sub_saturation.F90:139
scale_time::time_dtsec_atmos_phy_mp
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
scale_tracer::qqe
integer, public qqe
Definition: scale_tracer.F90:61
scale_tracer::qwe
integer, public qwe
Definition: scale_tracer.F90:64
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_atmos_saturation::cvovr_ice
real(rp), public cvovr_ice
Definition: scale_atmos_sub_saturation.F90:140
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_gridtrans::gtrans_j33g
real(rp), public gtrans_j33g
(3,3) element of Jacobian matrix * {G}^1/2
Definition: scale_gridtrans.F90:43
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_real::real_fz
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
Definition: scale_grid_real.F90:41
scale_tracer::qws
integer, public qws
Definition: scale_tracer.F90:63
scale_grid_real::real_cz
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
Definition: scale_grid_real.F90:40
scale_tracer::qis
integer, public qis
Definition: scale_tracer.F90:65
scale_tracer::tracer_type
character(len=h_short), public tracer_type
Definition: scale_tracer.F90:109
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_grid_index::ieb
integer, public ieb
Definition: scale_grid_index.F90:72
scale_atmos_phy_mp_common
module ATMOSPHERE / Physics Cloud Microphysics - Common
Definition: scale_atmos_phy_mp_common.F90:16
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_tracer::i_qv
integer, public i_qv
Definition: scale_tracer.F90:47
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_time
module TIME
Definition: scale_time.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_atmos_saturation::lovr_ice
real(rp), public lovr_ice
Definition: scale_atmos_sub_saturation.F90:142
scale_const::const_lhv
real(rp), public const_lhv
latent heat of vaporizaion for use
Definition: scale_const.F90:77
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_atmos_thermodyn::aq_cv
real(rp), dimension(:), allocatable, public aq_cv
CV for each hydrometeors [J/kg/K].
Definition: scale_atmos_sub_thermodyn.F90:143
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer
subroutine, public atmos_phy_mp_negative_fixer(DENS, RHOT, QTRC)
Negative fixer.
Definition: scale_atmos_phy_mp_common.F90:64
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:132
scale_const::const_lhf
real(rp), public const_lhf
latent heat of fusion for use
Definition: scale_const.F90:79
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_tracer::i_qi
integer, public i_qi
Definition: scale_tracer.F90:50
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_atmos_saturation::lovr_liq
real(rp), public lovr_liq
Definition: scale_atmos_sub_saturation.F90:141
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_tracer::qqs
integer, public qqs
Definition: scale_tracer.F90:60
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_grid_index::isb
integer, public isb
Definition: scale_grid_index.F90:71
scale_atmos_phy_mp_common::atmos_phy_mp_saturation_adjustment
subroutine, public atmos_phy_mp_saturation_adjustment(RHOE_t, QTRC_t, RHOE0, QTRC0, DENS0, flag_liquid)
Saturation adjustment.
Definition: scale_atmos_phy_mp_common.F90:129
scale_grid_index::jsb
integer, public jsb
Definition: scale_grid_index.F90:73
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:178
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_tracer::i_qc
integer, public i_qc
Definition: scale_tracer.F90:48
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13