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