SCALE-RM
Functions/Subroutines | Variables
scale_atmos_phy_mp_suzuki10 Module Reference

Module Spectran Bin Microphysics. More...

Functions/Subroutines

subroutine, public atmos_phy_mp_suzuki10_setup (MP_TYPE)
 Setup. More...
 
subroutine, public atmos_phy_mp_suzuki10 (DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
 Cloud Microphysics. More...
 
subroutine mp_suzuki10 (ijkmax, ijkmax_cold, ijkmax_warm, index_cold, index_warm, dens, pres, ccn, temp, qvap, ghyd, gaer, evaporate, dt)
 
subroutine r_collcoag (ijkmax, swgt, temp, gc, dtime)
 
subroutine, public atmos_phy_mp_suzuki10_cloudfraction (cldfrac, QTRC, mask_criterion)
 Calculate Cloud Fraction. More...
 
subroutine, public atmos_phy_mp_suzuki10_effectiveradius (Re, QTRC0, DENS0, TEMP0)
 Calculate Effective Radius. More...
 
subroutine, public atmos_phy_mp_suzuki10_mixingratio (Qe, QTRC0)
 Calculate mixing ratio of each category. More...
 
subroutine mkpara
 

Variables

real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
 

Detailed Description

Module Spectran Bin Microphysics.

Description:
This module contains subroutines for the Spectral Bin Model
Author
: Team SCALE
History: scale_atmos_phy_mp_suzuki10
  • ver.0.00 2012-06-14 (Y.Sato) [new] Import from version 4.1 of original code
  • ver.0.01 2012-09-14 (Y.Sato) [mod] add a stochastic method (Sato et al. 2009)
  • ver.0.02 2013-02-12 (Y.Sato) [mod] modified for latest version
  • ver.0.03 2013-12-26 (Y.Sato) [mod] mearge all version of Bin scheme
  • ver.0.04 2015-09-02 (Y.Sato) [mod] Tuning for K and FX10
  • ver.0.05 2015-09-08 (Y.Sato) [mod] Add evaporated cloud number concentration
  • ver.0.06 2016-01-22 (Y.Sato) [mod] Modify several bugs for using with MSTRNX
NAMELIST
  • PARAM_ATMOS_PHY_MP
    nametypedefault valuecomment
    MP_NTMAX_SEDIMENTATION integer 1 maxinum fractional step
    MP_DOAUTOCONVERSION logical .true. apply collision process ?
    MP_DOPRECIPITATION logical .true. apply sedimentation of hydrometeor ?
    MP_DONEGATIVE_FIXER logical .true. apply negative fixer?
    MP_COUPLE_AEROSOL logical .false. apply CCN effect?

  • PARAM_ATMOS_PHY_MP_SUZUKI10
    nametypedefault valuecomment
    RHO_AERO real(RP) — density of aerosol
    R_MIN real(RP) — minimum radius of aerosol (um)
    R_MAX real(RP) — maximum radius of aerosol (um)
    R0_AERO real(RP) — center radius of aerosol (um)
    S10_EMAER real(RP) — moleculer weight of aerosol
    S10_FLAG_REGENE logical — flag of regeneration
    S10_FLAG_NUCLEAT logical — flag of regeneration
    S10_FLAG_ICENUCLEAT logical — flag of regeneration
    S10_FLAG_SFAERO logical — flag of surface flux of aeorol
    S10_RNDM_FLGP integer — flag of surface flux of aeorol
    S10_RNDM_MSPC integer
    S10_RNDM_MBIN integer
    C_CCN real(RP) 100.E+6_RP N0 of Nc = N0*s^kappa
    KAPPA real(RP) 0.462_RP kappa of Nc = N0*s^kappa
    SIGMA real(RP) 7.5E-02_RP water surface tension N/m2
    VHFCT real(RP) 2.0_RP van't hoff factor (i in eq.(A.11) of Suzuki (2004))

History Output
namedescriptionunitvariable
QC Mixing ratio of QC kg/kg QHYD_out
QG Mixing ratio of QG kg/kg QHYD_out
QH Mixing ratio of QH kg/kg QHYD_out
QI Mixing ratio of QI kg/kg QHYD_out
QR Mixing ratio of QR kg/kg QHYD_out
QS Mixing ratio of QS kg/kg QHYD_out

Function/Subroutine Documentation

◆ atmos_phy_mp_suzuki10_setup()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_setup ( character(len=*), intent(in)  MP_TYPE)

Setup.

Definition at line 239 of file scale_atmos_phy_mp_suzuki10.F90.

References atmos_phy_mp_dens, scale_comm::comm_datatype, scale_comm::comm_world, scale_const::const_dice, scale_const::const_dwatr, scale_grid::grid_cdz, scale_grid::grid_cz, scale_grid_index::ia, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_get_available_fid(), scale_stdio::io_l, scale_grid_index::ja, scale_grid_index::ka, dc_log::log(), mkpara(), scale_process::prc_ismaster, scale_process::prc_masterrank, scale_process::prc_mpistop(), scale_tracer::qa, and scale_time::time_dtsec_atmos_phy_mp.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

239  use scale_process, only: &
240  prc_mpistop, &
241  prc_masterrank, &
243  use scale_const, only: &
244  const_dwatr, &
245  const_dice
246  use scale_comm, only: &
247  comm_world, &
249  use scale_grid, only: &
250  cdz => grid_cdz, &
251  cz => grid_cz
252  use scale_time, only: &
254  use scale_tracer, only: &
255  qad => qa
256  implicit none
257 
258  character(len=*), intent(in) :: mp_type
259 
260  real(RP) :: rho_aero !--- density of aerosol
261  real(RP) :: r0_aero !--- center radius of aerosol (um)
262  real(RP) :: r_min !--- minimum radius of aerosol (um)
263  real(RP) :: r_max !--- maximum radius of aerosol (um)
264  real(RP) :: s10_emaer !--- moleculer weight of aerosol
265  logical :: s10_flag_regene !--- flag of regeneration
266  logical :: s10_flag_nucleat !--- flag of regeneration
267  logical :: s10_flag_icenucleat !--- flag of regeneration
268  logical :: s10_flag_sfaero !--- flag of surface flux of aeorol
269  integer :: s10_rndm_flgp !--- flag of surface flux of aeorol
270  integer :: s10_rndm_mspc
271  integer :: s10_rndm_mbin
272 
273  namelist / param_atmos_phy_mp / &
274  mp_ntmax_sedimentation, &
275  mp_doautoconversion, &
276  mp_doprecipitation, &
277  mp_donegative_fixer, &
278  mp_couple_aerosol
279 
280  namelist / param_atmos_phy_mp_suzuki10 / &
281  rho_aero, &
282  r_min, &
283  r_max, &
284  r0_aero, &
285  s10_emaer, &
286  s10_flag_regene, &
287  s10_flag_nucleat, &
288  s10_flag_icenucleat, &
289  s10_flag_sfaero, &
290  s10_rndm_flgp, &
291  s10_rndm_mspc, &
292  s10_rndm_mbin, &
293  c_ccn, kappa, &
294  sigma, vhfct
295 
296  real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
297  integer :: nstep_max
298  integer :: nnspc, nnbin
299  integer :: nn, mm, mmyu, nnyu
300  integer :: myu, nyu, i, j, k, n, ierr
301  !---------------------------------------------------------------------------
302 
303  !--- allocation
304  allocate( xctr( nbin ) )
305  allocate( xbnd( nbin+1 ) )
306  allocate( radc( nbin ) )
307  allocate( cctr( nbin,nspc_mk ) )
308  allocate( cbnd( nbin+1,nspc_mk ) )
309  allocate( ck( nspc_mk,nspc_mk,nbin,nbin ) )
310  allocate( vt( nspc_mk,nbin ) )
311  allocate( br( nspc_mk,nbin ) )
312  allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
313  allocate( expxctr( nbin ) )
314  allocate( expxbnd( nbin+1 ) )
315  allocate( rexpxctr( nbin ) )
316  allocate( rexpxbnd( nbin+1 ) )
317  if ( nccn /= 0 ) then
318  allocate( xactr( nccn ) )
319  allocate( xabnd( nccn+1 ) )
320  allocate( rada( nccn ) )
321  allocate( expxactr( nccn ) )
322  allocate( expxabnd( nccn+1 ) )
323  allocate( rexpxactr( nccn ) )
324  allocate( rexpxabnd( nccn+1 ) )
325  endif
326 
327  mbin = nbin/2
328  mspc = nspc_mk*nspc_mk
329 
330  if( io_l ) write(io_fid_log,*)
331  if( io_l ) write(io_fid_log,*) '+++ Module[Cloud Microphisics]/Categ[ATMOS]'
332  if( io_l ) write(io_fid_log,*) '*** Wrapper for SBM (warm cloud)'
333 
334  if ( mp_type /= 'SUZUKI10' ) then
335  write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not SUZUKI10. Check!'
336  call prc_mpistop
337  endif
338 
339  rho_aero = rhoa
340  s10_emaer = emaer
341  r_min = rasta
342  r_max = raend
343  r0_aero = r0a
344  s10_flag_regene = flg_regeneration
345  s10_flag_nucleat = flg_nucl
346  s10_flag_icenucleat = flg_icenucl
347  s10_flag_sfaero = flg_sf_aero
348  s10_rndm_flgp = rndm_flgp
349  s10_rndm_mspc = mspc
350  s10_rndm_mbin = mbin
351 
352  rewind(io_fid_conf)
353  read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
354 
355  if ( ierr < 0 ) then !--- missing
356  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
357  elseif( ierr > 0 ) then !--- fatal error
358  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP, Check!'
359  call prc_mpistop
360  endif
361  if( io_l ) write(io_fid_log,nml=param_atmos_phy_mp)
362 
363  rewind(io_fid_conf)
364  read(io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
365 
366  if ( ierr < 0 ) then !--- missing
367  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
368  elseif( ierr > 0 ) then !--- fatal error
369  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!'
370  call prc_mpistop
371  endif
372  if( io_l ) write(io_fid_log,nml=param_atmos_phy_mp_suzuki10)
373 
374  if ( nspc /= 1 .AND. nspc /= 7 ) then
375  write(*,*) 'xxx nspc should be set as 1(warm rain) or 7(mixed phase) check!'
376  call prc_mpistop
377  endif
378 
379  rhoa = rho_aero
380  emaer = s10_emaer
381  rasta = r_min
382  raend = r_max
383  r0a = r0_aero
384  flg_regeneration = s10_flag_regene
385  flg_nucl = s10_flag_nucleat
386  flg_icenucl = s10_flag_icenucleat
387  flg_sf_aero = s10_flag_sfaero
388  rndm_flgp = s10_rndm_flgp
389  mspc = s10_rndm_mspc
390  mbin = s10_rndm_mbin
391 
392  !--- read micpara.dat (microphysical parameter) and broad cast
393  if ( prc_ismaster ) then
394 
395  fid_micpara = io_get_available_fid()
396  !--- open parameter of cloud microphysics
397  open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
398 
399  !--- micpara.dat does not exist
400  if ( ierr == 0 ) then
401 
402  read( fid_micpara,* ) nnspc, nnbin
403 
404  if ( nnbin /= nbin ) then
405  write(*,*) 'xxx nbin in inc_tracer and nbin in micpara.dat is different check!'
406  call prc_mpistop
407  endif
408 
409  ! grid parameter
410  if( io_l ) write(io_fid_log,*) '*** Radius of cloud ****'
411  do n = 1, nbin
412  read( fid_micpara,* ) nn, xctr( n ), radc( n )
413  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
414  "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
415  enddo
416  do n = 1, nbin+1
417  read( fid_micpara,* ) nn, xbnd( n )
418  enddo
419  read( fid_micpara,* ) dxmic
420  if( io_l ) write(io_fid_log,*) '*** Width of Cloud SDF= ', dxmic
421 
422  ! capacity
423  do myu = 1, nspc_mk
424  do n = 1, nbin
425 ! read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
426  read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
427  enddo
428  do n = 1, nbin+1
429 ! read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
430  read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
431  enddo
432  enddo
433 
434  ! collection kernel
435  do myu = 1, nspc_mk
436  do nyu = 1, nspc_mk
437  do i = 1, nbin
438  do j = 1, nbin
439  read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
440  enddo
441  enddo
442  enddo
443  enddo
444 
445  ! terminal velocity
446  do myu = 1, nspc_mk
447  do n = 1, nbin
448  read( fid_micpara,* ) mmyu, nn, vt( myu,n )
449  enddo
450  enddo
451 
452  ! bulk density
453  do myu = 1, nspc_mk
454  do n = 1, nbin
455  read( fid_micpara,* ) mmyu, nn, br( myu,n )
456  enddo
457  enddo
458 
459  close ( fid_micpara )
460 
461  !--- micpara.dat does not exist
462  else
463 
464  if ( io_l ) write(io_fid_log,*) 'micpara.dat is created'
465  call mkpara
466 
467  fid_micpara = io_get_available_fid()
468  !--- open parameter of cloud microphysics
469  open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
470 
471  read( fid_micpara,* ) nnspc, nnbin
472 
473  if ( nnbin /= nbin ) then
474  write(*,*) 'xxx nbin in inc_tracer and nbin in micpara.dat is different check!'
475  call prc_mpistop
476  endif
477 
478  ! grid parameter
479  if( io_l ) write(io_fid_log,*) '*** Radius of cloud ****'
480  do n = 1, nbin
481  read( fid_micpara,* ) nn, xctr( n ), radc( n )
482  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
483  "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
484  enddo
485  do n = 1, nbin+1
486  read( fid_micpara,* ) nn, xbnd( n )
487  enddo
488  read( fid_micpara,* ) dxmic
489  if( io_l ) write(io_fid_log,*) '*** Width of Cloud SDF= ', dxmic
490 
491  ! capacity
492  do myu = 1, nspc_mk
493  do n = 1, nbin
494 ! read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
495  read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
496  enddo
497  do n = 1, nbin+1
498 ! read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
499  read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
500  enddo
501  enddo
502 
503  ! collection kernel
504  do myu = 1, nspc_mk
505  do nyu = 1, nspc_mk
506  do i = 1, nbin
507  do j = 1, nbin
508  read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
509  enddo
510  enddo
511  enddo
512  enddo
513 
514  ! terminal velocity
515  do myu = 1, nspc_mk
516  do n = 1, nbin
517  read( fid_micpara,* ) mmyu, nn, vt( myu,n )
518  enddo
519  enddo
520 
521  ! bulk density
522  do myu = 1, nspc_mk
523  do n = 1, nbin
524  read( fid_micpara,* ) mmyu, nn, br( myu,n )
525  enddo
526  enddo
527 
528  close ( fid_micpara )
529 
530  endif
531 
532  endif
533 
534  call mpi_bcast( radc, nbin, comm_datatype, prc_masterrank, comm_world, ierr )
535  call mpi_bcast( xctr, nbin, comm_datatype, prc_masterrank, comm_world, ierr )
536  call mpi_bcast( dxmic,1, comm_datatype, prc_masterrank, comm_world, ierr )
537  call mpi_bcast( xbnd, nbin+1, comm_datatype, prc_masterrank, comm_world, ierr )
538  call mpi_bcast( cctr, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
539  call mpi_bcast( cbnd, (nbin+1)*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
540  call mpi_bcast( ck, nspc_mk*nspc_mk*nbin*nbin, comm_datatype, prc_masterrank, comm_world, ierr )
541  call mpi_bcast( br, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
542  call mpi_bcast( vt, nbin*nspc_mk, comm_datatype, prc_masterrank, comm_world, ierr )
543 
544  !--- aerosol ( CCN ) (not supported)
545  if ( nccn /= 0 ) then
546 
547  allocate ( ncld( 1:nccn ) )
548  xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
549  xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
550 
551  dxaer = ( xaend-xasta )/nccn
552 
553  do n = 1, nccn+1
554  xabnd( n ) = xasta + dxaer*( n-1 )
555  enddo
556  do n = 1, nccn
557  xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
558  rada( n ) = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
559  if( io_l ) write(io_fid_log,'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
560  "Radius of ", n, "th aerosol bin (bin center)= ", rada( n ) , "[m]"
561  enddo
562 
563  if ( flg_sf_aero ) then
564  write( *,* ) "flg_sf_aero=true is not supported stop!! "
565  call prc_mpistop
566 ! if ( CZ(KS) >= 10.0_RP ) then
567 ! R10M1 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
568 ! R10M2 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
569 ! R10H1 = 1.0_RP * 0.50_RP
570 ! R10H2 = 1.0_RP * 0.0_RP
571 ! R10E1 = 1.0_RP * 0.50_RP
572 ! R10E2 = 1.0_RP * 0.50_RP
573 ! K10_1 = KS
574 ! K10_2 = KS
575 ! else
576 ! k = 1
577 ! do while ( CZ(k) < 10.0_RP )
578 ! k = k + 1
579 ! K10_1 = k
580 ! K10_2 = k + 1
581 ! R10M1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
582 ! R10M2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
583 ! R10H1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
584 ! R10H2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
585 ! R10E1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
586 ! R10E2 = ( 10.0_RP - CZ(k) ) / CDZ(k)
587 ! enddo
588 ! endif
589  endif
590 
591  endif
592 
593  atmos_phy_mp_dens(i_mp_qc) = const_dwatr
594  atmos_phy_mp_dens(i_mp_qp) = const_dice
595  atmos_phy_mp_dens(i_mp_qcl) = const_dice
596  atmos_phy_mp_dens(i_mp_qd) = const_dice
597  atmos_phy_mp_dens(i_mp_qs) = const_dice
598  atmos_phy_mp_dens(i_mp_qg) = const_dice
599  atmos_phy_mp_dens(i_mp_qh) = const_dice
600 
601  !--- determine nbnd
602  do n = 1, nbin
603  if( radc( n ) > rbnd ) then
604  nbnd = n
605  exit
606  endif
607  enddo
608  if( io_l ) write(io_fid_log,'(A,ES15.7,A)') '*** Radius between cloud and rain is ', radc(nbnd), '[m]'
609 
610  !--- random number setup for stochastic method
611  if ( rndm_flgp > 0 ) then
612  call random_setup( ia*ja*ka )
613  endif
614 
615  if ( mp_couple_aerosol .AND. nccn /=0 ) then
616  write(*,*) 'xxx nccn should be 0 when MP_couple_aerosol = .true. !! stop'
617  call prc_mpistop
618  endif
619 
620  if ( nccn /= 0 ) then
621  do n = 1, nccn
622  expxactr( n ) = exp( xactr( n ) )
623  rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
624  enddo
625  do n = 1, nccn+1
626  expxabnd( n ) = exp( xabnd( n ) )
627  rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
628  enddo
629  endif
630 
631  allocate( vterm(ka,ia,ja,qad) )
632  vterm(:,:,:,:) = 0.0_rp
633  do myu = 1, nspc
634  do n = 1, nbin
635  vterm(:,:,:,i_qv+(myu-1)*nbin+n) = -vt( myu,n )
636  enddo
637  enddo
638  do n = 1, nbin
639  expxctr( n ) = exp( xctr( n ) )
640  rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
641  enddo
642  do n = 1, nbin+1
643  expxbnd( n ) = exp( xbnd( n ) )
644  rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
645  enddo
646 
647  allocate( kindx(nbin,nbin) )
648  call getrule( ifrsl,kindx )
649 
650  if ( const_thermodyn_type == 'EXACT' ) then
651  flg_thermodyn = 1.0_rp
652  elseif( const_thermodyn_type == 'SIMPLE' ) then
653  flg_thermodyn = 0.0_rp
654  endif
655  rtem00 = 1.0_rp / const_tem00
656 
657  nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
658  mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
659 
660  mp_nstep_sedimentation = mp_ntmax_sedimentation
661  mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
662  mp_dtsec_sedimentation = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
663 
664  if ( io_l ) write(io_fid_log,*)
665  if ( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, ' step'
666  if ( io_l ) write(io_fid_log,*) '*** DT of sedimentation is : ', mp_dtsec_sedimentation, '[s]'
667 
668  return
669 
integer, public comm_datatype
datatype of variable
Definition: scale_comm.F90:117
logical, public prc_ismaster
master process in local communicator?
subroutine, public prc_mpistop
Abort MPI.
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:87
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
real(rp), dimension(:), allocatable, public grid_cz
center coordinate [m]: z, local=global
subroutine, public random_setup
Setup.
integer, public qa
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:88
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer, public ka
of z whole cells (local, with HALO)
real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
integer, public comm_world
communication world ID
Definition: scale_comm.F90:118
module COMMUNICATION
Definition: scale_comm.F90:23
module TIME
Definition: scale_time.F90:15
module PROCESS
subroutine, public log(type, message)
Definition: dc_log.f90:133
module CONSTANT
Definition: scale_const.F90:14
integer, parameter, public prc_masterrank
master process in each communicator
module GRID (cartesian)
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
integer, parameter, public rp
integer, public ja
of y whole cells (local, with HALO)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10 ( real(rp), dimension(ka,ia,ja), intent(inout)  DENS,
real(rp), dimension(ka,ia,ja), intent(inout)  MOMZ,
real(rp), dimension(ka,ia,ja), intent(inout)  MOMX,
real(rp), dimension(ka,ia,ja), intent(inout)  MOMY,
real(rp), dimension(ka,ia,ja), intent(inout)  RHOT,
real(rp), dimension(ka,ia,ja,qad), intent(inout)  QTRC,
real(rp), dimension (ka,ia,ja), intent(in)  CCN,
real(rp), dimension(ka,ia,ja), intent(out)  EVAPORATE,
real(rp), dimension(ia,ja), intent(out)  SFLX_rain,
real(rp), dimension(ia,ja), intent(out)  SFLX_snow 
)

Cloud Microphysics.

Definition at line 685 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer(), scale_atmos_phy_mp_common::atmos_phy_mp_precipitation(), scale_const::const_tem00, scale_grid_index::ie, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, mp_suzuki10(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), scale_tracer::qa, and scale_time::time_dtsec_atmos_phy_mp.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

685  use scale_grid_index
686  use scale_time, only: &
688  use scale_tracer, only: &
689  qad => qa
690  use scale_atmos_thermodyn, only: &
691  thermodyn_rhoe => atmos_thermodyn_rhoe, &
692  thermodyn_rhot => atmos_thermodyn_rhot, &
693  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
694  use scale_atmos_phy_mp_common, only: &
695  mp_negative_fixer => atmos_phy_mp_negative_fixer, &
696  mp_precipitation => atmos_phy_mp_precipitation
697  use scale_history, only: &
698  hist_in
699  implicit none
700 
701  real(RP), intent(inout) :: dens(ka,ia,ja)
702  real(RP), intent(inout) :: momz(ka,ia,ja)
703  real(RP), intent(inout) :: momx(ka,ia,ja)
704  real(RP), intent(inout) :: momy(ka,ia,ja)
705  real(RP), intent(inout) :: rhot(ka,ia,ja)
706  real(RP), intent(inout) :: qtrc(ka,ia,ja,qad)
707  real(RP), intent(in) :: ccn (ka,ia,ja)
708  real(RP), intent(out) :: evaporate(ka,ia,ja) !---- evaporated cloud number concentration [/m3]
709  real(RP), intent(out) :: sflx_rain(ia,ja)
710  real(RP), intent(out) :: sflx_snow(ia,ja)
711 
712  real(RP) :: rhoe (ka,ia,ja)
713  real(RP) :: pott (ka,ia,ja)
714  real(RP) :: temp (ka,ia,ja)
715  real(RP) :: pres (ka,ia,ja)
716  real(RP) :: qsat (ka,ia,ja)
717  real(RP) :: ssliq(ka,ia,ja)
718  real(RP) :: ssice(ka,ia,ja)
719 
720  integer :: ijk_index (kijmax,3)
721  integer :: index_cld (kijmax)
722  integer :: index_cold(kijmax)
723  integer :: index_warm(kijmax)
724  integer :: ijkcount, ijkcount_cold, ijkcount_warm
725  integer :: ijk, indirect
726 
727  real(RP) :: dens_ijk(kijmax)
728  real(RP) :: pres_ijk(kijmax)
729  real(RP) :: temp_ijk(kijmax)
730  real(RP) :: qvap_ijk(kijmax)
731  real(RP) :: ccn_ijk(kijmax)
732  real(RP) :: evaporate_ijk(kijmax)
733  real(RP) :: ghyd_ijk(nbin,nspc,kijmax)
734  real(RP) :: gaer_ijk(nccn1 ,kijmax)
735  real(RP) :: cldsum
736  integer :: countbin
737 
738 ! logical, save :: ofirst_sdfa = .true.
739 ! real(RP) :: VELX (IA,JA)
740 ! real(RP) :: VELY (IA,JA)
741 ! real(RP) :: SFLX_AERO(IA,JA,nccn)
742 ! real(RP) :: Uabs, bparam
743 ! real(RP) :: AMR(KA,IA,JA)
744 
745  real(RP) :: flx_rain (ka,ia,ja)
746  real(RP) :: flx_snow (ka,ia,ja)
747  real(RP) :: wflux_rain(ka,ia,ja)
748  real(RP) :: wflux_snow(ka,ia,ja)
749  real(RP) :: qhyd_out(ka,ia,ja,6)
750  integer :: step
751  integer :: k, i, j, m, n, iq
752  !---------------------------------------------------------------------------
753 
754  if ( nspc == 1 ) then
755  if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(SBM Liquid water only)'
756  elseif( nspc > 1 ) then
757  if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(SBM Mixed phase)'
758  endif
759 
760  if ( mp_donegative_fixer ) then
761  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
762  rhot(:,:,:), & ! [INOUT]
763  qtrc(:,:,:,:) ) ! [INOUT]
764  endif
765 
766  call prof_rapstart('MP_ijkconvert', 3)
767 
768  ! Clear EVAPORATE
769  do j = js, je
770  do i = is, ie
771  do k = ks, ke
772  evaporate(k,i,j) = 0.0_rp
773  enddo
774  enddo
775  enddo
776 
777  ijk = 0
778  do j = js, je
779  do i = is, ie
780  do k = ks, ke
781  ijk = ijk + 1
782  ijk_index(ijk,1) = i
783  ijk_index(ijk,2) = j
784  ijk_index(ijk,3) = k
785  enddo
786  enddo
787  enddo
788 
789  call thermodyn_temp_pres( temp(:,:,:), & ! [OUT]
790  pres(:,:,:), & ! [OUT]
791  dens(:,:,:), & ! [IN]
792  rhot(:,:,:), & ! [IN]
793  qtrc(:,:,:,:) ) ! [IN]
794 
795  call saturation_pres2qsat_liq( qsat(:,:,:), & ! [OUT]
796  temp(:,:,:), & ! [IN]
797  pres(:,:,:) ) ! [IN]
798 
799  do j = js, je
800  do i = is, ie
801  do k = ks, ke
802  ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat(k,i,j) - 1.0_rp
803  enddo
804  enddo
805  enddo
806 
807  if ( nspc == 1 ) then
808  ssice(:,:,:) = 0.0_rp
809  else
810  call saturation_pres2qsat_ice( qsat(:,:,:), & ! [OUT]
811  temp(:,:,:), & ! [IN]
812  pres(:,:,:) ) ! [IN]
813 
814  do j = js, je
815  do i = is, ie
816  do k = ks, ke
817  ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat(k,i,j) - 1.0_rp
818  enddo
819  enddo
820  enddo
821  endif
822 
823 !--- store initial SDF of aerosol
824 !--- this option is not supported
825 ! if ( ofirst_sdfa ) then
826 ! allocate( marate( nccn ) )
827 ! do j = JS, JE
828 ! do i = IS, IE
829 ! do k = KS, KE
830 ! sum2 = 0.0_RP
831 ! do n = 1, nccn
832 ! marate( n ) = gdga(k,i,j,n)*rexpxactr( n )
833 ! sum2 = sum2 + gdga(k,i,j,n)*rexpxactr( n )
834 ! enddo
835 ! enddo
836 ! enddo
837 ! enddo
838 ! if ( sum2 /= 0.0_RP ) then
839 ! marate( 1:nccn ) = marate( 1:nccn )/sum2
840 ! ofirst_sdfa = .false.
841 ! endif
842 ! endif
843 
844  !--- Arrange array for microphysics
845  ijkcount = 0
846  ijkcount_cold = 0
847  ijkcount_warm = 0
848 
849  ijk = 0
850  do j = js, je
851  do i = is, ie
852  do k = ks, ke
853  ijk = ijk + 1
854 
855  cldsum = 0.0_rp
856  countbin = qqs
857  do m = 1, nspc
858  do n = 1, nbin
859  countbin = countbin + 1
860  cldsum = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
861  enddo
862  enddo
863 
864  if ( cldsum > cldmin &
865  .OR. ssliq(k,i,j) > 0.0_rp &
866  .OR. ssice(k,i,j) > 0.0_rp ) then
867 
868  ijkcount = ijkcount + 1
869 
870  index_cld(ijkcount) = ijk
871 
872  dens_ijk(ijkcount) = dens(k,i,j)
873  pres_ijk(ijkcount) = pres(k,i,j)
874  temp_ijk(ijkcount) = temp(k,i,j)
875  ccn_ijk(ijkcount) = ccn(k,i,j)
876  qvap_ijk(ijkcount) = qtrc(k,i,j,i_qv)
877 
878  countbin = qqs
879  do m = 1, nspc
880  do n = 1, nbin
881  countbin = countbin + 1
882  ghyd_ijk(n,m,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
883  enddo
884  enddo
885 
886  do n = 1, nccn
887  countbin = countbin + 1
888  gaer_ijk(n,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
889  enddo
890 
891  if ( temp(k,i,j) < const_tem00 .AND. nspc > 1 ) then ! cold
892  ijkcount_cold = ijkcount_cold + 1
893  index_cold(ijkcount_cold) = ijkcount
894  else ! warm
895  ijkcount_warm = ijkcount_warm + 1
896  index_warm(ijkcount_warm) = ijkcount
897  endif
898  endif
899 
900  enddo
901  enddo
902  enddo
903 
904  call prof_rapend ('MP_ijkconvert', 3)
905 
906  ! tentative timername registration
907  call prof_rapstart('MP_suzuki10', 3)
908  call prof_rapend ('MP_suzuki10', 3)
909  call prof_rapstart('_SBM_Nucleat', 3)
910  call prof_rapend ('_SBM_Nucleat', 3)
911 ! call PROF_rapstart('_SBM_NucleatA', 3)
912 ! call PROF_rapend ('_SBM_NucleatA', 3)
913  call prof_rapstart('_SBM_Liqphase', 3)
914  call prof_rapend ('_SBM_Liqphase', 3)
915  call prof_rapstart('_SBM_Icephase', 3)
916  call prof_rapend ('_SBM_Icephase', 3)
917  call prof_rapstart('_SBM_Mixphase', 3)
918  call prof_rapend ('_SBM_Mixphase', 3)
919  call prof_rapstart('_SBM_AdvLiq', 3)
920  call prof_rapend ('_SBM_AdvLiq', 3)
921  call prof_rapstart('_SBM_AdvIce', 3)
922  call prof_rapend ('_SBM_AdvIce', 3)
923  call prof_rapstart('_SBM_AdvMix', 3)
924  call prof_rapend ('_SBM_AdvMix', 3)
925 ! call PROF_rapstart('_SBM_FAero', 3)
926 ! call PROF_rapend ('_SBM_FAero', 3)
927  call prof_rapstart('_SBM_Freezing', 3)
928  call prof_rapend ('_SBM_Freezing', 3)
929  call prof_rapstart('_SBM_IceNucleat', 3)
930  call prof_rapend ('_SBM_IceNucleat', 3)
931  call prof_rapstart('_SBM_Melting', 3)
932  call prof_rapend ('_SBM_Melting', 3)
933  call prof_rapstart('_SBM_CollCoag', 3)
934  call prof_rapend ('_SBM_CollCoag', 3)
935 ! call PROF_rapstart('_SBM_CollCoagR', 3)
936 ! call PROF_rapend ('_SBM_CollCoagR', 3)
937 
938  if ( ijkcount > 0 ) then
939 
940  call prof_rapstart('MP_suzuki10', 3)
941 
942  call mp_suzuki10( ijkcount, & ! [IN]
943  ijkcount_cold, & ! [IN]
944  ijkcount_warm, & ! [IN]
945  index_cold( 1:ijkcount), & ! [IN]
946  index_warm( 1:ijkcount), & ! [IN]
947  dens_ijk( 1:ijkcount), & ! [IN]
948  pres_ijk( 1:ijkcount), & ! [IN]
949  ccn_ijk( 1:ijkcount), & ! [IN]
950  temp_ijk( 1:ijkcount), & ! [INOUT]
951  qvap_ijk( 1:ijkcount), & ! [INOUT]
952  ghyd_ijk(:,:,1:ijkcount), & ! [INOUT]
953  gaer_ijk(:, 1:ijkcount), & ! [INOUT]
954  evaporate_ijk( 1:ijkcount), & ! [OUT]
955  dt ) ! [IN]
956 
957  call prof_rapend ('MP_suzuki10', 3)
958 
959 ! if ( flg_sf_aero ) then
960 ! do j = JS-2, JE+2
961 ! do i = IS-2, IE+1
962 ! VELX(i,j) = MOMX(K10_1,i,j) / ( DENS(K10_1,i+1,j)+DENS(K10_1,i,j) ) * R10M1 &
963 ! + MOMX(K10_2,i,j) / ( DENS(K10_2,i+1,j)+DENS(K10_2,i,j) ) * R10M2
964 ! enddo
965 ! enddo
966 !
967 ! do j = JS-2, JE+1
968 ! do i = IS-2, IE+2
969 ! VELY(i,j) = MOMY(K10_1,i,j) / ( DENS(K10_1,i,j+1)+DENS(K10_1,i,j) ) * R10M1 &
970 ! + MOMY(K10_2,i,j) / ( DENS(K10_2,i,j+1)+DENS(K10_2,i,j) ) * R10M2
971 ! enddo
972 ! enddo
973 ! endif
974 !
975 ! !--- SURFACE FLUX by Monahan et al. (1986)
976 ! if ( flg_sf_aero .AND. nccn /= 0 ) then
977 ! do j = JS, JE
978 ! do i = IS, IE
979 ! ijk = ( j - JS ) * KMAX * IMAX &
980 ! + ( i - IS ) * KMAX &
981 ! + ( KS - KS ) &
982 ! + 1
983 ! Uabs = sqrt( ( ( VELX(i,j) + VELX(i-1,j ) ) * 0.50_RP )**2 &
984 ! + ( ( VELY(i,j) + VELY(i ,j-1) ) * 0.50_RP )**2 )
985 ! do n = 1, nccn
986 ! if ( rada( n ) <= 2.0E-5_RP .AND. rada( n ) >= 3.0E-7_RP ) then
987 ! bparam = ( 0.38_RP - log( rada( n ) ) )/0.65_RP
988 ! SFLX_AERO(i,j,n) = 1.373_RP * Uabs**( 3.41_RP ) * rada( n )**( -3.0_RP ) &
989 ! * ( 1.0_RP + 0.057_RP * rada( n )**( 1.05_RP ) ) &
990 ! * 10.0_RP**( 1.19_RP * exp( -bparam*bparam ) )
991 ! ! convert from [#/m^2/um/s] -> [kg/m^3/unit log (m)]
992 ! SFLX_AERO(i,j,n) = SFLX_AERO(i,j,n) / DENS(KS,i,j) &
993 ! / GRID_CDZ(KS) * rada( n ) / 3.0_RP * dt * expxactr( n )
994 ! Gaer_ijk(n,ijk) = Gaer_ijk(n,ijk) + SFLX_AERO(i,j,n)/dxaer
995 ! endif
996 ! enddo
997 ! enddo
998 ! enddo
999 ! endif
1000 
1001  call prof_rapstart('MP_ijkconvert', 3)
1002 
1003  !---- return original array
1004  do ijk = 1, ijkcount
1005  indirect = index_cld(ijk)
1006  i = ijk_index(indirect,1)
1007  j = ijk_index(indirect,2)
1008  k = ijk_index(indirect,3)
1009 
1010  temp(k,i,j) = temp_ijk(ijk)
1011  qtrc(k,i,j,i_qv) = qvap_ijk(ijk)
1012  evaporate(k,i,j) = evaporate_ijk(ijk) / dt ! [#/m3/s]
1013 
1014  countbin = qqs
1015  do m = 1, nspc
1016  do n = 1, nbin
1017  countbin = countbin + 1
1018  qtrc(k,i,j,countbin) = ghyd_ijk(n,m,ijk) / dens(k,i,j) * dxmic
1019  enddo
1020  enddo
1021 
1022  do n = 1, nccn
1023  countbin = countbin + 1
1024  qtrc(k,i,j,countbin) = gaer_ijk(n,ijk) / dens(k,i,j) * dxaer
1025  enddo
1026  enddo
1027 
1028  do j = js, je
1029  do i = is, ie
1030  do k = ks, ke
1031  countbin = qqs
1032  do m = 1, nspc
1033  do n = 1, nbin
1034  countbin = countbin + 1
1035  if ( qtrc(k,i,j,countbin) < eps ) then
1036  qtrc(k,i,j,countbin) = 0.0_rp
1037  endif
1038  enddo
1039  enddo
1040  enddo
1041  enddo
1042  enddo
1043 
1044  call thermodyn_pott( pott(:,:,:), & ! [OUT]
1045  temp(:,:,:), & ! [IN]
1046  pres(:,:,:), & ! [IN]
1047  qtrc(:,:,:,:) ) ! [IN]
1048 
1049  do j = js, je
1050  do i = is, ie
1051  do k = ks, ke
1052  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
1053  enddo
1054  enddo
1055  enddo
1056 
1057 ! if ( nccn /= 0 ) then
1058 ! AMR(:,:,:) = 0.0_RP
1059 ! do j = JS, JE
1060 ! do i = IS, IE
1061 ! do k = KS, KE
1062 ! do n = 1, nccn
1063 ! AMR(k,i,j) = AMR(k,i,j) + QTRC(k,i,j,QQE-1+n)
1064 ! enddo
1065 ! enddo
1066 ! enddo
1067 ! enddo
1068 ! endif
1069 
1070  call prof_rapend ('MP_ijkconvert', 3)
1071 
1072  endif
1073 
1074  flx_rain(:,:,:) = 0.0_rp
1075  flx_snow(:,:,:) = 0.0_rp
1076 
1077  if ( mp_doprecipitation ) then
1078 
1079  call thermodyn_rhoe( rhoe(:,:,:), & ! [OUT]
1080  rhot(:,:,:), & ! [IN]
1081  qtrc(:,:,:,:) ) ! [IN]
1082 
1083  do step = 1, mp_nstep_sedimentation
1084 
1085  call thermodyn_temp_pres_e( temp(:,:,:), & ! [OUT]
1086  pres(:,:,:), & ! [OUT]
1087  dens(:,:,:), & ! [IN]
1088  rhoe(:,:,:), & ! [IN]
1089  qtrc(:,:,:,:) ) ! [IN]
1090 
1091  call mp_precipitation( wflux_rain(:,:,:), & ! [OUT]
1092  wflux_snow(:,:,:), & ! [OUT]
1093  dens(:,:,:), & ! [INOUT]
1094  momz(:,:,:), & ! [INOUT]
1095  momx(:,:,:), & ! [INOUT]
1096  momy(:,:,:), & ! [INOUT]
1097  rhoe(:,:,:), & ! [INOUT]
1098  qtrc(:,:,:,:), & ! [INOUT]
1099  vterm(:,:,:,:), & ! [IN]
1100  temp(:,:,:), & ! [IN]
1101  mp_dtsec_sedimentation ) ! [IN]
1102 
1103  do j = js, je
1104  do i = is, ie
1105  do k = ks-1, ke
1106  flx_rain(k,i,j) = flx_rain(k,i,j) + wflux_rain(k,i,j) * mp_rnstep_sedimentation
1107  flx_snow(k,i,j) = flx_snow(k,i,j) + wflux_snow(k,i,j) * mp_rnstep_sedimentation
1108  enddo
1109  enddo
1110  enddo
1111 
1112  enddo
1113 
1114  call thermodyn_rhot( rhot(:,:,:), & ! [OUT]
1115  rhoe(:,:,:), & ! [IN]
1116  qtrc(:,:,:,:) ) ! [IN]
1117  endif
1118 
1119  !##### END MP Main #####
1120 
1121  if ( mp_donegative_fixer ) then
1122  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
1123  rhot(:,:,:), & ! [INOUT]
1124  qtrc(:,:,:,:) ) ! [INOUT]
1125  endif
1126 
1127  qhyd_out(:,:,:,:) = 0.0_rp
1128  do j = js, je
1129  do i = is, ie
1130  do k = ks, ke
1131  do n = 1, nbnd
1132  qhyd_out(k,i,j,1) = qhyd_out(k,i,j,1) + qtrc(k,i,j,qqs+(il-1)*nbin+n)
1133  enddo
1134  do n = nbnd+1, nbin
1135  qhyd_out(k,i,j,2) = qhyd_out(k,i,j,2) + qtrc(k,i,j,qqs+(il-1)*nbin+n)
1136  enddo
1137  enddo
1138  enddo
1139  enddo
1140  if( nspc > 1 ) then
1141  do j = js, je
1142  do i = is, ie
1143  do k = ks, ke
1144  do m = ic, id
1145  do n = 1, nbin
1146  qhyd_out(k,i,j,3) = qhyd_out(k,i,j,3) + ghyd_ijk(n,m,ijk) / dens(k,i,j) * dxmic
1147  enddo
1148  enddo
1149  do n = 1, nbin
1150  qhyd_out(k,i,j,4) = qhyd_out(k,i,j,4) + ghyd_ijk(n,iss,ijk) / dens(k,i,j) * dxmic
1151  enddo
1152  do n = 1, nbin
1153  qhyd_out(k,i,j,5) = qhyd_out(k,i,j,5) + ghyd_ijk(n,ig,ijk) / dens(k,i,j) * dxmic
1154  enddo
1155  do n = 1, nbin
1156  qhyd_out(k,i,j,6) = qhyd_out(k,i,j,6) + ghyd_ijk(n,ih,ijk) / dens(k,i,j) * dxmic
1157  enddo
1158  enddo
1159  enddo
1160  enddo
1161  endif
1162 
1163  if( nspc > 1 ) then
1164 
1165  do j = js, je
1166  do i = is, ie
1167  do k = ks, ke
1168  do m = ic, id
1169  do n = 1, nbin
1170  qhyd_out(k,i,j,3) = qhyd_out(k,i,j,3) + qtrc(k,i,j,qqs+(m-1)*nbin+n)
1171  enddo
1172  enddo
1173  do n = 1, nbin
1174  qhyd_out(k,i,j,4) = qhyd_out(k,i,j,4) + qtrc(k,i,j,qqs+(iss-1)*nbin+n)
1175  enddo
1176  do n = 1, nbin
1177  qhyd_out(k,i,j,5) = qhyd_out(k,i,j,5) + qtrc(k,i,j,qqs+(ig-1)*nbin+n)
1178  enddo
1179  do n = 1, nbin
1180  qhyd_out(k,i,j,6) = qhyd_out(k,i,j,6) + qtrc(k,i,j,qqs+(ih-1)*nbin+n)
1181  enddo
1182  enddo
1183  enddo
1184  enddo
1185 
1186  endif
1187 
1188  call hist_in( qhyd_out(:,:,:,1), 'QC', 'Mixing ratio of QC', 'kg/kg' )
1189  call hist_in( qhyd_out(:,:,:,2), 'QR', 'Mixing ratio of QR', 'kg/kg' )
1190  if( nspc > 1 ) then
1191  call hist_in( qhyd_out(:,:,:,3), 'QI', 'Mixing ratio of QI', 'kg/kg' )
1192  call hist_in( qhyd_out(:,:,:,4), 'QS', 'Mixing ratio of QS', 'kg/kg' )
1193  call hist_in( qhyd_out(:,:,:,5), 'QG', 'Mixing ratio of QG', 'kg/kg' )
1194  call hist_in( qhyd_out(:,:,:,6), 'QH', 'Mixing ratio of QH', 'kg/kg' )
1195  endif
1196 
1197  sflx_rain(:,:) = flx_rain(ks-1,:,:)
1198  sflx_snow(:,:) = flx_snow(ks-1,:,:)
1199 
1200  return
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
subroutine, public atmos_phy_mp_precipitation(flux_rain, flux_snow, DENS, MOMZ, MOMX, MOMY, RHOE, QTRC, vterm, temp, dt)
precipitation transport
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
integer, public ke
end point of inner domain: z, local
integer, public qa
module ATMOSPHERE / Physics Cloud Microphysics - Common
module grid index
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer, public ka
of z whole cells (local, with HALO)
integer, public js
start point of inner domain: y, local
module TIME
Definition: scale_time.F90:15
integer, public ks
start point of inner domain: z, local
subroutine, public atmos_phy_mp_negative_fixer(DENS, RHOT, QTRC)
Negative fixer.
integer, public kijmax
of computational cells: z*x*y
integer, public ie
end point of inner domain: x, local
module ATMOSPHERE / Thermodynamics
module HISTORY
integer, public ja
of y whole cells (local, with HALO)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mp_suzuki10()

subroutine scale_atmos_phy_mp_suzuki10::mp_suzuki10 ( integer, intent(in)  ijkmax,
integer, intent(in)  ijkmax_cold,
integer, intent(in)  ijkmax_warm,
integer, dimension(ijkmax), intent(in)  index_cold,
integer, dimension(ijkmax), intent(in)  index_warm,
real(rp), dimension (ijkmax), intent(in)  dens,
real(rp), dimension (ijkmax), intent(in)  pres,
real(rp), dimension (ijkmax), intent(in)  ccn,
real(rp), dimension (ijkmax), intent(inout)  temp,
real(rp), dimension (ijkmax), intent(inout)  qvap,
real(rp), dimension (nbin,nspc,ijkmax), intent(inout)  ghyd,
real(rp), dimension (nccn1, ijkmax), intent(inout)  gaer,
real(rp), dimension (ijkmax), intent(out)  evaporate,
real(dp), intent(in)  dt 
)

Definition at line 1219 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_const::const_psat0, scale_const::const_tem00, scale_atmos_saturation::cpovr_ice, scale_atmos_saturation::cpovr_liq, dc_log::log(), scale_atmos_saturation::lovr_ice, scale_atmos_saturation::lovr_liq, scale_process::prc_mpistop(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), r_collcoag(), and scale_random::random_get().

Referenced by atmos_phy_mp_suzuki10().

1219  implicit none
1220 
1221  integer, intent(in) :: ijkmax
1222  integer, intent(in) :: ijkmax_cold
1223  integer, intent(in) :: ijkmax_warm
1224  integer , intent(in) :: index_cold(ijkmax)
1225  integer , intent(in) :: index_warm(ijkmax)
1226  real(RP), intent(in) :: dens (ijkmax) ! Density [kg/m3]
1227  real(RP), intent(in) :: pres (ijkmax) ! Pressure [Pa]
1228  real(RP), intent(in) :: ccn (ijkmax) ! Number concentration of CCN [#/m3]
1229  real(RP), intent(inout) :: temp (ijkmax) ! Temperature [K]
1230  real(RP), intent(inout) :: qvap (ijkmax) ! Specific humidity [kg/kg]
1231  real(RP), intent(inout) :: ghyd (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
1232  real(RP), intent(inout) :: gaer (nccn1, ijkmax) ! Mass size distribution function of aerosol
1233  real(RP), intent(out) :: evaporate (ijkmax) ! Number concentration of evaporated cloud [/m3/s]
1234  real(DP), intent(in) :: dt ! Time step interval
1235  !---------------------------------------------------------------------------
1236 
1237  if ( nccn /= 0 ) then
1238  if ( nspc == 1 ) then
1239  !---< warm rain only with aerosol tracer >---
1240 
1241  ! nucleation from aerosol
1242  call nucleata( ijkmax, & ! [IN]
1243  dens(:), & ! [IN]
1244  pres(:), & ! [IN]
1245  temp(:), & ! [INOUT]
1246  qvap(:), & ! [INOUT]
1247  ghyd(:,:,:), & ! [INOUT]
1248  gaer(:,:), & ! [INOUT]
1249  dt ) ! [IN]
1250 
1251  ! condensation / evaporation
1252  call cndevpsbla( ijkmax, & ! [IN]
1253  dens(:), & ! [IN]
1254  pres(:), & ! [IN]
1255  temp(:), & ! [INOUT]
1256  qvap(:), & ! [INOUT]
1257  ghyd(:,:,:), & ! [INOUT]
1258  gaer(:,:), & ! [INOUT]
1259  evaporate(:),& ! [OUT]
1260  dt ) ! [IN]
1261 
1262  if ( mp_doautoconversion ) then
1263  ! collision-coagulation
1264  call collmain( ijkmax, & ! [IN]
1265  temp(:), & ! [IN]
1266  ghyd(:,:,:), & ! [INOUT]
1267  dt ) ! [IN]
1268  endif
1269 
1270  elseif( nspc > 1 ) then
1271  !---< mixed phase rain only with aerosol tracer >---
1272 
1273  ! nucleation from aerosol
1274  call nucleata( ijkmax, & ! [IN]
1275  dens(:), & ! [IN]
1276  pres(:), & ! [IN]
1277  temp(:), & ! [INOUT]
1278  qvap(:), & ! [INOUT]
1279  ghyd(:,:,:), & ! [INOUT]
1280  gaer(:,:), & ! [INOUT]
1281  dt ) ! [IN]
1282 
1283  ! freezing / melting
1284  call freezing( ijkmax, & ! [IN]
1285  ijkmax_cold, & ! [IN]
1286  index_cold(:), & ! [IN]
1287  dens(:), & ! [IN]
1288  temp(:), & ! [INOUT]
1289  ghyd(:,:,:), & ! [INOUT]
1290  dt ) ! [IN]
1291 
1292 ! call ice_nucleat( ijkmax, & ! [IN]
1293 ! ijkmax_cold, & ! [IN]
1294 ! index_cold(:), & ! [IN]
1295 ! dens(:), & ! [IN]
1296 ! pres(:), & ! [IN]
1297 ! temp(:), & ! [INOUT]
1298 ! qvap(:), & ! [INOUT]
1299 ! ghyd(:,:,:), & ! [INOUT]
1300 ! dt ) ! [IN]
1301 
1302  call melting( ijkmax, & ! [IN]
1303  ijkmax_warm, & ! [IN]
1304  index_warm(:), & ! [IN]
1305  dens(:), & ! [IN]
1306  temp(:), & ! [INOUT]
1307  ghyd(:,:,:), & ! [INOUT]
1308  dt ) ! [IN]
1309 
1310  ! condensation / evaporation
1311  call cndevpsbla( ijkmax, & ! [IN]
1312  dens(:), & ! [IN]
1313  pres(:), & ! [IN]
1314  temp(:), & ! [INOUT]
1315  qvap(:), & ! [INOUT]
1316  ghyd(:,:,:), & ! [INOUT]
1317  gaer(:,:), & ! [INOUT]
1318  evaporate(:),& ! [OUT]
1319  dt ) ! [IN]
1320 
1321  if ( mp_doautoconversion ) then
1322  ! collision-coagulation
1323  call collmainf( ijkmax, & ! [IN]
1324  temp(:), & ! [IN]
1325  ghyd(:,:,:), & ! [INOUT]
1326  dt ) ! [IN]
1327  endif
1328 
1329  endif
1330 
1331  elseif( nccn == 0 ) then
1332 
1333  if ( nspc == 1 ) then
1334  !---< warm rain only without aerosol tracer >---
1335 
1336  ! nucleation
1337  call nucleat( ijkmax, & ! [IN]
1338  dens(:), & ! [IN]
1339  pres(:), & ! [IN]
1340  ccn(:), & ! [IN]
1341  temp(:), & ! [INOUT]
1342  qvap(:), & ! [INOUT]
1343  ghyd(:,:,:), & ! [INOUT]
1344  dt ) ! [IN]
1345 
1346  ! condensation / evaporation
1347  call cndevpsbl( ijkmax, & ! [IN]
1348  dens(:), & ! [IN]
1349  pres(:), & ! [IN]
1350  temp(:), & ! [INOUT]
1351  qvap(:), & ! [INOUT]
1352  ghyd(:,:,:), & ! [INOUT]
1353  evaporate(:),& ! [OUT]
1354  dt ) ! [IN]
1355 
1356  if ( mp_doautoconversion ) then
1357  ! collision-coagulation
1358  call collmain( ijkmax, & ! [IN]
1359  temp(:), & ! [IN]
1360  ghyd(:,:,:), & ! [INOUT]
1361  dt ) ! [IN]
1362  endif
1363 
1364  elseif( nspc > 1 ) then
1365  !---< mixed phase rain only without aerosol tracer >---
1366 
1367  ! nucleation
1368  call nucleat( ijkmax, & ! [IN]
1369  dens(:), & ! [IN]
1370  pres(:), & ! [IN]
1371  ccn(:), & ! [IN]
1372  temp(:), & ! [INOUT]
1373  qvap(:), & ! [INOUT]
1374  ghyd(:,:,:), & ! [INOUT]
1375  dt ) ! [IN]
1376 
1377  ! freezing / melting
1378  call freezing( ijkmax, & ! [IN]
1379  ijkmax_cold, & ! [IN]
1380  index_cold(:), & ! [IN]
1381  dens(:), & ! [IN]
1382  temp(:), & ! [INOUT]
1383  ghyd(:,:,:), & ! [INOUT]
1384  dt ) ! [IN]
1385 
1386  call ice_nucleat( ijkmax, & ! [IN]
1387  ijkmax_cold, & ! [IN]
1388  index_cold(:), & ! [IN]
1389  dens(:), & ! [IN]
1390  pres(:), & ! [IN]
1391  temp(:), & ! [INOUT]
1392  qvap(:), & ! [INOUT]
1393  ghyd(:,:,:), & ! [INOUT]
1394  dt ) ! [IN]
1395 
1396  call melting( ijkmax, & ! [IN]
1397  ijkmax_warm, & ! [IN]
1398  index_warm(:), & ! [IN]
1399  dens(:), & ! [IN]
1400  temp(:), & ! [INOUT]
1401  ghyd(:,:,:), & ! [INOUT]
1402  dt ) ! [IN]
1403 
1404  ! condensation / evaporation
1405  call cndevpsbl( ijkmax, & ! [IN]
1406  dens(:), & ! [IN]
1407  pres(:), & ! [IN]
1408  temp(:), & ! [INOUT]
1409  qvap(:), & ! [INOUT]
1410  ghyd(:,:,:), & ! [INOUT]
1411  evaporate(:),& ! [OUT]
1412  dt ) ! [IN]
1413 
1414  if ( mp_doautoconversion ) then
1415  ! collision-coagulation
1416  call collmainf( ijkmax, & ! [IN]
1417  temp(:), & ! [IN]
1418  ghyd(:,:,:), & ! [INOUT]
1419  dt ) ! [IN]
1420  endif
1421 
1422  endif
1423 
1424  endif
1425 
1426  return
Here is the call graph for this function:
Here is the caller graph for this function:

◆ r_collcoag()

subroutine scale_atmos_phy_mp_suzuki10::r_collcoag ( integer, intent(in)  ijkmax,
real(rp), intent(in)  swgt,
real(rp), dimension(ijkmax), intent(in)  temp,
real(rp), dimension (nbin,nspc,ijkmax), intent(inout)  gc,
real(dp), intent(in)  dtime 
)

Definition at line 3622 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_grid_index::ia, scale_grid_index::ja, scale_grid_index::ka, dc_log::log(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), and scale_random::random_get().

Referenced by mp_suzuki10().

3622  use scale_random, only: &
3623  random_get
3624  implicit none
3625 
3626  integer, intent(in) :: ijkmax
3627  real(RP), intent(in) :: swgt
3628  real(RP), intent(in) :: temp(ijkmax) ! Temperature [K]
3629  real(RP), intent(inout) :: gc (nbin,nspc,ijkmax) ! Mass size distribution function of hydrometeor
3630  real(DP), intent(in) :: dtime ! Time step interval
3631 
3632  integer :: i, j, k, l
3633  real(RP) :: xi, xj, xnew, dmpi, dmpj, frci, frcj
3634  real(RP) :: gprime, gprimk, wgt, crn, sum, flux
3635  integer, parameter :: ldeg = 2
3636  real(RP), parameter :: dmpmin = 1.e-01_rp, cmin = 1.e-10_rp
3637  real(RP) :: acoef( 0:ldeg )
3638  !
3639  !--- Y.sato added to use code6
3640  integer :: nums( mbin ), numl( mbin )
3641  real(RP), parameter :: gt = 1.0_rp
3642  integer :: s, det
3643  real(RP) :: nbinr, mbinr ! use to weight
3644 ! real(RP) :: beta
3645  real(RP) :: tmpi, tmpj
3646 
3647  integer :: ibnd( ijkmax )
3648  integer :: iflg( nspc,ijkmax )
3649  integer :: iexst( nbin,nspc,ijkmax )
3650  real(RP) :: csum( nspc,ijkmax )
3651  integer :: ijk, nn, mm, pp, qq, myu, n, isml, ilrg, irsl
3652  !---------------------------------------------------------------------------
3653 
3654  call prof_rapstart('_SBM_CollCoagR', 3)
3655 
3656  iflg( :,: ) = 0
3657  iexst( :,:,: ) = 0
3658  csum( :,: ) = 0.0_rp
3659  do ijk = 1, ijkmax
3660  !--- judgement of particle existence
3661  do n = 1, nbin
3662  csum( il,ijk ) = csum( il,ijk ) + gc( n,il,ijk )*dxmic
3663  csum( ic,ijk ) = csum( ic,ijk ) + gc( n,ic,ijk )*dxmic
3664  csum( ip,ijk ) = csum( ip,ijk ) + gc( n,ip,ijk )*dxmic
3665  csum( id,ijk ) = csum( id,ijk ) + gc( n,id,ijk )*dxmic
3666  csum( iss,ijk ) = csum( iss,ijk ) + gc( n,iss,ijk )*dxmic
3667  csum( ig,ijk ) = csum( ig,ijk ) + gc( n,ig,ijk )*dxmic
3668  csum( ih,ijk ) = csum( ih,ijk ) + gc( n,ih,ijk )*dxmic
3669  enddo
3670  if ( csum( il,ijk ) > cldmin ) iflg( il,ijk ) = 1
3671  if ( csum( ic,ijk ) > cldmin ) iflg( ic,ijk ) = 1
3672  if ( csum( ip,ijk ) > cldmin ) iflg( ip,ijk ) = 1
3673  if ( csum( id,ijk ) > cldmin ) iflg( id,ijk ) = 1
3674  if ( csum( iss,ijk ) > cldmin ) iflg( iss,ijk ) = 1
3675  if ( csum( ig,ijk ) > cldmin ) iflg( ig,ijk ) = 1
3676  if ( csum( ih,ijk ) > cldmin ) iflg( ih,ijk ) = 1
3677 
3678  if ( temp(ijk) < tcrit ) then
3679  ibnd(ijk) = 1
3680  else
3681  ibnd(ijk) = 2
3682  endif
3683 
3684  do myu = 1, nspc
3685  do n = 1, nbin
3686  if ( gc( n,myu,ijk ) > cldmin ) then
3687  iexst( n,myu,ijk ) = 1
3688  endif
3689  enddo
3690  enddo
3691 
3692  enddo
3693 
3694 !OCL PARALLEL
3695  do ijk = 1, ijkmax
3696  do isml = 1, nspc
3697  do nn = 1, iflg( isml,ijk )
3698 
3699  do ilrg = 1, nspc
3700  do mm = 1, iflg( ilrg,ijk )
3701  !--- rule of interaction
3702  irsl = ifrsl( ibnd(ijk),isml,ilrg )
3703 
3704  call random_get( rndm )
3705  det = int( rndm(1,1,1)*ia*ja*ka )
3706  nbinr = real( nbin )
3707  mbinr = real( mbin )
3708  nums( 1:mbin ) = bsml( det,1:mbin )
3709  numl( 1:mbin ) = blrg( det,1:mbin )
3710 
3711  do s = 1, mbin
3712  i = nums( s )
3713  j = numl( s )
3714 
3715  do pp = 1, iexst( i,isml,ijk )
3716  do qq = 1, iexst( j,ilrg,ijk )
3717 
3718  k = kindx( i,j )
3719  xi = expxctr( i )
3720  xj = expxctr( j )
3721  xnew = log( xi+xj )
3722 
3723  dmpi = ck( isml,ilrg,i,j )*gc( j,ilrg,ijk )/xj*dxmic*dtime
3724  dmpj = ck( ilrg,isml,i,j )*gc( i,isml,ijk )/xi*dxmic*dtime
3725 
3726  if ( dmpi <= dmpmin ) then
3727  frci = gc( i,isml,ijk )*dmpi
3728  else
3729  frci = gc( i,isml,ijk )*( 1.0_rp-exp( -dmpi ) )
3730  endif
3731 
3732  if ( dmpj <= dmpmin ) then
3733  frcj = gc( j,ilrg,ijk )*dmpj
3734  else
3735  frcj = gc( j,ilrg,ijk )*( 1.0_rp-exp( -dmpj ) )
3736  endif
3737  tmpi = gc( i,isml,ijk )
3738  tmpj = gc( j,ilrg,ijk )
3739 
3740  gc( i,isml,ijk ) = gc( i,isml,ijk )-frci*swgt
3741  gc( j,ilrg,ijk ) = gc( j,ilrg,ijk )-frcj*swgt
3742 
3743  if ( j /= k ) then
3744  gc( j,ilrg,ijk ) = max( gc( j,ilrg,ijk ), 0.0_rp )
3745  endif
3746  gc( i,isml,ijk ) = max( gc( i,isml,ijk ), 0.0_rp )
3747 
3748  frci = tmpi - gc( i,isml,ijk )
3749  frcj = tmpj - gc( j,ilrg,ijk )
3750 
3751  gprime = frci+frcj
3752 
3753  !-----------------------------------------------
3754  !--- Exponential Flux Method (Bott, 2000, JAS)
3755  !-----------------------------------------------
3756  ! if ( gprime <= 0.0_RP ) cycle !large
3757  ! if ( gprime > 0.0_RP .AND. k < nbin ) then
3758  ! gprimk = gc( (irsl-1)*nbin+k ) + gprime
3759  !
3760  ! beta = log( gc( (irsl-1)*nbin+k+1 )/gprimk+1.E-60_RP )
3761  ! crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
3762  !
3763  ! flux = ( gprime/beta )*( exp( beta*0.50_RP ) -exp( beta*( 0.50_RP-crn ) ) )
3764  ! flux = min( gprimk ,gprime )
3765  !
3766  ! gc( (irsl-1)*nbin+k ) = gprimk - flux
3767  ! gc( (irsl-1)*nbin+k+1 ) = gc( (irsl-1)*nbin+k+1 ) + flux
3768  ! endif
3769 
3770  !-----------------------------------------------
3771  !--- Flux Method (Bott, 1998, JAS)
3772  !-----------------------------------------------
3773 ! if ( gprime <= 0.0_RP ) cycle !large
3774  if ( gprime > 0.0_rp .AND. k < nbin ) then
3775 
3776  gprimk = gc( k,irsl,ijk ) + gprime
3777  wgt = gprime / gprimk
3778  crn = ( xnew-xctr( k ) )/( xctr( k+1 )-xctr( k ) )
3779 
3780  acoef( 0 ) = -( gc( k+1,irsl,ijk )-26.0_rp*gprimk+gc( k-1,irsl,ijk ) )/24.0_rp
3781  acoef( 1 ) = ( gc( k+1,irsl,ijk )-gc( k-1,irsl,ijk ) ) *0.5_rp
3782  acoef( 2 ) = ( gc( k+1,irsl,ijk )-2.0_rp*gprimk+gc( k-1,irsl,ijk ) ) *0.50_rp
3783 
3784  sum = 0.0_rp
3785  do l = 0, ldeg
3786  sum = sum + acoef( l )/( l+1 )/2.0_rp**( l+1 ) &
3787  *( 1.0_rp-( 1.0_rp-2.0_rp*crn )**( l+1 ) )
3788  enddo
3789 
3790  flux = wgt*sum
3791  flux = min( max( flux,0.0_rp ),gprime )
3792 
3793  gc( k,irsl,ijk ) = gprimk - flux
3794  gc( k+1,irsl,ijk ) = gc( k+1,irsl,ijk ) + flux
3795  endif
3796 
3797  enddo
3798  enddo
3799 
3800  enddo ! bin
3801  !
3802  enddo
3803  enddo
3804 
3805  enddo
3806  enddo
3807 
3808  enddo
3809 
3810  call prof_rapend ('_SBM_CollCoagR', 3)
3811 
3812  return
integer, public ia
of x whole cells (local, with HALO)
integer, public ka
of z whole cells (local, with HALO)
subroutine, public log(type, message)
Definition: dc_log.f90:133
subroutine, public random_get(var)
Get random number.
module RANDOM
integer, public ja
of y whole cells (local, with HALO)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_cloudfraction()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_cloudfraction ( real(rp), dimension(ka,ia,ja), intent(out)  cldfrac,
real(rp), dimension (ka,ia,ja,qad), intent(in)  QTRC,
real(rp), intent(in)  mask_criterion 
)

Calculate Cloud Fraction.

Definition at line 3821 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::mp_qa, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

3821  use scale_grid_index
3822  use scale_tracer, only: &
3823  qad => qa, &
3824  mp_qad => mp_qa
3825  use scale_tracer_suzuki10, only: &
3826  i_mp_qc
3827  implicit none
3828 
3829  real(RP), intent(out) :: cldfrac(ka,ia,ja)
3830  real(RP), intent(in) :: qtrc (ka,ia,ja,qad)
3831  real(RP), intent(in) :: mask_criterion
3832 
3833  real(RP) :: qhydro
3834  integer :: k, i, j, iq, ihydro
3835  !---------------------------------------------------------------------------
3836 
3837  if( nspc > 1 ) then
3838  do j = js, je
3839  do i = is, ie
3840  do k = ks, ke
3841  qhydro = 0.0_rp
3842  do ihydro = 1, mp_qa
3843  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
3844  qhydro = qhydro + qtrc(k,i,j,iq)
3845  enddo
3846  enddo
3847  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
3848  enddo
3849  enddo
3850  enddo
3851  elseif( nspc == 1 ) then
3852  do j = js, je
3853  do i = is, ie
3854  do k = ks, ke
3855  qhydro = 0.0_rp
3856  do ihydro = 1, i_mp_qc
3857  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
3858  qhydro = qhydro + qtrc(k,i,j,iq)
3859  enddo
3860  enddo
3861  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
3862  enddo
3863  enddo
3864  enddo
3865  endif
3866 
3867  return
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
integer, public ke
end point of inner domain: z, local
integer, public qa
module TRACER / suzuki10
module grid index
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer, public mp_qa
integer, public ka
of z whole cells (local, with HALO)
integer, public js
start point of inner domain: y, local
integer, public ks
start point of inner domain: z, local
integer, public ie
end point of inner domain: x, local
integer, public ja
of y whole cells (local, with HALO)
Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_effectiveradius()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_effectiveradius ( real(rp), dimension (ka,ia,ja,mp_qad), intent(out)  Re,
real(rp), dimension(ka,ia,ja,qad), intent(in)  QTRC0,
real(rp), dimension(ka,ia,ja), intent(in)  DENS0,
real(rp), dimension(ka,ia,ja), intent(in)  TEMP0 
)

Calculate Effective Radius.

Definition at line 3877 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_const::const_eps, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::mp_qa, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

3877  use scale_grid_index
3878  use scale_const, only: &
3879  eps => const_eps
3880  use scale_tracer, only: &
3881  qad => qa, &
3882  mp_qad => mp_qa
3883  use scale_tracer_suzuki10, only: &
3884  i_mp_qc
3885  implicit none
3886 
3887  real(RP), intent(out) :: re (ka,ia,ja,mp_qad) ! effective radius [cm]
3888  real(RP), intent(in) :: qtrc0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
3889  real(RP), intent(in) :: dens0(ka,ia,ja) ! density [kg/m3]
3890  real(RP), intent(in) :: temp0(ka,ia,ja) ! temperature [K]
3891 
3892  real(RP), parameter :: um2cm = 100.0_rp
3893 
3894  real(RP) :: sum2(ka,ia,ja,mp_qa), sum3(ka,ia,ja,mp_qa)
3895  integer :: i, j, k, iq, ihydro
3896  !---------------------------------------------------------------------------
3897 
3898  sum2(:,:,:,:) = 0.0_rp
3899  sum3(:,:,:,:) = 0.0_rp
3900 
3901  if( nspc > 1 ) then
3902  do ihydro = 1, mp_qa
3903  do k = ks, ke
3904  do j = js, je
3905  do i = is, ie
3906  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
3907  sum3(k,i,j,ihydro) = sum3(k,i,j,ihydro) + &
3908  ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
3909  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
3910  * radc( iq-(i_qv+nbin*(ihydro-1)) )**3.0_rp )
3911  sum2(k,i,j,ihydro) = sum2(k,i,j,ihydro) + &
3912  ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
3913  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
3914  * radc( iq-(i_qv+nbin*(ihydro-1)) )**2.0_rp )
3915  enddo
3916  sum2(k,i,j,ihydro) = 0.5_rp + sign(0.5_rp,sum2(k,i,j,ihydro)-eps)
3917  sum3(k,i,j,ihydro) = 0.5_rp + sign(0.5_rp,sum3(k,i,j,ihydro)-eps)
3918 
3919  if ( sum2(k,i,j,ihydro) /= 0.0_rp ) then
3920  re(k,i,j,ihydro) = sum3(k,i,j,ihydro) / sum2(k,i,j,ihydro) * um2cm
3921  else
3922  re(k,i,j,ihydro) = 0.0_rp
3923  endif
3924  enddo
3925  enddo
3926  enddo
3927  enddo
3928  elseif( nspc == 1 ) then
3929  re(:,:,:,:) = 0.0_rp
3930  do ihydro = 1, i_mp_qc
3931  do k = ks, ke
3932  do j = js, je
3933  do i = is, ie
3934  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
3935  sum3(k,i,j,ihydro) = sum3(k,i,j,ihydro) + &
3936  ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
3937  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
3938  * radc( iq-(i_qv+nbin*(ihydro-1)) )**3.0_rp )
3939  sum2(k,i,j,ihydro) = sum2(k,i,j,ihydro) + &
3940  ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
3941  * rexpxctr( iq-(i_qv+nbin*(ihydro-1)) ) & !--- mass -> number
3942  * radc( iq-(i_qv+nbin*(ihydro-1)) )**2.0_rp )
3943  enddo
3944  sum2(k,i,j,ihydro) = 0.5_rp + sign(0.5_rp,sum2(k,i,j,ihydro)-eps)
3945  sum3(k,i,j,ihydro) = 0.5_rp + sign(0.5_rp,sum3(k,i,j,ihydro)-eps)
3946 
3947  if ( sum2(k,i,j,ihydro) /= 0.0_rp ) then
3948  re(k,i,j,ihydro) = sum3(k,i,j,ihydro) / sum2(k,i,j,ihydro) * um2cm
3949  else
3950  re(k,i,j,ihydro) = 0.0_rp
3951  endif
3952  enddo
3953  enddo
3954  enddo
3955  enddo
3956  endif
3957 
3958  return
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
integer, public ke
end point of inner domain: z, local
integer, public qa
module TRACER / suzuki10
module grid index
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer, public mp_qa
integer, public ka
of z whole cells (local, with HALO)
integer, public js
start point of inner domain: y, local
module CONSTANT
Definition: scale_const.F90:14
integer, public ks
start point of inner domain: z, local
integer, public ie
end point of inner domain: x, local
real(rp), public const_eps
small number
Definition: scale_const.F90:36
integer, public ja
of y whole cells (local, with HALO)
Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_mixingratio()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_mixingratio ( real(rp), dimension (ka,ia,ja,mp_qad), intent(out)  Qe,
real(rp), dimension(ka,ia,ja,qad), intent(in)  QTRC0 
)

Calculate mixing ratio of each category.

Definition at line 3966 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_const::const_eps, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::mp_qa, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

3966  use scale_grid_index
3967  use scale_const, only: &
3968  eps => const_eps
3969  use scale_tracer, only: &
3970  qad => qa, &
3971  mp_qad => mp_qa
3972  use scale_tracer_suzuki10, only: &
3973  i_mp_qc
3974  implicit none
3975 
3976  real(RP), intent(out) :: qe (ka,ia,ja,mp_qad) ! mixing ratio of each cateory [kg/kg]
3977  real(RP), intent(in) :: qtrc0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
3978 
3979  real(RP) :: sum2
3980  integer :: i, j, k, iq, ihydro
3981  !---------------------------------------------------------------------------
3982 
3983  qe(:,:,:,:) = 0.0_rp
3984  if( nspc > 1 ) then
3985  do k = ks, ke
3986  do j = js, je
3987  do i = is, ie
3988  do ihydro = 1, mp_qa
3989  sum2 = 0.0_rp
3990  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
3991  sum2 = sum2 + qtrc0(k,i,j,iq)
3992  enddo
3993  qe(k,i,j,ihydro) = sum2
3994  enddo
3995  enddo
3996  enddo
3997  enddo
3998  elseif( nspc == 1 ) then
3999  do k = ks, ke
4000  do j = js, je
4001  do i = is, ie
4002  do ihydro = 1, i_mp_qc
4003  sum2 = 0.0_rp
4004  do iq = i_qv+nbin*(ihydro-1)+1, i_qv+nbin*ihydro
4005  sum2 = sum2 + qtrc0(k,i,j,iq)
4006  enddo
4007  qe(k,i,j,ihydro) = sum2
4008  enddo
4009  enddo
4010  enddo
4011  enddo
4012  endif
4013 
4014  return
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
integer, public ke
end point of inner domain: z, local
integer, public qa
module TRACER / suzuki10
module grid index
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer, public mp_qa
integer, public ka
of z whole cells (local, with HALO)
integer, public js
start point of inner domain: y, local
module CONSTANT
Definition: scale_const.F90:14
integer, public ks
start point of inner domain: z, local
integer, public ie
end point of inner domain: x, local
real(rp), public const_eps
small number
Definition: scale_const.F90:36
integer, public ja
of y whole cells (local, with HALO)
Here is the caller graph for this function:

◆ mkpara()

subroutine scale_atmos_phy_mp_suzuki10::mkpara ( )

Definition at line 4022 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_stdio::io_fid_log, scale_stdio::io_l, and dc_log::log().

Referenced by atmos_phy_mp_suzuki10_setup().

4022 
4023  implicit none
4024 
4025  integer :: i, j
4026  !-----------------------------------------------------------------------------
4027 
4028  allocate( radc_mk( nbin ) )
4029  allocate( xctr_mk( nbin ) )
4030  allocate( xbnd_mk( nbin+1 ) )
4031  allocate( cctr_mk( nspc_mk,nbin ) )
4032  allocate( cbnd_mk( nspc_mk,nbin+1 ) )
4033  allocate( ck_mk( nspc_mk,nspc_mk,nbin,nbin ) )
4034  allocate( vt_mk( nspc_mk,nbin ) )
4035  allocate( br_mk( nspc_mk,nbin ) )
4036 
4037  !--- file reading
4038  call rdkdat
4039 
4040  !--- grid setting
4041  call sdfgrid
4042 
4043  !--- capacity
4044  call getcp
4045 
4046  !--- collection kernel
4047  call getck
4048 
4049  !--- terminal velocity
4050  call getvt
4051 
4052  !--- bulk radius
4053  call getbr
4054 
4055  !--- output
4056  call paraout
4057 
4058  deallocate( radc_mk )
4059  deallocate( xctr_mk )
4060  deallocate( xbnd_mk )
4061  deallocate( cctr_mk )
4062  deallocate( cbnd_mk )
4063  deallocate( ck_mk )
4064  deallocate( vt_mk )
4065  deallocate( br_mk )
4066 
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ atmos_phy_mp_dens

real(rp), dimension(mp_qa), target, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_dens

Definition at line 99 of file scale_atmos_phy_mp_suzuki10.F90.

Referenced by atmos_phy_mp_suzuki10_setup().

99  real(RP), public, target :: atmos_phy_mp_dens(mp_qa) ! hydrometeor density [kg/m3]=[g/L]
integer, public mp_qa
real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens