SCALE-RM
Data Types | Functions/Subroutines
scale_atmos_phy_mp_common Module Reference

module ATMOSPHERE / Physics Cloud Microphysics - Common More...

Functions/Subroutines

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. More...
 
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. More...
 
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)
 
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)
 
subroutine, public atmos_phy_mp_precipitation_momentum (KA, KS, KE, DENS, MOMZ, U, V, mflx, RCDZ, RFDZ, MOMZ_t, RHOU_t, RHOV_t)
 

Detailed Description

module ATMOSPHERE / Physics Cloud Microphysics - Common

Description
Common module for Cloud Microphysics Sedimentation/Precipitation and Saturation adjustment
Author
Team SCALE

Function/Subroutine Documentation

◆ atmos_phy_mp_negative_fixer()

subroutine, public scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
integer, intent(in)  QLA,
integer, intent(in)  QIA,
real(rp), intent(in)  limit_negative,
real(rp), dimension (ka,ia,ja), intent(inout)  DENS,
real(rp), dimension (ka,ia,ja), intent(inout)  TEMP,
real(rp), dimension(ka,ia,ja), intent(inout)  CVtot,
real(rp), dimension(ka,ia,ja), intent(inout)  CPtot,
real(rp), dimension (ka,ia,ja), intent(inout)  QV,
real(rp), dimension (ka,ia,ja,qla+qia), intent(inout)  QTRC,
real(rp), dimension (ka,ia,ja), intent(out), optional  RHOH,
real(rp), dimension(ka,ia,ja), intent(out), optional  DENS_diff,
real(rp), dimension(ka,ia,ja), intent(out), optional  ENGI_diff 
)

ATMOS_PHY_MP_negative_fixer negative fixer.

Definition at line 69 of file scale_atmos_phy_mp_common.F90.

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

References scale_const::const_cpdry, scale_const::const_cvdry, scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_vapor, scale_atmos_hydrometeor::cp_water, scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_vapor, scale_atmos_hydrometeor::cv_water, scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_tracer::k, scale_atmos_grid_cartesc_index::ka, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, scale_prc::prc_abort(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_adjustment().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_saturation_adjustment_3d()

subroutine scale_atmos_phy_mp_common::atmos_phy_mp_saturation_adjustment_3d ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension(ka,ia,ja), intent(in)  DENS,
logical, intent(in)  flag_liquid,
real(rp), dimension (ka,ia,ja), intent(inout)  TEMP,
real(rp), dimension (ka,ia,ja), intent(inout)  QV,
real(rp), dimension (ka,ia,ja), intent(inout)  QC,
real(rp), dimension (ka,ia,ja), intent(inout)  QI,
real(rp), dimension(ka,ia,ja), intent(inout)  CPtot,
real(rp), dimension(ka,ia,ja), intent(inout)  CVtot,
real(rp), dimension(ka,ia,ja), intent(out)  RHOE_d 
)

Saturation adjustment.

Definition at line 283 of file scale_atmos_phy_mp_common.F90.

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

References scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_vapor, scale_atmos_hydrometeor::cp_water, scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_vapor, scale_atmos_hydrometeor::cv_water, scale_atmos_hydrometeor::lhf, scale_atmos_hydrometeor::lhv, scale_prc::prc_abort(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Here is the call graph for this function:

◆ atmos_phy_mp_precipitation_upwind()

subroutine, public scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_upwind ( integer, intent(in), value  KA,
integer, intent(in), value  KS,
integer, intent(in), value  KE,
integer, intent(in), value  QHA,
integer, intent(in), value  QLA,
integer, intent(in), value  QIA,
real(rp), dimension (ka), intent(in)  TEMP,
real(rp), dimension(ka,qha), intent(in)  vterm,
real(rp), dimension (ka), intent(in)  FDZ,
real(rp), dimension (ka), intent(in)  RCDZ,
real(rp), intent(in)  dt,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension (ka), intent(inout)  DENS,
real(rp), dimension (ka,qha), intent(inout)  RHOQ,
real(rp), dimension(ka), intent(inout)  CPtot,
real(rp), dimension(ka), intent(inout)  CVtot,
real(rp), dimension (ka), intent(inout)  RHOE,
real(rp), dimension (ka), intent(out)  mflx,
real(rp), dimension (2), intent(out)  sflx,
real(rp), intent(out)  esflx 
)

Definition at line 451 of file scale_atmos_phy_mp_common.F90.

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

References scale_const::const_grav, scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_water, scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_water, scale_tracer::k, scale_atmos_grid_cartesc_index::ke, and scale_atmos_grid_cartesc_index::ks.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

Here is the caller graph for this function:

◆ atmos_phy_mp_precipitation_semilag()

subroutine, public scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_semilag ( integer, intent(in), value  KA,
integer, intent(in), value  KS,
integer, intent(in), value  KE,
integer, intent(in), value  QHA,
integer, intent(in), value  QLA,
integer, intent(in), value  QIA,
real(rp), dimension (ka), intent(in)  TEMP,
real(rp), dimension(ka,qha), intent(in)  vterm,
real(rp), dimension (ka), intent(in)  FZ,
real(rp), dimension (ka), intent(in)  FDZ,
real(rp), dimension (ka), intent(in)  RCDZ,
real(rp), intent(in)  dt,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension (ka), intent(inout)  DENS,
real(rp), dimension (ka,qha), intent(inout)  RHOQ,
real(rp), dimension(ka), intent(inout)  CPtot,
real(rp), dimension(ka), intent(inout)  CVtot,
real(rp), dimension (ka), intent(inout)  RHOE,
real(rp), dimension (ka), intent(out)  mflx,
real(rp), dimension (2), intent(out)  sflx,
real(rp), intent(out)  esflx 
)

Definition at line 641 of file scale_atmos_phy_mp_common.F90.

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

References scale_const::const_grav, scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_water, scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_water, scale_tracer::k, scale_atmos_grid_cartesc_index::ke, and scale_atmos_grid_cartesc_index::ks.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

Here is the caller graph for this function:

◆ atmos_phy_mp_precipitation_momentum()

subroutine, public scale_atmos_phy_mp_common::atmos_phy_mp_precipitation_momentum ( integer, intent(in), value  KA,
integer, intent(in), value  KS,
integer, intent(in), value  KE,
real(rp), dimension (ka), intent(in)  DENS,
real(rp), dimension (ka), intent(in)  MOMZ,
real(rp), dimension (ka), intent(in)  U,
real(rp), dimension (ka), intent(in)  V,
real(rp), dimension (ka), intent(in)  mflx,
real(rp), dimension (ka), intent(in)  RCDZ,
real(rp), dimension (ka), intent(in)  RFDZ,
real(rp), dimension(ka), intent(out)  MOMZ_t,
real(rp), dimension(ka), intent(out)  RHOU_t,
real(rp), dimension(ka), intent(out)  RHOV_t 
)

Definition at line 852 of file scale_atmos_phy_mp_common.F90.

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

References scale_tracer::k, scale_atmos_grid_cartesc_index::ke, and scale_atmos_grid_cartesc_index::ks.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

Here is the caller graph for this function:
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
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_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:149
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:61
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_hydrometeor::lhf
real(rp), public lhf
latent heat of fusion for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:146
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
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