SCALE-RM
mod_mkinit.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
30 !-------------------------------------------------------------------------------
31 module mod_mkinit
32  !-----------------------------------------------------------------------------
33  !
34  !++ used modules
35  !
36  use scale_precision
37  use scale_stdio
38  use scale_prof
40  use scale_tracer
41 
42  use scale_process, only: &
44  use scale_const, only: &
45  pi => const_pi, &
46  grav => const_grav, &
47  pstd => const_pstd, &
48  rdry => const_rdry, &
49  cpdry => const_cpdry, &
50  lhv => const_lhv, &
51  p00 => const_pre00, &
52  i_sw => const_i_sw, &
53  i_lw => const_i_lw
54  use scale_random, only: &
56  use scale_comm, only: &
57  comm_vars8, &
58  comm_wait
59  use scale_grid, only: &
60  grid_cz, &
61  grid_cx, &
62  grid_cy, &
63  grid_fz, &
64  grid_fx, &
65  grid_fy, &
66  grid_cxg
67  use scale_grid_real, only: &
68  real_cz, &
69  real_fz
70  use scale_atmos_profile, only: &
71  profile_isa => atmos_profile_isa
72  use scale_atmos_hydrostatic, only: &
73  hydrostatic_buildrho => atmos_hydrostatic_buildrho, &
74  hydrostatic_buildrho_atmos => atmos_hydrostatic_buildrho_atmos, &
75  hydrostatic_buildrho_bytemp => atmos_hydrostatic_buildrho_bytemp
76  use scale_atmos_saturation, only: &
77  saturation_pres2qsat_all => atmos_saturation_pres2qsat_all, &
78  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq
79  use mod_atmos_vars, only: &
80  dens, &
81  momx, &
82  momy, &
83  momz, &
84  rhot, &
85  qtrc
86  use mod_atmos_phy_ae_vars, only: &
87  ccn => atmos_phy_ae_ccn
88  !-----------------------------------------------------------------------------
89  implicit none
90  private
91  !-----------------------------------------------------------------------------
92  !
93  !++ Public procedure
94  !
95  public :: mkinit_setup
96  public :: mkinit
97 
98  !-----------------------------------------------------------------------------
99  !
100  !++ Public parameters & variables
101  !
102  integer, public :: mkinit_type = -1
103  integer, public, parameter :: i_ignore = 0
104 
105  integer, public, parameter :: i_planestate = 1
106  integer, public, parameter :: i_tracerbubble = 2
107  integer, public, parameter :: i_coldbubble = 3
108 
109  integer, public, parameter :: i_lambwave = 4
110  integer, public, parameter :: i_gravitywave = 5
111  integer, public, parameter :: i_khwave = 6
112  integer, public, parameter :: i_turbulence = 7
113  integer, public, parameter :: i_mountainwave = 8
114 
115  integer, public, parameter :: i_warmbubble = 9
116  integer, public, parameter :: i_supercell = 10
117  integer, public, parameter :: i_squallline = 11
118  integer, public, parameter :: i_wk1982 = 12
119  integer, public, parameter :: i_dycoms2_rf01 = 13
120  integer, public, parameter :: i_dycoms2_rf02 = 14
121  integer, public, parameter :: i_rico = 15
122 
123  integer, public, parameter :: i_interporation = 16
124 
125  integer, public, parameter :: i_landcouple = 17
126  integer, public, parameter :: i_oceancouple = 18
127  integer, public, parameter :: i_urbancouple = 19
128  integer, public, parameter :: i_triplecouple = 20
129  integer, public, parameter :: i_bubblecouple = 21
130 
131  integer, public, parameter :: i_seabreeze = 22
132  integer, public, parameter :: i_heatisland = 23
133 
134  integer, public, parameter :: i_dycoms2_rf02_dns = 24
135 
136  integer, public, parameter :: i_real = 25
137 
138  integer, public, parameter :: i_grayzone = 26
139  integer, public, parameter :: i_boxaero = 27
140  integer, public, parameter :: i_warmbubbleaero = 28
141 
142  integer, public, parameter :: i_cavityflow = 29
143 
144  !-----------------------------------------------------------------------------
145  !
146  !++ Private procedure
147  !
148  private :: bubble_setup
149  private :: sbmaero_setup
150  private :: aerosol_setup
151 
152  private :: mkinit_planestate
153  private :: mkinit_tracerbubble
154  private :: mkinit_coldbubble
155  private :: mkinit_lambwave
156  private :: mkinit_gravitywave
157  private :: mkinit_khwave
158  private :: mkinit_turbulence
159  private :: mkinit_cavityflow
160  private :: mkinit_mountainwave
161 
162  private :: mkinit_warmbubble
163  private :: mkinit_supercell
164  private :: mkinit_squallline
165  private :: mkinit_wk1982
166  private :: mkinit_dycoms2_rf01
167  private :: mkinit_dycoms2_rf02
168  private :: mkinit_rico
169 
170  private :: mkinit_interporation
171 
172  private :: mkinit_landcouple
173  private :: mkinit_oceancouple
174  private :: mkinit_urbancouple
175  private :: mkinit_seabreeze
176  private :: mkinit_heatisland
177 
178  private :: mkinit_dycoms2_rf02_dns
179 
180  private :: mkinit_real
181 
182  private :: mkinit_grayzone
183 
184  private :: mkinit_boxaero
185  private :: mkinit_warmbubbleaero
186 
187  !-----------------------------------------------------------------------------
188  !
189  !++ Private parameters & variables
190  !
191  real(RP), private, parameter :: thetastd = 300.0_rp ! [K]
192 
193  real(RP), private, allocatable :: pres (:,:,:) ! pressure [Pa]
194  real(RP), private, allocatable :: temp (:,:,:) ! temperature [K]
195  real(RP), private, allocatable :: pott (:,:,:) ! potential temperature [K]
196  real(RP), private, allocatable :: qsat (:,:,:) ! satulated water vapor [kg/kg]
197  real(RP), private, allocatable :: qv (:,:,:) ! water vapor [kg/kg]
198  real(RP), private, allocatable :: qc (:,:,:) ! cloud water [kg/kg]
199  real(RP), private, allocatable :: velx (:,:,:) ! velocity u [m/s]
200  real(RP), private, allocatable :: vely (:,:,:) ! velocity v [m/s]
201 
202  real(RP), private, allocatable :: pres_sfc(:,:,:) ! surface pressure [Pa]
203  real(RP), private, allocatable :: temp_sfc(:,:,:) ! surface temperature [K]
204  real(RP), private, allocatable :: pott_sfc(:,:,:) ! surface potential temperature [K]
205  real(RP), private, allocatable :: qsat_sfc(:,:,:) ! surface satulated water vapor [kg/kg]
206  real(RP), private, allocatable :: qv_sfc (:,:,:) ! surface water vapor [kg/kg]
207  real(RP), private, allocatable :: qc_sfc (:,:,:) ! surface cloud water [kg/kg]
208 
209  real(RP), private, allocatable :: rndm (:,:,:) ! random number (0-1)
210  real(RP), private, allocatable, target :: bubble (:,:,:) ! bubble factor (0-1)
211  real(RP), private, allocatable, target :: rect (:,:,:) ! rectangle factor (0-1)
212  real(RP), private, allocatable :: gan (:) ! gamma factor (0-1)
213 
214  logical, private :: flg_intrp = .false.
215 
216  !-----------------------------------------------------------------------------
217 contains
218  !-----------------------------------------------------------------------------
220  subroutine mkinit_setup
221  implicit none
222 
223  character(len=H_SHORT) :: MKINIT_initname = 'NONE'
224 
225  namelist / param_mkinit / &
226  mkinit_initname, &
227  flg_intrp
228 
229  integer :: ierr
230  !---------------------------------------------------------------------------
231 
232  if( io_l ) write(io_fid_log,*)
233  if( io_l ) write(io_fid_log,*) '++++++ Module[make init] / Categ[preprocess] / Origin[SCALE-RM]'
234 
235  !--- read namelist
236  rewind(io_fid_conf)
237  read(io_fid_conf,nml=param_mkinit,iostat=ierr)
238  if( ierr < 0 ) then !--- missing
239  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
240  elseif( ierr > 0 ) then !--- fatal error
241  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT. Check!'
242  call prc_mpistop
243  endif
244  if( io_lnml ) write(io_fid_log,nml=param_mkinit)
245 
246  allocate( pres(ka,ia,ja) )
247  allocate( temp(ka,ia,ja) )
248  allocate( pott(ka,ia,ja) )
249  allocate( qsat(ka,ia,ja) )
250  allocate( qv(ka,ia,ja) )
251  allocate( qc(ka,ia,ja) )
252  allocate( velx(ka,ia,ja) )
253  allocate( vely(ka,ia,ja) )
254 
255  allocate( pres_sfc(1,ia,ja) )
256  allocate( temp_sfc(1,ia,ja) )
257  allocate( pott_sfc(1,ia,ja) )
258  allocate( qsat_sfc(1,ia,ja) )
259  allocate( qv_sfc(1,ia,ja) )
260  allocate( qc_sfc(1,ia,ja) )
261 
262  allocate( rndm(ka,ia,ja) )
263  allocate( bubble(ka,ia,ja) )
264  allocate( rect(ka,ia,ja) )
265 
266  select case(trim(mkinit_initname))
267  case('NONE')
269  case('PLANESTATE')
271  case('TRACERBUBBLE')
273  case('COLDBUBBLE')
275  call bubble_setup
276  case('LAMBWAVE')
278  call bubble_setup
279  case('GRAVITYWAVE')
281  call bubble_setup
282  case('KHWAVE')
284  case('TURBULENCE')
286  case('MOUNTAINWAVE')
288  call bubble_setup
289  case('WARMBUBBLE')
291  call bubble_setup
292  case('SUPERCELL')
294  call bubble_setup
295  case('SQUALLLINE')
297  case('WK1982')
299  call bubble_setup
300  case('DYCOMS2_RF01')
302  case('DYCOMS2_RF02')
304  case('RICO')
306  case('INTERPORATION')
308  case('LANDCOUPLE')
310  case('OCEANCOUPLE')
312  case('URBANCOUPLE')
314  case('TRIPLECOUPLE')
316  case('BUBBLECOUPLE')
318  call bubble_setup
319  case('SEABREEZE')
321  case('HEATISLAND')
323  case('DYCOMS2_RF02_DNS')
325  case('REAL')
327  case('GRAYZONE')
329  case('BOXAERO')
331  case('WARMBUBBLEAERO')
333  call bubble_setup
334  case('CAVITYFLOW')
336  case default
337  write(*,*) ' xxx Unsupported TYPE:', trim(mkinit_initname)
338  call prc_mpistop
339  endselect
340 
341  return
342  end subroutine mkinit_setup
343 
344  !-----------------------------------------------------------------------------
346  subroutine mkinit
347  use scale_const, only: &
349  use scale_landuse, only: &
351  use mod_atmos_driver, only: &
353  use mod_atmos_vars, only: &
354  atmos_sw_restart => atmos_restart_output, &
356  use mod_ocean_driver, only: &
358  use mod_ocean_vars, only: &
359  ocean_sw_restart => ocean_restart_output, &
361  use mod_land_driver, only: &
363  use mod_land_vars, only: &
364  land_sw_restart => land_restart_output, &
366  use mod_urban_driver, only: &
368  use mod_urban_vars, only: &
369  urban_sw_restart => urban_restart_output, &
371  implicit none
372 
373  integer :: iq
374  !---------------------------------------------------------------------------
375 
376  if ( mkinit_type == i_ignore ) then
377  if( io_l ) write(io_fid_log,*)
378  if( io_l ) write(io_fid_log,*) '++++++ SKIP MAKING INITIAL DATA ++++++'
379  else
380  if( io_l ) write(io_fid_log,*)
381  if( io_l ) write(io_fid_log,*) '++++++ START MAKING INITIAL DATA ++++++'
382 
383  !--- Initialize variables
384 #ifndef DRY
385  do iq = 2, qa
386  qtrc(:,:,:,iq) = 0.0_rp
387  enddo
388 #endif
389 
390  pres(:,:,:) = const_undef8
391  temp(:,:,:) = const_undef8
392  pott(:,:,:) = const_undef8
393  qsat(:,:,:) = const_undef8
394  qv(:,:,:) = const_undef8
395  qc(:,:,:) = const_undef8
396  velx(:,:,:) = const_undef8
397  vely(:,:,:) = const_undef8
398 
399  rndm(:,:,:) = const_undef8
400 
401  pres_sfc(:,:,:) = const_undef8
402  temp_sfc(:,:,:) = const_undef8
403  pott_sfc(:,:,:) = const_undef8
404  qsat_sfc(:,:,:) = const_undef8
405  qv_sfc(:,:,:) = const_undef8
406  qc_sfc(:,:,:) = const_undef8
407 
408  if ( tracer_type == 'SUZUKI10' ) then
409  if( io_l ) write(io_fid_log,*) '*** Aerosols for SBM are included ***'
410  call sbmaero_setup
411  endif
412 
413  select case(mkinit_type)
414  case(i_planestate)
415  call mkinit_planestate
416  case(i_tracerbubble)
417  call mkinit_tracerbubble
418  case(i_coldbubble)
419  call mkinit_coldbubble
420  case(i_lambwave)
421  call mkinit_lambwave
422  case(i_gravitywave)
423  call mkinit_gravitywave
424  case(i_khwave)
425  call mkinit_khwave
426  case(i_turbulence)
427  call mkinit_turbulence
428  case(i_mountainwave)
429  call mkinit_mountainwave
430  case(i_warmbubble)
431  call mkinit_warmbubble
432  case(i_supercell)
433  call mkinit_supercell
434  case(i_squallline)
435  call mkinit_squallline
436  case(i_wk1982)
437  call mkinit_wk1982
438  case(i_dycoms2_rf01)
439  call mkinit_dycoms2_rf01
440  case(i_dycoms2_rf02)
441  call mkinit_dycoms2_rf02
442  case(i_rico)
443  call mkinit_rico
444  case(i_interporation)
445  call mkinit_interporation
446  case(i_oceancouple)
447  call mkinit_planestate
448  call mkinit_oceancouple
449  case(i_landcouple)
450  call mkinit_planestate
451  call mkinit_landcouple
452  case(i_urbancouple)
453  call mkinit_planestate
454  call mkinit_landcouple ! tentative
455  call mkinit_urbancouple
456  case(i_triplecouple)
457  call mkinit_planestate
458  call mkinit_oceancouple
459  call mkinit_landcouple
460  call mkinit_urbancouple
461  case(i_bubblecouple)
462  call mkinit_planestate
463  call mkinit_warmbubble
464  call mkinit_oceancouple
465  call mkinit_landcouple
466  call mkinit_urbancouple
467  case(i_seabreeze)
468  call mkinit_planestate
469  call mkinit_seabreeze
470  case(i_heatisland)
471  call mkinit_planestate
472  call mkinit_heatisland
473  case(i_dycoms2_rf02_dns)
474  call mkinit_dycoms2_rf02_dns
475  case(i_real)
476  call mkinit_real
477  case(i_grayzone)
478  call mkinit_grayzone
479  case(i_boxaero)
480  call mkinit_boxaero
481  case(i_warmbubbleaero)
482  call mkinit_warmbubbleaero
483  case(i_cavityflow)
484  call mkinit_cavityflow
485  case default
486  write(*,*) ' xxx Unsupported TYPE:', mkinit_type
487  call prc_mpistop
488  endselect
489 
490  if( io_l ) write(io_fid_log,*) '++++++ END MAKING INITIAL DATA ++++++'
491 
492  ! setup surface condition
493  call ocean_surface_set( countup = .false. )
494  call land_surface_set ( countup = .false. )
495  call urban_surface_set( countup = .false. )
496  call atmos_surface_get
497 
498  ! output boundary file
499  call landuse_write
500 
501  ! output restart file
502  if( atmos_sw_restart ) call atmos_vars_restart_write
503  if( ocean_sw_restart ) call ocean_vars_restart_write
504  if( land_sw_restart ) call land_vars_restart_write
505  if( urban_sw_restart ) call urban_vars_restart_write
506  endif
507 
508  return
509  end subroutine mkinit
510 
511  !-----------------------------------------------------------------------------
513  subroutine bubble_setup
514  use scale_const, only: &
516  implicit none
517 
518  ! Bubble
519  logical :: BBL_eachnode = .false. ! Arrange bubble at each node? [kg/kg]
520  real(RP) :: BBL_CZ = 2.e3_rp ! center location [m]: z
521  real(RP) :: BBL_CX = 2.e3_rp ! center location [m]: x
522  real(RP) :: BBL_CY = 2.e3_rp ! center location [m]: y
523  real(RP) :: BBL_RZ = 0.0_rp ! bubble radius [m]: z
524  real(RP) :: BBL_RX = 0.0_rp ! bubble radius [m]: x
525  real(RP) :: BBL_RY = 0.0_rp ! bubble radius [m]: y
526 
527  namelist / param_bubble / &
528  bbl_eachnode, &
529  bbl_cz, &
530  bbl_cx, &
531  bbl_cy, &
532  bbl_rz, &
533  bbl_rx, &
534  bbl_ry
535 
536  real(RP) :: CZ_offset
537  real(RP) :: CX_offset
538  real(RP) :: CY_offset
539  real(RP) :: dist
540 
541  integer :: ierr
542  integer :: k, i, j
543  !---------------------------------------------------------------------------
544 
545  if( io_l ) write(io_fid_log,*)
546  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit bubble] / Categ[preprocess] / Origin[SCALE-RM]'
547 
548  !--- read namelist
549  rewind(io_fid_conf)
550  read(io_fid_conf,nml=param_bubble,iostat=ierr)
551 
552  if( ierr < 0 ) then !--- missing
553  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
554  elseif( ierr > 0 ) then !--- fatal error
555  write(*,*) 'xxx Not appropriate names in namelist PARAM_BUBBLE. Check!'
556  call prc_mpistop
557  endif
558  if( io_lnml ) write(io_fid_log,nml=param_bubble)
559 
560  if ( abs(bbl_rz*bbl_rx*bbl_ry) <= 0.0_rp ) then
561  if( io_l ) write(io_fid_log,*) '*** no bubble'
562  bubble(:,:,:) = 0.0_rp
563  else
564 
565  bubble(:,:,:) = const_undef8
566 
567  if ( bbl_eachnode ) then
568  cz_offset = grid_cz(ks)
569  cx_offset = grid_cx(is)
570  cy_offset = grid_cy(js)
571  else
572  cz_offset = 0.0_rp
573  cx_offset = 0.0_rp
574  cy_offset = 0.0_rp
575  endif
576 
577  do j = 1, ja
578  do i = 1, ia
579  do k = ks, ke
580 
581  ! make tracer bubble
582  dist = ( (grid_cz(k)-cz_offset-bbl_cz)/bbl_rz )**2 &
583  + ( (grid_cx(i)-cx_offset-bbl_cx)/bbl_rx )**2 &
584  + ( (grid_cy(j)-cy_offset-bbl_cy)/bbl_ry )**2
585 
586  bubble(k,i,j) = cos( 0.5_rp*pi*sqrt( min(dist,1.0_rp) ) )**2
587 
588  enddo
589  enddo
590  enddo
591  endif
592 
593  return
594  end subroutine bubble_setup
595 
596  !-----------------------------------------------------------------------------
598  subroutine rect_setup
599  use scale_const, only: &
601  implicit none
602 
603  ! Bubble
604  logical :: RCT_eachnode = .false. ! Arrange rectangle at each node? [kg/kg]
605  real(RP) :: RCT_CZ = 2.e3_rp ! center location [m]: z
606  real(RP) :: RCT_CX = 2.e3_rp ! center location [m]: x
607  real(RP) :: RCT_CY = 2.e3_rp ! center location [m]: y
608  real(RP) :: RCT_RZ = 2.e3_rp ! rectangle z width [m]: z
609  real(RP) :: RCT_RX = 2.e3_rp ! rectangle x width [m]: x
610  real(RP) :: RCT_RY = 2.e3_rp ! rectangle y width [m]: y
611 
612  namelist / param_rect / &
613  rct_eachnode, &
614  rct_cz, &
615  rct_cx, &
616  rct_cy, &
617  rct_rz, &
618  rct_rx, &
619  rct_ry
620 
621  real(RP) :: CZ_offset
622  real(RP) :: CX_offset
623  real(RP) :: CY_offset
624  real(RP) :: dist
625 
626  integer :: ierr
627  integer :: k, i, j
628  !---------------------------------------------------------------------------
629 
630  if( io_l ) write(io_fid_log,*)
631  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit rectangle] / Categ[preprocess] / Origin[SCALE-RM]'
632 
633  !--- read namelist
634  rewind(io_fid_conf)
635  read(io_fid_conf,nml=param_rect,iostat=ierr)
636 
637  if( ierr < 0 ) then !--- missing
638  if( io_l ) write(io_fid_log,*) 'xxx Not found namelist. Check!'
639  call prc_mpistop
640  elseif( ierr > 0 ) then !--- fatal error
641  write(*,*) 'xxx Not appropriate names in namelist PARAM_RECT. Check!'
642  call prc_mpistop
643  endif
644  if( io_lnml ) write(io_fid_log,nml=param_rect)
645 
646  rect(:,:,:) = const_undef8
647 
648  if ( rct_eachnode ) then
649  cz_offset = grid_cz(ks)
650  cx_offset = grid_cx(is)
651  cy_offset = grid_cy(js)
652  else
653  cz_offset = 0.0_rp
654  cx_offset = 0.0_rp
655  cy_offset = 0.0_rp
656  endif
657 
658  do j = 1, ja
659  do i = 1, ia
660  do k = ks, ke
661 
662  ! make tracer rectangle
663  dist = 2.0_rp * max( &
664  abs(grid_cz(k) - cz_offset - rct_cz)/rct_rz, &
665  abs(grid_cx(i) - cx_offset - rct_cx)/rct_rx, &
666  abs(grid_cy(j) - cy_offset - rct_cy)/rct_ry &
667  & )
668  if ( dist <= 1.0_rp ) then
669  rect(k,i,j) = 1.0_rp
670  else
671  rect(k,i,j) = 0.0_rp
672  end if
673  enddo
674  enddo
675  enddo
676 
677  return
678  end subroutine rect_setup
679 
680  !-----------------------------------------------------------------------------
682  subroutine aerosol_setup
683 
684  use scale_tracer
685  use scale_const, only: &
686  pi => const_pi, &
687  const_cvdry, &
688  const_rdry, &
689  const_rvap, &
691  use scale_atmos_thermodyn, only: &
692  aq_cv
693  implicit none
694 
695  integer, parameter :: ia_m0 = 1 !1. number conc [#/m3]
696  integer, parameter :: ia_m2 = 2 !2. 2nd mom conc [m2/m3]
697  integer, parameter :: ia_m3 = 3 !3. 3rd mom conc [m3/m3]
698  integer, parameter :: ia_ms = 4 !4. mass conc [ug/m3]
699  integer, parameter :: ia_kp = 5
700  integer, parameter :: ik_out = 1
701 
702  real(RP) :: m0_init = 0.0_rp ! initial total num. conc. of modes (Atk,Acm,Cor) [#/m3]
703  real(RP) :: dg_init = 80.e-9_rp ! initial number equivalen diameters of modes [m]
704  real(RP) :: sg_init = 1.6_rp ! initial standard deviation [-]
705 
706  real(RP), parameter :: d_min_def = 1.e-9_rp ! default lower bound of 1st size bin
707  real(RP), parameter :: d_max_def = 1.e-5_rp ! upper bound of last size bin
708  integer, parameter :: n_kap_def = 1 ! number of kappa bins
709  real(RP), parameter :: k_min_def = 0.e0_rp ! lower bound of 1st kappa bin
710  real(RP), parameter :: k_max_def = 1.e0_rp ! upper bound of last kappa bin
711  real(RP) :: c_kappa = 0.3_rp ! hygroscopicity of condensable mass [-]
712  real(RP), parameter :: cleannumber = 1.e-3_rp
713 
714  real(RP) :: d_min_inp(3) = d_min_def
715  real(RP) :: d_max_inp(3) = d_max_def
716  real(RP) :: k_min_inp(3) = k_min_def
717  real(RP) :: k_max_inp(3) = k_max_def
718  integer :: n_kap_inp(3) = n_kap_def
719 
720  namelist / param_aero / &
721  m0_init, &
722  dg_init, &
723  sg_init, &
724  d_min_inp, &
725  d_max_inp, &
726  k_min_inp, &
727  k_max_inp, &
728  n_kap_inp
729 
730  ! local variables
731  real(RP), parameter :: rhod_ae = 1.83_rp ! particle density [g/cm3] sulfate assumed
732  real(RP), parameter :: conv_vl_ms = rhod_ae/1.e-12_rp ! M3(volume)[m3/m3] to mass[m3/m3]
733 
734  integer :: ia0, ik, is0, ic, k, i, j, it, iq
735  integer :: ierr, n_trans
736  real(RP) :: m0t, dgt, sgt, m2t, m3t, mst
737  real(RP),allocatable :: aerosol_procs(:,:,:,:) !(n_atr,n_siz_max,n_kap_max,n_ctg)
738  real(RP),allocatable :: aerosol_activ(:,:,:,:) !(n_atr,n_siz_max,n_kap_max,n_ctg)
739  real(RP),allocatable :: emis_procs(:,:,:,:) !(n_atr,n_siz_max,n_kap_max,n_ctg)
740  !--- gas
741  real(RP) :: conc_gas(gas_ctg) !concentration [ug/m3]
742  integer :: n_siz_max, n_kap_max, n_ctg
743 ! integer, allocatable :: it_procs2trans(:,:,:,:) !procs to trans conversion
744 ! integer, allocatable :: ia_trans2procs(:) !trans to procs conversion
745 ! integer, allocatable :: is_trans2procs(:) !trans to procs conversion
746 ! integer, allocatable :: ik_trans2procs(:) !trans to procs conversion
747 ! integer, allocatable :: ic_trans2procs(:)
748  !--- bin settings (lower bound, center, upper bound)
749  real(RP),allocatable :: d_lw(:,:), d_ct(:,:), d_up(:,:) !diameter [m]
750  real(RP),allocatable :: k_lw(:,:), k_ct(:,:), k_up(:,:) !kappa [-]
751  real(RP) :: dlogd, dk !delta log(D), delta K
752  real(RP),allocatable :: d_min(:) !lower bound of 1st size bin (n_ctg)
753  real(RP),allocatable :: d_max(:) !upper bound of last size bin (n_ctg)
754  integer, allocatable :: n_kap(:) !number of kappa bins (n_ctg)
755  real(RP),allocatable :: k_min(:) !lower bound of 1st kappa bin (n_ctg)
756  real(RP),allocatable :: k_max(:)
757 
758  real(RP) :: pott, qdry, pres
759  real(RP) :: temp, cpa, cva, qsat_tmp, ssliq, Rmoist
760  real(RP) :: pi6
761 
762  if( io_l ) write(io_fid_log,*)
763  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit aerosol] / Categ[preprocess] / Origin[SCALE-RM]'
764 
765  pi6 = pi / 6._rp
766  n_ctg = ae_ctg
767 
768  allocate( d_min(n_ctg) )
769  allocate( d_max(n_ctg) )
770  allocate( n_kap(n_ctg) )
771  allocate( k_min(n_ctg) )
772  allocate( k_max(n_ctg) )
773 
774  d_min(1:n_ctg) = d_min_def ! lower bound of 1st size bin
775  d_max(1:n_ctg) = d_max_def ! upper bound of last size bin
776  n_kap(1:n_ctg) = n_kap_def ! number of kappa bins
777  k_min(1:n_ctg) = k_min_def ! lower bound of 1st kappa bin
778  k_max(1:n_ctg) = k_max_def
779  d_min_inp(1:n_ctg) = d_min(1:n_ctg) !def value used if not defined in namelist
780  d_max_inp(1:n_ctg) = d_max(1:n_ctg) !def value used if not defined in namelist
781  n_kap_inp(1:n_ctg) = n_kap(1:n_ctg) !def value used if not defined in namelist
782  k_min_inp(1:n_ctg) = k_min(1:n_ctg) !def value used if not defined in namelist
783  k_max_inp(1:n_ctg) = k_max(1:n_ctg) !def value used if not defined in namelist
784 
785  !--- read namelist
786  rewind(io_fid_conf)
787  read(io_fid_conf,nml=param_aero,iostat=ierr)
788 
789  if( ierr < 0 ) then !--- missing
790  if( io_l ) write(io_fid_log,*) 'xxx Not found namelist. Default used!'
791 ! call PRC_MPIstop
792  elseif( ierr > 0 ) then !--- fatal error
793  write(*,*) 'xxx Not appropriate names in namelist PARAM_AERO. Check!'
794  call prc_mpistop
795  else
796  d_min(1:n_ctg) = d_min_inp(1:n_ctg) ! lower bound of 1st size bin
797  d_max(1:n_ctg) = d_max_inp(1:n_ctg) ! upper bound of last size bin
798  n_kap(1:n_ctg) = n_kap_inp(1:n_ctg) ! number of kappa bins
799  k_min(1:n_ctg) = k_min_inp(1:n_ctg) ! lower bound of 1st kappa bin
800  k_max(1:n_ctg) = k_max_inp(1:n_ctg) ! upper bound of last kappa bin
801  endif
802  if( io_lnml ) write(io_fid_log,nml=param_aero)
803 
804  !--- diagnose parameters (n_trans, n_siz_max, n_kap_max)
805  n_trans = 0
806  n_siz_max = 0
807  n_kap_max = 0
808  do ic = 1, n_ctg
809  n_trans = n_trans + nsiz(ic) * nkap(ic) * n_atr
810  n_siz_max = max(n_siz_max, nsiz(ic))
811  n_kap_max = max(n_kap_max, nkap(ic))
812  enddo
813 
814  allocate( aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg) )
815  allocate( aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg) )
816  allocate( emis_procs(n_atr,n_siz_max,n_kap_max,n_ctg) )
817  aerosol_procs(:,:,:,:) = 0._rp
818  aerosol_activ(:,:,:,:) = 0._rp
819  emis_procs(:,:,:,:) = 0._rp
820 
821 ! allocate( it_procs2trans(N_ATR,n_siz_max,n_kap_max,n_ctg) )
822 ! allocate( ia_trans2procs(n_trans) )
823 ! allocate( is_trans2procs(n_trans) )
824 ! allocate( ik_trans2procs(n_trans) )
825 ! allocate( ic_trans2procs(n_trans) )
826 
827  !bin setting
828  allocate(d_lw(n_siz_max,n_ctg))
829  allocate(d_ct(n_siz_max,n_ctg))
830  allocate(d_up(n_siz_max,n_ctg))
831  allocate(k_lw(n_kap_max,n_ctg))
832  allocate(k_ct(n_kap_max,n_ctg))
833  allocate(k_up(n_kap_max,n_ctg))
834  d_lw(:,:) = 0._rp
835  d_ct(:,:) = 0._rp
836  d_up(:,:) = 0._rp
837  k_lw(:,:) = 0._rp
838  k_ct(:,:) = 0._rp
839  k_up(:,:) = 0._rp
840 
841  do ic = 1, n_ctg
842 
843  dlogd = (log(d_max(ic)) - log(d_min(ic)))/real(NSIZ(ic),kind=rp)
844 
845  do is0 = 1, nsiz(ic) !size bin
846  d_lw(is0,ic) = exp(log(d_min(ic))+dlogd* real(is0-1,kind=RP) )
847  d_ct(is0,ic) = exp(log(d_min(ic))+dlogd*(real(is0 ,kind=rp)-0.5_rp))
848  d_up(is0,ic) = exp(log(d_min(ic))+dlogd* real(is0 ,kind=RP) )
849  enddo !is (1:n_siz(ic))
850  dk = (k_max(ic) - k_min(ic))/real(n_kap(ic),kind=rp)
851  do ik = 1, n_kap(ic) !size bin
852  k_lw(ik,ic) = k_min(ic) + dk * real(ik-1,kind=rp)
853  k_ct(ik,ic) = k_min(ic) + dk *(real(ik ,kind=rp)-0.5_rp)
854  k_up(ik,ic) = k_min(ic) + dk * real(ik ,kind=rp)
855  enddo !ik (1:n_kap(ic))
856 
857  enddo !ic (1:n_ctg)
858 ! ik = 1 !only one kappa bin
859 
860  m0t = m0_init !total M0 [#/m3]
861  dgt = dg_init ![m]
862  sgt = sg_init ![-]
863 
864  if ( m0t <= cleannumber ) then
865  m0t = cleannumber
866  dgt = 0.1e-6_rp
867  sgt = 1.3_rp
868  if ( io_l ) then
869  write(io_fid_log,*) '*** WARNING! Initial aerosol number is set as ', cleannumber, '[#/m3]'
870  endif
871  endif
872 
873  m2t = m0t*dgt**(2.d0) *dexp(2.0d0 *(dlog(real(sgt,kind=dp))**2.d0)) !total M2 [m2/m3]
874  m3t = m0t*dgt**(3.d0) *dexp(4.5d0 *(dlog(real(sgt,kind=dp))**2.d0)) !total M3 [m3/m3]
875  mst = m3t*pi6*conv_vl_ms !total Ms [ug/m3]
876 
877  do ic = 1, n_ctg
878  !aerosol_procs initial condition
879  do ik = 1, n_kap(ic) !kappa bin
880  do is0 = 1, nsiz(ic)
881  if (dgt >= d_lw(is0,ic) .and. dgt < d_up(is0,ic)) then
882  aerosol_procs(ia_m0,is0,ik,ic) = aerosol_procs(ia_m0,is0,ik,ic) + m0t ![#/m3]
883  aerosol_procs(ia_m2,is0,ik,ic) = aerosol_procs(ia_m2,is0,ik,ic) + m2t ![m2/m3]
884  aerosol_procs(ia_m3,is0,ik,ic) = aerosol_procs(ia_m3,is0,ik,ic) + m3t ![m3/m3]
885  aerosol_procs(ia_ms,is0,ik,ic) = aerosol_procs(ia_ms,is0,ik,ic) + mst*1.e-9_rp ![kg/m3]
886  elseif (dgt < d_lw(1,ic)) then
887  aerosol_procs(ia_m0,1 ,ik,ic) = aerosol_procs(ia_m0,1 ,ik,ic) + m0t ![#/m3]
888  aerosol_procs(ia_m2,1 ,ik,ic) = aerosol_procs(ia_m2,1 ,ik,ic) + m2t ![m2/m3]
889  aerosol_procs(ia_m3,1 ,ik,ic) = aerosol_procs(ia_m3,1 ,ik,ic) + m3t ![m3/m3]
890  aerosol_procs(ia_ms,1 ,ik,ic) = aerosol_procs(ia_ms,1 ,ik,ic) + mst*1.e-9_rp ![kg/m3]
891  elseif (dgt >= d_up(nsiz(ic),ic)) then
892  aerosol_procs(ia_m0,nsiz(ic),ik,ic) = aerosol_procs(ia_m0,nsiz(ic),ik,ic) + m0t ![#/m3]
893  aerosol_procs(ia_m2,nsiz(ic),ik,ic) = aerosol_procs(ia_m2,nsiz(ic),ik,ic) + m2t ![m2/m3]
894  aerosol_procs(ia_m3,nsiz(ic),ik,ic) = aerosol_procs(ia_m3,nsiz(ic),ik,ic) + m3t ![m3/m3]
895  aerosol_procs(ia_ms,nsiz(ic),ik,ic) = aerosol_procs(ia_ms,nsiz(ic),ik,ic) + mst*1.e-9_rp ![kg/m3]
896  endif
897  enddo
898  enddo
899  enddo
900 
901  ccn(:,:,:) = 0.0_rp
902  do j = js, je
903  do i = is, ie
904 
905  do k = ks, ke
906  pott = rhot(k,i,j) / dens(k,i,j)
907  qdry = 1.0_rp
908  do iq = qqs, qqe
909  qdry = qdry - qtrc(k,i,j,iq)
910  enddo
911  cva = qdry * const_cvdry
912  do iq = qqs, qqe
913  cva = cva + qtrc(k,i,j,iq) * aq_cv(iq)
914  enddo
915  rmoist = const_rdry * qdry + const_rvap * qtrc(k,i,j,i_qv)
916  cpa = cva + rmoist
917  pres = const_pre00 &
918  * ( dens(k,i,j)*rmoist*pott/const_pre00 )**( cpa/(cpa-rmoist) )
919  temp = pres / ( dens(k,i,j) * rmoist )
920  !--- calculate super saturation of water
921  call saturation_pres2qsat_liq( qsat_tmp,temp,pres )
922  ssliq = qtrc(k,i,j,i_qv)/qsat_tmp - 1.0_rp
924  (c_kappa, ssliq, temp, ia_m0, ia_m2, ia_m3, &
925  n_atr,n_siz_max,n_kap_max,n_ctg,nsiz,n_kap, &
926  d_ct,aerosol_procs, aerosol_activ)
927 
928  do is0 = 1, nsiz(ic_mix)
929  ccn(k,i,j) = ccn(k,i,j) + aerosol_activ(ia_m0,is0,ik_out,ic_mix)
930  enddo
931 
932  enddo
933 
934  enddo
935  enddo
936 
937  conc_gas(:) = 0.0_rp
938  do k = ks, ke
939  do i = is, ie
940  do j = js, je
941  it = 0
942  do ic = 1, n_ctg !category
943  do ik = 1, nkap(ic) !kappa bin
944  do is0 = 1, nsiz(ic) !size bin
945  do ia0 = 1, n_atr !attributes
946  it = it + 1
947  qtrc(k,i,j,qaes-1+it) = aerosol_procs(ia0,is0,ik,ic)/dens(k,i,j) !#,m2,m3,kg/m3 -> #,m2,m3,kg/kg
948  enddo !ia0 (1:N_ATR )
949  enddo !is (1:n_siz(ic) )
950  enddo !ik (1:n_kap(ic) )
951  enddo !ic (1:n_ctg )
952  qtrc(k,i,j,qaee-gas_ctg+1:qaee-gas_ctg+gas_ctg) = conc_gas(1:gas_ctg)/dens(k,i,j)*1.e-9_rp !mixing ratio [kg/kg]
953  enddo
954  enddo
955  enddo
956 
957  return
958 
959  end subroutine aerosol_setup
960  !----------------------------------------------------------------------------------------
961  ! subroutine 4. aerosol_activation
962  ! Abdul-Razzak et al., JGR, 1998 [AR98]
963  ! Abdul-Razzak and Ghan, JGR, 2000 [AR00]
964  !----------------------------------------------------------------------------------------
965  subroutine aerosol_activation_kajino13(c_kappa, super, temp_k, &
966  ia_m0, ia_m2, ia_m3, &
967  n_atr,n_siz_max,n_kap_max,n_ctg,n_siz,n_kap, &
968  d_ct, aerosol_procs, aerosol_activ)
969  use scale_const, only: &
970  const_mvap, & ! mean molecular weight for water vapor [g/mol]
971  const_dwatr, & ! water density [kg/m3]
972  const_r, & ! universal gas constant [J/mol/K]
973  const_tem00
974  implicit none
975  !i/o variables
976  real(RP),intent(in) :: super ! supersaturation [-]
977  real(RP),intent(in) :: c_kappa ! hygroscopicity of condensable mass [-]
978  real(RP),intent(in) :: temp_k ! temperature
979  integer, intent(in) :: ia_m0, ia_m2, ia_m3
980  integer, intent(in) :: n_atr
981  integer, intent(in) :: n_siz_max
982  integer, intent(in) :: n_kap_max
983  integer, intent(in) :: n_ctg
984  real(RP),intent(in) :: d_ct(n_siz_max,n_ctg)
985  real(RP) :: aerosol_procs(n_atr,n_siz_max,n_kap_max,n_ctg)
986  real(RP) :: aerosol_activ(n_atr,n_siz_max,n_kap_max,n_ctg)
987  integer, intent(in) :: n_siz(n_ctg), n_kap(n_ctg)
988  !local variables
989  real(RP),parameter :: two3 = 2._rp/3._rp
990  real(RP),parameter :: rt2 = sqrt(2._rp)
991  real(RP),parameter :: twort2 = rt2 ! 2/sqrt(2) = sqrt(2)
992  real(RP),parameter :: thrrt2 = 3._rp/rt2 ! 3/sqrt(2)
993  real(RP) :: smax_inv ! inverse
994  real(RP) :: am,scrit_am,aa,tc,st,bb,ac
995  real(RP) :: m0t,m2t,m3t,dgt,sgt,dm2
996  real(RP) :: d_crit ! critical diameter
997  real(RP) :: tmp1, tmp2, tmp3 !
998  real(RP) :: ccn_frc,cca_frc,ccv_frc ! activated number,area,volume
999  integer :: is0, ik, ic
1000 
1001  aerosol_activ(:,:,:,:) = 0._rp
1002 
1003  if (super.le.0._rp) return
1004 
1005  smax_inv = 1._rp / super
1006 
1007  !--- surface tension of water
1008  tc = temp_k - const_tem00
1009  if (tc >= 0._rp ) then
1010  st = 75.94_rp-0.1365_rp*tc-0.3827e-3_rp*tc**2._rp !Gittens, JCIS, 69, Table 4.
1011  else !t[deg C].lt.0.
1012  st = 75.93_rp +0.115_rp*tc & !Eq.5-12, pp.130, PK97
1013  + 6.818e-2_rp*tc**2._rp+6.511e-3_rp*tc**3._rp &
1014  + 2.933e-4_rp*tc**4._rp+6.283e-6_rp*tc**5._rp &
1015  + 5.285e-8_rp*tc**6._rp
1016  endif
1017  st = st * 1.e-3_rp ![J/m2]
1018 
1019  !-- Kelvin effect
1020  ! [J m-2] [kg mol-1] [m3 kg-1] [mol K J-1] [K-1]
1021  aa = 2._rp * st * const_mvap * 1.e-3_rp / (const_dwatr * const_r * temp_k ) ![m] Eq.5 in AR98
1022 
1023  do ic = 1, n_ctg
1024  do ik = 1, n_kap(ic)
1025  do is0 = 1, n_siz(ic)
1026  m0t = aerosol_procs(ia_m0,is0,ik,ic)
1027  m2t = aerosol_procs(ia_m2,is0,ik,ic)
1028  m3t = aerosol_procs(ia_m3,is0,ik,ic)
1029  call diag_ds(m0t,m2t,m3t,dgt,sgt,dm2)
1030  if (dgt <= 0._rp) dgt = d_ct(is0,ic)
1031  if (sgt <= 0._rp) sgt = 1.3_rp
1032  am = dgt * 0.5_rp !geometric dry mean radius [m]
1033  bb = c_kappa
1034  if (bb > 0._rp .and. am > 0._rp ) then
1035  scrit_am = 2._rp/sqrt(bb)*(aa/(3._rp*am))**1.5_rp !AR00 Eq.9
1036  else
1037  scrit_am = 0._rp
1038  endif
1039  ac = am * (scrit_am * smax_inv)**two3 !AR00 Eq.12
1040  d_crit = ac * 2._rp
1041  tmp1 = log(d_crit) - log(dgt)
1042  tmp2 = 1._rp/(rt2*log(sgt))
1043  tmp3 = log(sgt)
1044  ccn_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2))
1045  cca_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2-twort2*tmp3))
1046  ccv_frc= 0.5_rp*(1._rp-erf(tmp1*tmp2-thrrt2*tmp3))
1047  aerosol_activ(ia_m0,is0,ik,ic) = ccn_frc * aerosol_procs(ia_m0,is0,ik,ic)
1048  aerosol_activ(ia_m2,is0,ik,ic) = cca_frc * aerosol_procs(ia_m2,is0,ik,ic)
1049  aerosol_activ(ia_m3,is0,ik,ic) = ccv_frc * aerosol_procs(ia_m3,is0,ik,ic)
1050  enddo !is(1:n_siz(ic))
1051  enddo !ik(1:n_kap(ic))
1052  enddo !ic(1:n_ctg)
1053 
1054  return
1055  end subroutine aerosol_activation_kajino13
1056  !----------------------------------------------------------------------------------------
1057  subroutine diag_ds(m0,m2,m3, & !i
1058  dg,sg,dm2) !o
1059  implicit none
1060  real(RP) :: m0,m2,m3,dg,sg,m3_bar,m2_bar
1061  real(RP) :: m2_new,m2_old,dm2
1062  real(RP), parameter :: sgmax=2.5_rp
1063  real(RP), parameter :: rk1=2._rp
1064  real(RP), parameter :: rk2=3._rp
1065  real(RP), parameter :: ratio =rk1/rk2
1066  real(RP), parameter :: rk1_hat=1._rp/(ratio*(rk2-rk1))
1067  real(RP), parameter :: rk2_hat=ratio/(rk1-rk2)
1068 
1069  dm2=0._rp
1070 
1071  if (m0 <= 0._rp .or. m2 <= 0._rp .or. m3 <= 0._rp) then
1072  m0=0._rp
1073  m2=0._rp
1074  m3=0._rp
1075  dg=-1._rp
1076  sg=-1._rp
1077  return
1078  endif
1079 
1080  m2_old = m2
1081  m3_bar = m3/m0
1082  m2_bar = m2/m0
1083  dg = m2_bar**rk1_hat*m3_bar**rk2_hat
1084 
1085  if (m2_bar/m3_bar**ratio < 1._rp) then !stdev > 1.
1086  sg = exp(sqrt(2._rp/(rk1*(rk1-rk2)) &
1087  * log(m2_bar/m3_bar**ratio) ))
1088  endif
1089 
1090  if (sg > sgmax) then
1091  print *,'sg=',sg
1092  sg = sgmax
1093  call diag_d2(m0,m3,sg, & !i
1094  m2,dg ) !o
1095  ! print *,'warning! modified as sg exceeded sgmax (diag_ds)'
1096  endif
1097 
1098  if (m2_bar/m3_bar**ratio >= 1._rp) then !stdev < 1.
1099  sg = 1._rp
1100  call diag_d2(m0,m3,sg, & !i
1101  m2,dg ) !o
1102  ! print *,'warning! modified as sg < 1. (diag_ds)'
1103  endif
1104 
1105  m2_new = m2
1106  dm2 = m2_old - m2_new !m2_pres - m2_diag
1107 
1108  return
1109  end subroutine diag_ds
1110  !----------------------------------------------------------------------------------------
1111  ! subroutine 6. diag_d2
1112  !----------------------------------------------------------------------------------------
1113  subroutine diag_d2(m0,m3,sg, & !i
1114  m2,dg ) !o
1115  implicit none
1116  real(RP) :: dg,sg,m0,m2,m3,aaa
1117  real(RP), parameter :: one3=1._rp/3._rp
1118  aaa = m0 * exp( 4.5_rp * (log(sg)**2._rp) )
1119  dg =(m3/aaa)**one3
1120  m2 = m0 * dg ** 2._rp * exp( 2.0_rp * (log(sg)**2._rp) )
1121 
1122  return
1123  end subroutine diag_d2
1124  !--------------------------------------------------------------
1126  subroutine sbmaero_setup
1127  use scale_const, only: &
1128  pi => const_pi
1129  use scale_tracer_suzuki10, only: &
1130  nccn, nbin
1131  implicit none
1132 
1133 #ifndef DRY
1134  real(RP) :: xasta, xaend, dxaer
1135  real(RP), allocatable :: xabnd( : ), xactr( : )
1136 
1137  real(RP) :: F0_AERO = 1.e+7_rp
1138  real(RP) :: R0_AERO = 1.e-7_rp
1139  real(RP) :: R_MAX = 1.e-06_rp
1140  real(RP) :: R_MIN = 1.e-08_rp
1141  real(RP) :: A_ALPHA = 3.0_rp
1142  real(RP) :: RHO_AERO = 2.25e+03_rp
1143 
1144  namelist / param_sbmaero / &
1145  f0_aero, &
1146  r0_aero, &
1147  r_max, &
1148  r_min, &
1149  a_alpha, &
1150  rho_aero
1151 
1152  integer :: ierr
1153  integer :: iq, i, j, k
1154  !---------------------------------------------------------------------------
1155 
1156  if( io_l ) write(io_fid_log,*)
1157  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit aerobin] / Categ[preprocess] / Origin[SCALE-RM]'
1158 
1159  !--- read namelist
1160  rewind(io_fid_conf)
1161  read(io_fid_conf,nml=param_sbmaero,iostat=ierr)
1162 
1163  if( ierr < 0 ) then !--- missing
1164  if( io_l ) write(io_fid_log,*) 'xxx Not found namelist. default value used'
1165  elseif( ierr > 0 ) then !--- fatal error
1166  write(*,*) 'xxx Not appropriate names in namelist SBMAERO. Check!'
1167  call prc_mpistop
1168  endif
1169  if( io_lnml ) write(io_fid_log,nml=param_sbmaero)
1170 
1171  if( nccn /= 0 ) then
1172  allocate( gan( nccn ) )
1173  allocate( xactr(nccn) )
1174  allocate( xabnd(nccn+1) )
1175 
1176  xasta = log( rho_aero*4.0_rp/3.0_rp*pi * ( r_min )**3 )
1177  xaend = log( rho_aero*4.0_rp/3.0_rp*pi * ( r_max )**3 )
1178  dxaer = ( xaend-xasta )/nccn
1179  do iq = 1, nccn+1
1180  xabnd( iq ) = xasta + dxaer*( iq-1 )
1181  enddo
1182  do iq = 1, nccn
1183  xactr( iq ) = ( xabnd( iq )+xabnd( iq+1 ) )*0.5_rp
1184  enddo
1185  do iq = 1, nccn
1186  gan( iq ) = faero( f0_aero,r0_aero,xactr( iq ), a_alpha, rho_aero )*exp( xactr(iq) )
1187  enddo
1188  endif
1189 
1190  !--- Hydrometeor is zero at initial time for Bin method
1191  do iq = 2, qqa
1192  do j = js, je
1193  do i = is, ie
1194  do k = ks, ke
1195  qtrc(k,i,j,iq) = 0.0_rp
1196  enddo
1197  enddo
1198  enddo
1199  enddo
1200 
1201  !-- Aerosol distribution
1202  if( nccn /= 0 ) then
1203  do iq = qqa+1, qa
1204  do j = js, je
1205  do i = is, ie
1206  do k = ks, ke
1207  qtrc(k,i,j,iq) = gan(iq-qqa) !/ DENS(k,i,j)
1208  enddo
1209  enddo
1210  enddo
1211  enddo
1212 
1213  deallocate( xactr )
1214  deallocate( xabnd )
1215  endif
1216 
1217 #endif
1218  return
1219  end subroutine sbmaero_setup
1220 
1221  !-----------------------------------------------------------------------------
1222  function faero( f0,r0,x,alpha,rhoa )
1223  use scale_const, only: &
1224  pi => const_pi
1225  implicit none
1226 
1227  real(RP), intent(in) :: x, f0, r0, alpha, rhoa
1228  real(RP) :: faero
1229  real(RP) :: rad
1230  !---------------------------------------------------------------------------
1231 
1232  rad = ( exp(x) * 3.0_rp / 4.0_rp / pi / rhoa )**(1.0_rp/3.0_rp)
1233 
1234  faero = f0 * (rad/r0)**(-alpha)
1235 
1236  return
1237  end function faero
1238 
1239  !-----------------------------------------------------------------------------
1241  subroutine flux_setup
1243  sflx_rain => atmos_phy_mp_sflx_rain, &
1244  sflx_snow => atmos_phy_mp_sflx_snow
1245  use mod_atmos_phy_rd_vars, only: &
1246  sflx_lw_up => atmos_phy_rd_sflx_lw_up, &
1247  sflx_lw_dn => atmos_phy_rd_sflx_lw_dn, &
1248  sflx_sw_up => atmos_phy_rd_sflx_sw_up, &
1249  sflx_sw_dn => atmos_phy_rd_sflx_sw_dn, &
1250  toaflx_lw_up => atmos_phy_rd_toaflx_lw_up, &
1251  toaflx_lw_dn => atmos_phy_rd_toaflx_lw_dn, &
1252  toaflx_sw_up => atmos_phy_rd_toaflx_sw_up, &
1253  toaflx_sw_dn => atmos_phy_rd_toaflx_sw_dn, &
1254  sflx_rad_dn => atmos_phy_rd_sflx_downall
1255  implicit none
1256  ! Flux from Atmosphere
1257  real(RP) :: FLX_rain = 0.0_rp ! surface rain flux [kg/m2/s]
1258  real(RP) :: FLX_snow = 0.0_rp ! surface snow flux [kg/m2/s]
1259  real(RP) :: FLX_LW_dn = 0.0_rp ! surface downwad long-wave radiation flux [J/m2/s]
1260  real(RP) :: FLX_SW_dn = 0.0_rp ! surface downwad short-wave radiation flux [J/m2/s]
1261 
1262  namelist / param_mkinit_flux / &
1263  flx_rain, &
1264  flx_snow, &
1265  flx_lw_dn, &
1266  flx_sw_dn
1267 
1268  integer :: i, j
1269  integer :: ierr
1270  !---------------------------------------------------------------------------
1271 
1272  !--- read namelist
1273  rewind(io_fid_conf)
1274  read(io_fid_conf,nml=param_mkinit_flux,iostat=ierr)
1275  if( ierr < 0 ) then !--- missing
1276  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1277  elseif( ierr > 0 ) then !--- fatal error
1278  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_FLUX. Check!'
1279  call prc_mpistop
1280  endif
1281  if( io_lnml ) write(io_fid_log,nml=param_mkinit_flux)
1282 
1283  do j = js, je
1284  do i = is, ie
1285  sflx_rain(i,j) = flx_rain
1286  sflx_snow(i,j) = flx_snow
1287 
1288  sflx_lw_up(i,j) = 0.0_rp
1289  sflx_lw_dn(i,j) = flx_lw_dn
1290  sflx_sw_up(i,j) = 0.0_rp
1291  sflx_sw_dn(i,j) = flx_sw_dn
1292 
1293  toaflx_lw_up(i,j) = 0.0_rp
1294  toaflx_lw_dn(i,j) = 0.0_rp
1295  toaflx_sw_up(i,j) = 0.0_rp
1296  toaflx_sw_dn(i,j) = 0.0_rp
1297 
1298  sflx_rad_dn(i,j,1,1) = flx_sw_dn
1299  sflx_rad_dn(i,j,1,2) = 0.0_rp
1300  sflx_rad_dn(i,j,2,1) = flx_lw_dn
1301  sflx_rad_dn(i,j,2,2) = 0.0_rp
1302  enddo
1303  enddo
1304 
1305  return
1306  end subroutine flux_setup
1307 
1308  !-----------------------------------------------------------------------------
1310  subroutine land_setup
1311  use mod_land_vars, only: &
1312  land_temp, &
1313  land_water, &
1314  land_sfc_temp, &
1316  implicit none
1317  ! Land state
1318  real(RP) :: LND_TEMP ! soil temperature [K]
1319  real(RP) :: LND_WATER = 0.15_rp ! soil moisture [m3/m3]
1320  real(RP) :: SFC_TEMP ! land skin temperature [K]
1321  real(RP) :: SFC_albedo_LW = 0.01_rp ! land surface albedo for LW [0-1]
1322  real(RP) :: SFC_albedo_SW = 0.20_rp ! land surface albedo for SW [0-1]
1323 
1324  integer :: i, j
1325  integer :: ierr
1326 
1327  namelist /param_mkinit_land/ &
1328  lnd_temp, &
1329  lnd_water, &
1330  sfc_temp, &
1331  sfc_albedo_lw, &
1332  sfc_albedo_sw
1333 
1334  lnd_temp = thetastd
1335  sfc_temp = thetastd
1336 
1337  !--- read namelist
1338  rewind(io_fid_conf)
1339  read(io_fid_conf,nml=param_mkinit_land,iostat=ierr)
1340  if( ierr < 0 ) then !--- missing
1341  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1342  elseif( ierr > 0 ) then !--- fatal error
1343  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_LAND. Check!'
1344  call prc_mpistop
1345  endif
1346  if( io_lnml ) write(io_fid_log,nml=param_mkinit_land)
1347 
1348  do j = js, je
1349  do i = is, ie
1350  land_temp(:,i,j) = lnd_temp
1351  land_water(:,i,j) = lnd_water
1352  land_sfc_temp(i,j) = sfc_temp
1353  land_sfc_albedo(i,j,i_lw) = sfc_albedo_lw
1354  land_sfc_albedo(i,j,i_sw) = sfc_albedo_sw
1355  end do
1356  end do
1357 
1358  return
1359  end subroutine land_setup
1360 
1361  !-----------------------------------------------------------------------------
1363  subroutine ocean_setup
1364  use mod_ocean_vars, only: &
1365  ocean_temp, &
1366  ocean_sfc_temp, &
1367  ocean_sfc_albedo, &
1368  ocean_sfc_z0m, &
1369  ocean_sfc_z0h, &
1371  implicit none
1372  ! Ocean state
1373  real(RP) :: OCN_TEMP ! ocean temperature [K]
1374  real(RP) :: SFC_TEMP ! ocean skin temperature [K]
1375  real(RP) :: SFC_albedo_LW = 0.04_rp ! ocean surface albedo for LW [0-1]
1376  real(RP) :: SFC_albedo_SW = 0.05_rp ! ocean surface albedo for SW [0-1]
1377  real(RP) :: SFC_Z0M = 1.0e-4_rp ! ocean surface roughness length (momentum) [m]
1378  real(RP) :: SFC_Z0H = 1.0e-4_rp ! ocean surface roughness length (heat) [m]
1379  real(RP) :: SFC_Z0E = 1.0e-4_rp ! ocean surface roughness length (vapor) [m]
1380 
1381  integer :: i, j
1382  integer :: ierr
1383 
1384  namelist /param_mkinit_ocean/ &
1385  ocn_temp, &
1386  sfc_temp, &
1387  sfc_albedo_lw, &
1388  sfc_albedo_sw, &
1389  sfc_z0m, &
1390  sfc_z0h, &
1391  sfc_z0e
1392 
1393  ocn_temp = thetastd
1394  sfc_temp = thetastd
1395 
1396  !--- read namelist
1397  rewind(io_fid_conf)
1398  read(io_fid_conf,nml=param_mkinit_ocean,iostat=ierr)
1399  if( ierr < 0 ) then !--- missing
1400  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1401  elseif( ierr > 0 ) then !--- fatal error
1402  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_OCEAN. Check!'
1403  call prc_mpistop
1404  endif
1405  if( io_lnml ) write(io_fid_log,nml=param_mkinit_ocean)
1406 
1407 
1408  do j = js, je
1409  do i = is, ie
1410  ocean_temp(i,j) = ocn_temp
1411  ocean_sfc_temp(i,j) = sfc_temp
1412  ocean_sfc_albedo(i,j,i_lw) = sfc_albedo_lw
1413  ocean_sfc_albedo(i,j,i_sw) = sfc_albedo_sw
1414  ocean_sfc_z0m(i,j) = sfc_z0m
1415  ocean_sfc_z0h(i,j) = sfc_z0h
1416  ocean_sfc_z0e(i,j) = sfc_z0e
1417  end do
1418  end do
1419 
1420  return
1421  end subroutine ocean_setup
1422 
1423  !-----------------------------------------------------------------------------
1425  subroutine urban_setup
1426  use mod_urban_vars, only: &
1427  urban_tr, &
1428  urban_tb, &
1429  urban_tg, &
1430  urban_tc, &
1431  urban_qc, &
1432  urban_uc, &
1433  urban_trl, &
1434  urban_tbl, &
1435  urban_tgl, &
1436  urban_rainr, &
1437  urban_rainb, &
1438  urban_raing, &
1439  urban_roff, &
1440  urban_sfc_temp, &
1442  implicit none
1443  ! urban state
1444  real(RP) :: URB_ROOF_TEMP ! Surface temperature of roof [K]
1445  real(RP) :: URB_BLDG_TEMP ! Surface temperature of building [K
1446  real(RP) :: URB_GRND_TEMP ! Surface temperature of ground [K]
1447  real(RP) :: URB_CNPY_TEMP ! Diagnostic canopy air temperature
1448  real(RP) :: URB_CNPY_HMDT = 0.0_rp ! Diagnostic canopy humidity [-]
1449  real(RP) :: URB_CNPY_WIND = 0.0_rp ! Diagnostic canopy wind [m/s]
1450  real(RP) :: URB_ROOF_LAYER_TEMP ! temperature in layer of roof [K]
1451  real(RP) :: URB_BLDG_LAYER_TEMP ! temperature in layer of building [
1452  real(RP) :: URB_GRND_LAYER_TEMP ! temperature in layer of ground [K]
1453  real(RP) :: URB_ROOF_RAIN = 0.0_rp ! temperature in layer of roof [K]
1454  real(RP) :: URB_BLDG_RAIN = 0.0_rp ! temperature in layer of building [
1455  real(RP) :: URB_GRND_RAIN = 0.0_rp ! temperature in layer of ground [K]
1456  real(RP) :: URB_RUNOFF = 0.0_rp ! temperature in layer of ground [K]
1457  real(RP) :: URB_SFC_TEMP ! Grid average of surface temperature [K]
1458  real(RP) :: URB_ALB_LW = 0.0_rp ! Grid average of surface albedo for LW [0-1]
1459  real(RP) :: URB_ALB_SW = 0.0_rp ! Grid average of surface albedo for SW [0-1]
1460 
1461  integer :: i, j
1462  integer :: ierr
1463 
1464  namelist /param_mkinit_urban/ &
1465  urb_roof_temp, &
1466  urb_bldg_temp, &
1467  urb_grnd_temp, &
1468  urb_cnpy_temp, &
1469  urb_cnpy_hmdt, &
1470  urb_cnpy_wind, &
1471  urb_roof_layer_temp, &
1472  urb_bldg_layer_temp, &
1473  urb_grnd_layer_temp, &
1474  urb_roof_rain, &
1475  urb_bldg_rain, &
1476  urb_grnd_rain, &
1477  urb_runoff, &
1478  urb_sfc_temp, &
1479  urb_alb_lw, &
1480  urb_alb_sw
1481 
1482  urb_roof_temp = thetastd
1483  urb_bldg_temp = thetastd
1484  urb_grnd_temp = thetastd
1485  urb_cnpy_temp = thetastd
1486  urb_roof_layer_temp = thetastd
1487  urb_bldg_layer_temp = thetastd
1488  urb_grnd_layer_temp = thetastd
1489 
1490  urb_sfc_temp = thetastd
1491 
1492  !--- read namelist
1493  rewind(io_fid_conf)
1494  read(io_fid_conf,nml=param_mkinit_urban,iostat=ierr)
1495  if( ierr < 0 ) then !--- missing
1496  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1497  elseif( ierr > 0 ) then !--- fatal error
1498  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_URBAN. Check!'
1499  call prc_mpistop
1500  endif
1501  if( io_lnml ) write(io_fid_log,nml=param_mkinit_urban)
1502 
1503 
1504  do j = js, je
1505  do i = is, ie
1506  urban_tr(i,j) = urb_roof_temp
1507  urban_tb(i,j) = urb_bldg_temp
1508  urban_tg(i,j) = urb_grnd_temp
1509  urban_tc(i,j) = urb_cnpy_temp
1510  urban_qc(i,j) = urb_cnpy_hmdt
1511  urban_uc(i,j) = urb_cnpy_wind
1512  urban_trl(:,i,j) = urb_roof_layer_temp
1513  urban_tbl(:,i,j) = urb_bldg_layer_temp
1514  urban_tgl(:,i,j) = urb_grnd_layer_temp
1515  urban_rainr(i,j) = urb_roof_rain
1516  urban_rainb(i,j) = urb_bldg_rain
1517  urban_raing(i,j) = urb_grnd_rain
1518  urban_roff(i,j) = urb_runoff
1519  urban_sfc_temp(i,j) = urb_sfc_temp
1520  urban_sfc_albedo(i,j,i_lw) = urb_alb_lw
1521  urban_sfc_albedo(i,j,i_sw) = urb_alb_sw
1522  end do
1523  end do
1524 
1525  return
1526  end subroutine urban_setup
1527 
1528  !-----------------------------------------------------------------------------
1530  subroutine read_sounding( &
1531  DENS, VELX, VELY, POTT, QV )
1532  implicit none
1533  real(RP), intent(out) :: DENS(ka)
1534  real(RP), intent(out) :: VELX(ka)
1535  real(RP), intent(out) :: VELY(ka)
1536  real(RP), intent(out) :: POTT(ka)
1537  real(RP), intent(out) :: QV (ka)
1538 
1539  real(RP) :: TEMP(ka)
1540  real(RP) :: PRES(ka)
1541  real(RP) :: QC (ka)
1542 
1543  character(len=H_LONG) :: ENV_IN_SOUNDING_file = ''
1544 
1545  integer, parameter :: EXP_klim = 100
1546  integer :: EXP_kmax
1547 
1548  real(RP) :: SFC_THETA ! surface potential temperature [K]
1549  real(RP) :: SFC_PRES ! surface pressure [hPa]
1550  real(RP) :: SFC_QV ! surface watervapor [g/kg]
1551 
1552  real(RP) :: EXP_z (exp_klim+1) ! height [m]
1553  real(RP) :: EXP_pott(exp_klim+1) ! potential temperature [K]
1554  real(RP) :: EXP_qv (exp_klim+1) ! water vapor [g/kg]
1555  real(RP) :: EXP_u (exp_klim+1) ! velocity u [m/s]
1556  real(RP) :: EXP_v (exp_klim+1) ! velocity v [m/s]
1557 
1558  real(RP) :: fact1, fact2
1559  integer :: k, kref
1560  integer :: fid
1561  integer :: ierr
1562 
1563  namelist /param_mkinit_sounding/ &
1564  env_in_sounding_file
1565 
1566  !--- read namelist
1567  rewind(io_fid_conf)
1568  read(io_fid_conf,nml=param_mkinit_sounding,iostat=ierr)
1569 
1570  if( ierr < 0 ) then !--- missing
1571  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1572  elseif( ierr > 0 ) then !--- fatal error
1573  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_SOUNDING. Check!'
1574  call prc_mpistop
1575  endif
1576  if( io_lnml ) write(io_fid_log,nml=param_mkinit_sounding)
1577 
1578  !--- prepare sounding profile
1579  if( io_l ) write(io_fid_log,*) '+++ Input sounding file:', trim(env_in_sounding_file)
1580  fid = io_get_available_fid()
1581  open( fid, &
1582  file = trim(env_in_sounding_file), &
1583  form = 'formatted', &
1584  status = 'old', &
1585  iostat = ierr )
1586 
1587  if ( ierr /= 0 ) then
1588  if( io_l ) write(*,*) 'xxx Input file not found!'
1589  endif
1590 
1591  !--- read sounding file till end
1592  read(fid,*) sfc_pres, sfc_theta, sfc_qv
1593 
1594  if( io_l ) write(io_fid_log,*) '+++ Surface pressure [hPa]', sfc_pres
1595  if( io_l ) write(io_fid_log,*) '+++ Surface pot. temp [K]', sfc_theta
1596  if( io_l ) write(io_fid_log,*) '+++ Surface water vapor [g/kg]', sfc_qv
1597 
1598  do k = 2, exp_klim
1599  read(fid,*,iostat=ierr) exp_z(k), exp_pott(k), exp_qv(k), exp_u(k), exp_v(k)
1600  if ( ierr /= 0 ) exit
1601  enddo
1602 
1603  exp_kmax = k - 1
1604  close(fid)
1605 
1606  ! Boundary
1607  exp_z(1) = 0.0_rp
1608  exp_pott(1) = sfc_theta
1609  exp_qv(1) = sfc_qv
1610  exp_u(1) = exp_u(2)
1611  exp_v(1) = exp_v(2)
1612  exp_z(exp_kmax+1) = 100.e3_rp
1613  exp_pott(exp_kmax+1) = exp_pott(exp_kmax)
1614  exp_qv(exp_kmax+1) = exp_qv(exp_kmax)
1615  exp_u(exp_kmax+1) = exp_u(exp_kmax)
1616  exp_v(exp_kmax+1) = exp_v(exp_kmax)
1617 
1618  do k = 1, exp_kmax+1
1619  exp_qv(k) = exp_qv(k) * 1.e-3_rp ! [g/kg]->[kg/kg]
1620  enddo
1621 
1622  ! calc in dry condition
1623  pres_sfc = sfc_pres * 1.e2_rp ! [hPa]->[Pa]
1624  pott_sfc = sfc_theta
1625  qv_sfc = sfc_qv * 1.e-3_rp ! [g/kg]->[kg/kg]
1626  qc_sfc = 0.0_rp
1627 
1628  !--- linear interpolate to model grid
1629  do k = ks, ke
1630  qc(k) = 0.0_rp
1631 
1632  do kref = 2, exp_kmax+1
1633  if ( grid_cz(k) > exp_z(kref-1) &
1634  .AND. grid_cz(k) <= exp_z(kref ) ) then
1635 
1636  fact1 = ( exp_z(kref) - grid_cz(k) ) / ( exp_z(kref)-exp_z(kref-1) )
1637  fact2 = ( grid_cz(k) - exp_z(kref-1) ) / ( exp_z(kref)-exp_z(kref-1) )
1638 
1639  pott(k) = exp_pott(kref-1) * fact1 &
1640  + exp_pott(kref ) * fact2
1641  qv(k) = exp_qv(kref-1) * fact1 &
1642  + exp_qv(kref ) * fact2
1643  velx(k) = exp_u(kref-1) * fact1 &
1644  + exp_u(kref ) * fact2
1645  vely(k) = exp_v(kref-1) * fact1 &
1646  + exp_v(kref ) * fact2
1647  endif
1648  enddo
1649  enddo
1650 
1651  ! make density & pressure profile in moist condition
1652  call hydrostatic_buildrho( dens(:), & ! [OUT]
1653  temp(:), & ! [OUT]
1654  pres(:), & ! [OUT]
1655  pott(:), & ! [IN]
1656  qv(:), & ! [IN]
1657  qc(:), & ! [IN]
1658  temp_sfc(1,1,1), & ! [OUT]
1659  pres_sfc(1,1,1), & ! [IN]
1660  pott_sfc(1,1,1), & ! [IN]
1661  qv_sfc(1,1,1), & ! [IN]
1662  qc_sfc(1,1,1) ) ! [IN]
1663 
1664  return
1665  end subroutine read_sounding
1666 
1667  !-----------------------------------------------------------------------------
1669  subroutine mkinit_planestate
1670  implicit none
1671 
1672  ! Surface state
1673  real(RP) :: SFC_THETA ! surface potential temperature [K]
1674  real(RP) :: SFC_PRES ! surface pressure [Pa]
1675  real(RP) :: SFC_RH = 0.0_rp ! surface relative humidity [%]
1676  ! Environment state
1677  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1678  real(RP) :: ENV_TLAPS = 0.0_rp ! Lapse rate of THETA [K/m]
1679  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
1680  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1681  real(RP) :: ENV_RH = 0.0_rp ! relative humidity of environment [%]
1682  ! Disturbance
1683  real(RP) :: RANDOM_THETA = 0.0_rp ! amplitude of random disturbance theta
1684  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
1685  real(RP) :: RANDOM_V = 0.0_rp ! amplitude of random disturbance v
1686  real(RP) :: RANDOM_RH = 0.0_rp ! amplitude of random disturbance RH
1687 
1688  namelist / param_mkinit_planestate / &
1689  sfc_theta, &
1690  sfc_pres, &
1691  sfc_rh, &
1692  env_theta, &
1693  env_tlaps, &
1694  env_u, &
1695  env_v, &
1696  env_rh, &
1697  random_theta, &
1698  random_u, &
1699  random_v, &
1700  random_rh
1701 
1702  integer :: ierr
1703  integer :: k, i, j
1704  !---------------------------------------------------------------------------
1705 
1706  if( io_l ) write(io_fid_log,*)
1707  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit Horiz_UNIFORM] / Categ[preprocess] / Origin[SCALE-RM]'
1708 
1709  sfc_theta = thetastd
1710  sfc_pres = pstd
1711  env_theta = thetastd
1712 
1713  !--- read namelist
1714  rewind(io_fid_conf)
1715  read(io_fid_conf,nml=param_mkinit_planestate,iostat=ierr)
1716 
1717  if( ierr < 0 ) then !--- missing
1718  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1719  elseif( ierr > 0 ) then !--- fatal error
1720  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_PLANESTATE. Check!'
1721  call prc_mpistop
1722  endif
1723  if( io_lnml ) write(io_fid_log,nml=param_mkinit_planestate)
1724 
1725  ! calc in dry condition
1726  do j = js, je
1727  do i = is, ie
1728  pott_sfc(1,i,j) = sfc_theta
1729  pres_sfc(1,i,j) = sfc_pres
1730  qv_sfc(1,i,j) = 0.0_rp
1731  qc_sfc(1,i,j) = 0.0_rp
1732 
1733  do k = ks, ke
1734  qv(k,i,j) = 0.0_rp
1735  qc(k,i,j) = 0.0_rp
1736  enddo
1737  enddo
1738  enddo
1739 
1740  if ( env_theta < 0.0_rp ) then ! use isa profile
1741 
1742  call profile_isa( ka, ks, ke, & ! [IN]
1743  ia, is, ie, & ! [IN]
1744  ja, js, je, & ! [IN]
1745  pott_sfc(1,:,:), & ! [IN]
1746  pres_sfc(1,:,:), & ! [IN]
1747  real_cz(:,:,:), & ! [IN]
1748  pott(:,:,:) ) ! [OUT]
1749 
1750  else
1751 
1752  do j = js, je
1753  do i = is, ie
1754  do k = ks, ke
1755  pott(k,i,j) = env_theta + env_tlaps * real_cz(k,i,j)
1756  enddo
1757  enddo
1758  enddo
1759 
1760  endif
1761 
1762  ! make density & pressure profile in moist condition
1763  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
1764  temp(:,:,:), & ! [OUT]
1765  pres(:,:,:), & ! [OUT]
1766  pott(:,:,:), & ! [IN]
1767  qv(:,:,:), & ! [IN]
1768  qc(:,:,:), & ! [IN]
1769  temp_sfc(:,:,:), & ! [OUT]
1770  pres_sfc(:,:,:), & ! [IN]
1771  pott_sfc(:,:,:), & ! [IN]
1772  qv_sfc(:,:,:), & ! [IN]
1773  qc_sfc(:,:,:) ) ! [IN]
1774 
1775  ! calc QV from RH
1776  call saturation_pres2qsat_all( qsat_sfc(1,:,:), temp_sfc(1,:,:), pres_sfc(1,:,:) )
1777  call saturation_pres2qsat_all( qsat(:,:,:), temp(:,:,:), pres(:,:,:) )
1778 
1779  call random_get(rndm) ! make random
1780  do j = js, je
1781  do i = is, ie
1782  qv_sfc(1,i,j) = ( sfc_rh + rndm(ks-1,i,j) * random_rh ) * 1.e-2_rp * qsat_sfc(1,i,j)
1783 
1784  do k = ks, ke
1785  qv(k,i,j) = ( env_rh + rndm(k,i,j) * random_rh ) * 1.e-2_rp * qsat(k,i,j)
1786  enddo
1787  enddo
1788  enddo
1789 
1790  call random_get(rndm) ! make random
1791  do j = js, je
1792  do i = is, ie
1793  pott_sfc(1,i,j) = pott_sfc(1,i,j) + rndm(ks-1,i,j) * random_theta
1794 
1795  do k = ks, ke
1796  pott(k,i,j) = pott(k,i,j) + rndm(k,i,j) * random_theta
1797  enddo
1798  enddo
1799  enddo
1800 
1801  ! make density & pressure profile in moist condition
1802  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
1803  temp(:,:,:), & ! [OUT]
1804  pres(:,:,:), & ! [OUT]
1805  pott(:,:,:), & ! [IN]
1806  qv(:,:,:), & ! [IN]
1807  qc(:,:,:), & ! [IN]
1808  temp_sfc(:,:,:), & ! [OUT]
1809  pres_sfc(:,:,:), & ! [IN]
1810  pott_sfc(:,:,:), & ! [IN]
1811  qv_sfc(:,:,:), & ! [IN]
1812  qc_sfc(:,:,:) ) ! [IN]
1813 
1814  call comm_vars8( dens(:,:,:), 1 )
1815  call comm_wait ( dens(:,:,:), 1 )
1816 
1817  call random_get(rndm) ! make random
1818  do j = js, je
1819  do i = is, ie
1820  do k = ks, ke
1821  momx(k,i,j) = ( env_u + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u ) &
1822  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
1823  enddo
1824  enddo
1825  enddo
1826 
1827  call random_get(rndm) ! make random
1828  do j = js, je
1829  do i = is, ie
1830  do k = ks, ke
1831  momy(k,i,j) = ( env_v + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_v ) &
1832  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
1833  enddo
1834  enddo
1835  enddo
1836 
1837  do j = js, je
1838  do i = is, ie
1839  do k = ks, ke
1840  momz(k,i,j) = 0.0_rp
1841  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
1842 
1843  qtrc(k,i,j,i_qv) = qv(k,i,j)
1844  enddo
1845  enddo
1846  enddo
1847 
1848  call flux_setup
1849 
1850  call ocean_setup
1851 
1852  return
1853  end subroutine mkinit_planestate
1854 
1855  !-----------------------------------------------------------------------------
1857  subroutine mkinit_tracerbubble
1858  implicit none
1859 
1860 #ifndef DRY
1861  ! Surface state
1862  real(RP) :: SFC_THETA ! surface potential temperature [K]
1863  real(RP) :: SFC_PRES ! surface pressure [Pa]
1864  ! Environment state
1865  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1866  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
1867  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1868  ! Bubble
1869  character(len=H_SHORT) :: SHAPE_NC = 'BUBBLE' ! BUBBLE or RECT
1870  real(RP) :: BBL_NC = 1.0_rp ! extremum of NC in bubble [kg/kg]
1871 
1872  namelist / param_mkinit_tracerbubble / &
1873  sfc_theta, &
1874  sfc_pres, &
1875  env_theta, &
1876  env_u, &
1877  env_v, &
1878  shape_nc, &
1879  bbl_nc
1880 
1881  real(RP), pointer :: shapeFac(:,:,:) => null()
1882 
1883  integer :: k, i, j, iq
1884  integer :: ierr
1885  !---------------------------------------------------------------------------
1886 
1887  if( io_l ) write(io_fid_log,*)
1888  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit TRACERBUBBLE] / Categ[preprocess] / Origin[SCALE-RM]'
1889 
1890  sfc_theta = thetastd
1891  sfc_pres = pstd
1892  env_theta = thetastd
1893 
1894  !--- read namelist
1895  rewind(io_fid_conf)
1896  read(io_fid_conf,nml=param_mkinit_tracerbubble,iostat=ierr)
1897 
1898  if( ierr < 0 ) then !--- missing
1899  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
1900  elseif( ierr > 0 ) then !--- fatal error
1901  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_TRACERBUBBLE. Check!'
1902  call prc_mpistop
1903  endif
1904  if( io_lnml ) write(io_fid_log,nml=param_mkinit_tracerbubble)
1905 
1906  ! calc in dry condition
1907  pres_sfc(1,1,1) = sfc_pres
1908  pott_sfc(1,1,1) = sfc_theta
1909  qv_sfc(1,1,1) = 0.0_rp
1910  qc_sfc(1,1,1) = 0.0_rp
1911 
1912  do k = ks, ke
1913  pott(k,1,1) = env_theta
1914  qv(k,1,1) = 0.0_rp
1915  qc(k,1,1) = 0.0_rp
1916  enddo
1917 
1918  ! make density & pressure profile in dry condition
1919  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
1920  temp(:,1,1), & ! [OUT]
1921  pres(:,1,1), & ! [OUT]
1922  pott(:,1,1), & ! [IN]
1923  qv(:,1,1), & ! [IN]
1924  qc(:,1,1), & ! [IN]
1925  temp_sfc(1,1,1), & ! [OUT]
1926  pres_sfc(1,1,1), & ! [IN]
1927  pott_sfc(1,1,1), & ! [IN]
1928  qv_sfc(1,1,1), & ! [IN]
1929  qc_sfc(1,1,1) ) ! [IN]
1930 
1931  do j = js, je
1932  do i = is, ie
1933  do k = ks, ke
1934  dens(k,i,j) = dens(k,1,1)
1935  momz(k,i,j) = 0.0_rp
1936  momx(k,i,j) = env_u * dens(k,1,1)
1937  momy(k,i,j) = env_v * dens(k,1,1)
1938  rhot(k,i,j) = pott(k,1,1) * dens(k,1,1)
1939 
1940  do iq = 1, qa
1941  qtrc(k,i,j,iq) = 0.0_rp
1942  enddo
1943  enddo
1944  enddo
1945  enddo
1946 
1947  ! make tracer bubble
1948  if ( i_nc > 0 ) then
1949 
1950  select case(shape_nc)
1951  case('BUBBLE')
1952  call bubble_setup
1953  shapefac => bubble
1954  case('RECT')
1955  call rect_setup
1956  shapefac => rect
1957  case default
1958  write(*,*) 'xxx SHAPE_NC=', trim(shape_nc), ' cannot be used on advect. Check!'
1959  call prc_mpistop
1960  end select
1961 
1962  do j = js, je
1963  do i = is, ie
1964  do k = ks, ke
1965  qtrc(k,i,j,i_nc) = bbl_nc * shapefac(k,i,j)
1966  enddo
1967  enddo
1968  enddo
1969  else
1970  write(*,*) 'xxx tracer I_NC is not defined. Check!'
1971  call prc_mpistop
1972  endif
1973 
1974  if ( tracer_type == 'SUZUKI10' ) then
1975  write(*,*) 'xxx SBM cannot be used on tracerbubble. Check!'
1976  call prc_mpistop
1977  endif
1978 #endif
1979 
1980  return
1981  end subroutine mkinit_tracerbubble
1982 
1983  !-----------------------------------------------------------------------------
1993  subroutine mkinit_coldbubble
1994  implicit none
1995 
1996  ! Surface state
1997  real(RP) :: SFC_THETA ! surface potential temperature [K]
1998  real(RP) :: SFC_PRES ! surface pressure [Pa]
1999  ! Environment state
2000  real(RP) :: ENV_THETA ! potential temperature of environment [K]
2001  ! Bubble
2002  real(RP) :: BBL_TEMP = -15.0_rp ! extremum of temperature in bubble [K]
2003 
2004  namelist / param_mkinit_coldbubble / &
2005  sfc_theta, &
2006  sfc_pres, &
2007  env_theta, &
2008  bbl_temp
2009 
2010  real(RP) :: RovCP
2011 
2012  integer :: ierr
2013  integer :: k, i, j, iq
2014  !---------------------------------------------------------------------------
2015 
2016  if( io_l ) write(io_fid_log,*)
2017  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit COLDBUBBLE] / Categ[preprocess] / Origin[SCALE-RM]'
2018 
2019  sfc_theta = thetastd
2020  sfc_pres = pstd
2021  env_theta = thetastd
2022 
2023  !--- read namelist
2024  rewind(io_fid_conf)
2025  read(io_fid_conf,nml=param_mkinit_coldbubble,iostat=ierr)
2026 
2027  if( ierr < 0 ) then !--- missing
2028  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2029  elseif( ierr > 0 ) then !--- fatal error
2030  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_COLDBUBBLE. Check!'
2031  call prc_mpistop
2032  endif
2033  if( io_lnml ) write(io_fid_log,nml=param_mkinit_coldbubble)
2034 
2035  rovcp = rdry / cpdry
2036 
2037  ! calc in dry condition
2038  pres_sfc(1,1,1) = sfc_pres
2039  pott_sfc(1,1,1) = sfc_theta
2040  qv_sfc(1,1,1) = 0.0_rp
2041  qc_sfc(1,1,1) = 0.0_rp
2042 
2043  do k = ks, ke
2044  pott(k,1,1) = env_theta
2045  qv(k,1,1) = 0.0_rp
2046  qc(k,1,1) = 0.0_rp
2047  enddo
2048 
2049  ! make density & pressure profile in dry condition
2050  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2051  temp(:,1,1), & ! [OUT]
2052  pres(:,1,1), & ! [OUT]
2053  pott(:,1,1), & ! [IN]
2054  qv(:,1,1), & ! [IN]
2055  qc(:,1,1), & ! [IN]
2056  temp_sfc(1,1,1), & ! [OUT]
2057  pres_sfc(1,1,1), & ! [IN]
2058  pott_sfc(1,1,1), & ! [IN]
2059  qv_sfc(1,1,1), & ! [IN]
2060  qc_sfc(1,1,1) ) ! [IN]
2061 
2062  do j = 1, ja
2063  do i = 1, ia
2064  do k = ks, ke
2065  dens(k,i,j) = dens(k,1,1)
2066  momz(k,i,j) = 0.0_rp
2067  momx(k,i,j) = 0.0_rp
2068  momy(k,i,j) = 0.0_rp
2069 
2070  ! make cold bubble
2071  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) &
2072  + bbl_temp * ( p00/pres(k,1,1) )**rovcp * bubble(k,i,j) )
2073 
2074  do iq = 1, qa
2075  qtrc(k,i,j,iq) = 0.0_rp
2076  enddo
2077  enddo
2078  enddo
2079  enddo
2080 
2081  if ( tracer_type == 'SUZUKI10' ) then
2082  write(*,*) 'xxx SBM cannot be used on coldbubble. Check!'
2083  call prc_mpistop
2084  endif
2085 
2086  return
2087  end subroutine mkinit_coldbubble
2088 
2089  !-----------------------------------------------------------------------------
2091  subroutine mkinit_lambwave
2092  implicit none
2093 
2094  ! Surface state
2095  real(RP) :: SFC_PRES ! surface pressure [Pa]
2096  ! Environment state
2097  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
2098  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2099  real(RP) :: ENV_TEMP = 300.0_rp ! temperature of environment [K]
2100  ! Bubble
2101  real(RP) :: BBL_PRES = 100._rp ! extremum of pressure in bubble [Pa]
2102 
2103  namelist / param_mkinit_lambwave / &
2104  sfc_pres, &
2105  env_u, &
2106  env_v, &
2107  env_temp, &
2108  bbl_pres
2109 
2110  real(RP) :: RovCP
2111 
2112  integer :: ierr
2113  integer :: k, i, j, iq
2114  !---------------------------------------------------------------------------
2115 
2116  if( io_l ) write(io_fid_log,*)
2117  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit LAMBWAVE] / Categ[preprocess] / Origin[SCALE-RM]'
2118 
2119  sfc_pres = pstd
2120 
2121  !--- read namelist
2122  rewind(io_fid_conf)
2123  read(io_fid_conf,nml=param_mkinit_lambwave,iostat=ierr)
2124 
2125  if( ierr < 0 ) then !--- missing
2126  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2127  elseif( ierr > 0 ) then !--- fatal error
2128  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_LAMBWAVE. Check!'
2129  call prc_mpistop
2130  endif
2131  if( io_lnml ) write(io_fid_log,nml=param_mkinit_lambwave)
2132 
2133  rovcp = rdry / cpdry
2134 
2135  do j = js, je
2136  do i = is, ie
2137  do k = ks, ke
2138  dens(k,i,j) = sfc_pres/(rdry*env_temp) * exp( - grav/(rdry*env_temp) * grid_cz(k) )
2139  momz(k,i,j) = 0.0_rp
2140  momx(k,i,j) = env_u * dens(k,i,j)
2141  momy(k,i,j) = env_v * dens(k,i,j)
2142 
2143  ! make pressure bubble
2144  pres(k,i,j) = dens(k,i,j) * env_temp * rdry + bbl_pres * bubble(k,i,j)
2145 
2146  rhot(k,i,j) = dens(k,i,j) * env_temp * ( p00/pres(k,i,j) )**rovcp
2147 
2148  do iq = 1, qa
2149  qtrc(k,i,j,iq) = 0.0_rp
2150  enddo
2151  enddo
2152  enddo
2153  enddo
2154 
2155  if ( tracer_type == 'SUZUKI10' ) then
2156  write(*,*) 'xxx SBM cannot be used on lambwave. Check!'
2157  call prc_mpistop
2158  endif
2159 
2160  return
2161  end subroutine mkinit_lambwave
2162 
2163  !-----------------------------------------------------------------------------
2166  subroutine mkinit_gravitywave
2167  implicit none
2168 
2169  ! Surface state
2170  real(RP) :: SFC_THETA ! surface potential temperature [K]
2171  real(RP) :: SFC_PRES ! surface pressure [Pa]
2172  ! Environment state
2173  real(RP) :: ENV_U = 20.0_rp ! velocity u of environment [m/s]
2174  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2175  real(RP) :: ENV_BVF = 0.01_rp ! Brunt Vaisala frequencies of environment [1/s]
2176  ! Bubble
2177  real(RP) :: BBL_THETA = 0.01_rp ! extremum of potential temperature in bubble [K]
2178 
2179  namelist / param_mkinit_gravitywave / &
2180  sfc_theta, &
2181  sfc_pres, &
2182  env_u, &
2183  env_v, &
2184  env_bvf, &
2185  bbl_theta
2186 
2187  integer :: ierr
2188  integer :: k, i, j, iq
2189  !---------------------------------------------------------------------------
2190 
2191  if( io_l ) write(io_fid_log,*)
2192  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit GRAVITYWAVE] / Categ[preprocess] / Origin[SCALE-RM]'
2193 
2194  sfc_theta = thetastd
2195  sfc_pres = pstd
2196 
2197  !--- read namelist
2198  rewind(io_fid_conf)
2199  read(io_fid_conf,nml=param_mkinit_gravitywave,iostat=ierr)
2200 
2201  if( ierr < 0 ) then !--- missing
2202  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2203  elseif( ierr > 0 ) then !--- fatal error
2204  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_GRAVITYWAVE. Check!'
2205  call prc_mpistop
2206  endif
2207  if( io_lnml ) write(io_fid_log,nml=param_mkinit_gravitywave)
2208 
2209  ! calc in dry condition
2210  pres_sfc(1,1,1) = sfc_pres
2211  pott_sfc(1,1,1) = sfc_theta
2212  qv_sfc(1,1,1) = 0.0_rp
2213  qc_sfc(1,1,1) = 0.0_rp
2214 
2215  do k = ks, ke
2216  pott(k,1,1) = sfc_theta * exp( env_bvf*env_bvf / grav * grid_cz(k) )
2217  qv(k,1,1) = 0.0_rp
2218  qc(k,1,1) = 0.0_rp
2219  enddo
2220 
2221  ! make density & pressure profile in dry condition
2222  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2223  temp(:,1,1), & ! [OUT]
2224  pres(:,1,1), & ! [OUT]
2225  pott(:,1,1), & ! [IN]
2226  qv(:,1,1), & ! [IN]
2227  qc(:,1,1), & ! [IN]
2228  temp_sfc(1,1,1), & ! [OUT]
2229  pres_sfc(1,1,1), & ! [IN]
2230  pott_sfc(1,1,1), & ! [IN]
2231  qv_sfc(1,1,1), & ! [IN]
2232  qc_sfc(1,1,1) ) ! [IN]
2233 
2234  do j = js, je
2235  do i = is, ie
2236  do k = ks, ke
2237  dens(k,i,j) = dens(k,1,1)
2238  momz(k,i,j) = 0.0_rp
2239  momx(k,i,j) = env_u * dens(k,1,1)
2240  momy(k,i,j) = env_v * dens(k,1,1)
2241 
2242  ! make warm bubble
2243  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
2244 
2245  do iq = 1, qa
2246  qtrc(k,i,j,iq) = 0.0_rp
2247  enddo
2248  enddo
2249  enddo
2250  enddo
2251 
2252  if ( tracer_type == 'SUZUKI10' ) then
2253  write(*,*) 'xxx SBM cannot be used on gravitywave. Check!'
2254  call prc_mpistop
2255  endif
2256 
2257  return
2258  end subroutine mkinit_gravitywave
2259 
2260  !-----------------------------------------------------------------------------
2262  subroutine mkinit_khwave
2263  implicit none
2264 
2265  ! Surface state
2266  real(RP) :: SFC_THETA ! surface potential temperature [K]
2267  real(RP) :: SFC_PRES ! surface pressure [Pa]
2268  ! Environment state
2269  real(RP) :: ENV_L1_ZTOP = 1900.0_rp ! top height of the layer1 (low THETA) [m]
2270  real(RP) :: ENV_L3_ZBOTTOM = 2100.0_rp ! bottom height of the layer3 (high THETA) [m]
2271  real(RP) :: ENV_L1_THETA = 300.0_rp ! THETA in the layer1 (low THETA) [K]
2272  real(RP) :: ENV_L3_THETA = 301.0_rp ! THETA in the layer3 (high THETA) [K]
2273  real(RP) :: ENV_L1_U = 0.0_rp ! velocity u in the layer1 (low THETA) [K]
2274  real(RP) :: ENV_L3_U = 20.0_rp ! velocity u in the layer3 (high THETA) [K]
2275  ! Disturbance
2276  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
2277 
2278  namelist / param_mkinit_khwave / &
2279  sfc_theta, &
2280  sfc_pres, &
2281  env_l1_ztop, &
2282  env_l3_zbottom, &
2283  env_l1_theta, &
2284  env_l3_theta, &
2285  env_l1_u, &
2286  env_l3_u, &
2287  random_u
2288 
2289  real(RP) :: fact
2290 
2291  integer :: ierr
2292  integer :: k, i, j, iq
2293  !---------------------------------------------------------------------------
2294 
2295  if( io_l ) write(io_fid_log,*)
2296  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit KHWAVE] / Categ[preprocess] / Origin[SCALE-RM]'
2297 
2298  sfc_theta = thetastd
2299  sfc_pres = pstd
2300 
2301  !--- read namelist
2302  rewind(io_fid_conf)
2303  read(io_fid_conf,nml=param_mkinit_khwave,iostat=ierr)
2304 
2305  if( ierr < 0 ) then !--- missing
2306  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2307  elseif( ierr > 0 ) then !--- fatal error
2308  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_KHWAVE. Check!'
2309  call prc_mpistop
2310  endif
2311  if( io_lnml ) write(io_fid_log,nml=param_mkinit_khwave)
2312 
2313  ! calc in dry condition
2314  pres_sfc(1,1,1) = sfc_pres
2315  pott_sfc(1,1,1) = sfc_theta
2316  qv_sfc(1,1,1) = 0.0_rp
2317  qc_sfc(1,1,1) = 0.0_rp
2318 
2319  do k = ks, ke
2320  fact = ( grid_cz(k)-env_l1_ztop ) / ( env_l3_zbottom-env_l1_ztop )
2321  fact = max( min( fact, 1.0_rp ), 0.0_rp )
2322 
2323  pott(k,1,1) = env_l1_theta * ( 1.0_rp - fact ) &
2324  + env_l3_theta * ( fact )
2325 
2326  qv(k,1,1) = 0.0_rp
2327  qc(k,1,1) = 0.0_rp
2328  enddo
2329 
2330  ! make density & pressure profile in dry condition
2331  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2332  temp(:,1,1), & ! [OUT]
2333  pres(:,1,1), & ! [OUT]
2334  pott(:,1,1), & ! [IN]
2335  qv(:,1,1), & ! [IN]
2336  qc(:,1,1), & ! [IN]
2337  temp_sfc(1,1,1), & ! [OUT]
2338  pres_sfc(1,1,1), & ! [IN]
2339  pott_sfc(1,1,1), & ! [IN]
2340  qv_sfc(1,1,1), & ! [IN]
2341  qc_sfc(1,1,1) ) ! [IN]
2342 
2343  ! calc QV from RH
2344  call saturation_pres2qsat_all( qsat_sfc(1,1,1), temp_sfc(1,1,1), pres_sfc(1,1,1) )
2345  call saturation_pres2qsat_all( qsat(:,1,1), temp(:,1,1), pres(:,1,1) )
2346 
2347  do j = js, je
2348  do i = is, ie
2349  do k = ks, ke
2350  dens(k,i,j) = dens(k,1,1)
2351  momz(k,i,j) = 0.0_rp
2352  momy(k,i,j) = 0.0_rp
2353  rhot(k,i,j) = dens(k,1,1) * pott(k,1,1)
2354 
2355  do iq = 1, qa
2356  qtrc(k,i,j,iq) = 0.0_rp
2357  enddo
2358  enddo
2359  enddo
2360  enddo
2361 
2362  call random_get(rndm) ! make random
2363  do j = js, je
2364  do i = is, ie
2365  do k = ks, ke
2366  fact = ( grid_cz(k)-env_l1_ztop ) / ( env_l3_zbottom-env_l1_ztop )
2367  fact = max( min( fact, 1.0_rp ), 0.0_rp )
2368 
2369  momx(k,i,j) = ( env_l1_u * ( 1.0_rp - fact ) &
2370  + env_l3_u * ( fact ) &
2371  + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u &
2372  ) * dens(k,i,j)
2373  enddo
2374  enddo
2375  enddo
2376 
2377  if ( tracer_type == 'SUZUKI10' ) then
2378  write(*,*) 'xxx SBM cannot be used on khwave. Check!'
2379  call prc_mpistop
2380  endif
2381 
2382  return
2383  end subroutine mkinit_khwave
2384 
2385  !-----------------------------------------------------------------------------
2387  subroutine mkinit_turbulence
2388  implicit none
2389 
2390  ! Surface state
2391  real(RP) :: SFC_THETA ! surface potential temperature [K]
2392  real(RP) :: SFC_PRES ! surface pressure [Pa]
2393  real(RP) :: SFC_RH = 0.0_rp ! surface relative humidity [%]
2394  ! Environment state
2395  real(RP) :: ENV_THETA ! potential temperature of environment [K]
2396  real(RP) :: ENV_TLAPS = 4.e-3_rp ! Lapse rate of THETA [K/m]
2397  real(RP) :: ENV_U = 5.0_rp ! velocity u of environment [m/s]
2398  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2399  real(RP) :: ENV_RH = 0.0_rp ! relative humidity of environment [%]
2400  ! Disturbance
2401  real(RP) :: RANDOM_THETA = 1.0_rp ! amplitude of random disturbance theta
2402  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
2403  real(RP) :: RANDOM_V = 0.0_rp ! amplitude of random disturbance v
2404  real(RP) :: RANDOM_RH = 0.0_rp ! amplitude of random disturbance RH
2405 
2406  namelist / param_mkinit_turbulence / &
2407  sfc_theta, &
2408  sfc_pres, &
2409  sfc_rh, &
2410  env_theta, &
2411  env_tlaps, &
2412  env_u, &
2413  env_v, &
2414  env_rh, &
2415  random_theta, &
2416  random_u, &
2417  random_v, &
2418  random_rh
2419 
2420  integer :: ierr
2421  integer :: k, i, j, iq
2422  !---------------------------------------------------------------------------
2423 
2424  if( io_l ) write(io_fid_log,*)
2425  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit TURBULENCE] / Categ[preprocess] / Origin[SCALE-RM]'
2426 
2427  sfc_theta = thetastd
2428  sfc_pres = pstd
2429  env_theta = thetastd
2430 
2431  !--- read namelist
2432  rewind(io_fid_conf)
2433  read(io_fid_conf,nml=param_mkinit_turbulence,iostat=ierr)
2434 
2435  if( ierr < 0 ) then !--- missing
2436  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2437  elseif( ierr > 0 ) then !--- fatal error
2438  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_TURBULENCE. Check!'
2439  call prc_mpistop
2440  endif
2441  if( io_lnml ) write(io_fid_log,nml=param_mkinit_turbulence)
2442 
2443  ! calc in dry condition
2444  pres_sfc(1,1,1) = sfc_pres
2445  pott_sfc(1,1,1) = sfc_theta
2446  qv_sfc(1,1,1) = 0.0_rp
2447  qc_sfc(1,1,1) = 0.0_rp
2448 
2449  do k = ks, ke
2450  pott(k,1,1) = env_theta + env_tlaps * grid_cz(k)
2451  qv(k,1,1) = 0.0_rp
2452  qc(k,1,1) = 0.0_rp
2453  enddo
2454 
2455  ! make density & pressure profile in dry condition
2456  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2457  temp(:,1,1), & ! [OUT]
2458  pres(:,1,1), & ! [OUT]
2459  pott(:,1,1), & ! [IN]
2460  qv(:,1,1), & ! [IN]
2461  qc(:,1,1), & ! [IN]
2462  temp_sfc(1,1,1), & ! [OUT]
2463  pres_sfc(1,1,1), & ! [IN]
2464  pott_sfc(1,1,1), & ! [IN]
2465  qv_sfc(1,1,1), & ! [IN]
2466  qc_sfc(1,1,1) ) ! [IN]
2467 
2468  ! calc QV from RH
2469  call saturation_pres2qsat_all( qsat_sfc(1,1,1), temp_sfc(1,1,1), pres_sfc(1,1,1) )
2470  call saturation_pres2qsat_all( qsat(:,1,1), temp(:,1,1), pres(:,1,1) )
2471 
2472  call random_get(rndm) ! make random
2473  do j = js, je
2474  do i = is, ie
2475  qv_sfc(1,i,j) = ( sfc_rh + rndm(ks-1,i,j) * random_rh ) * 1.e-2_rp * qsat_sfc(1,1,1)
2476  qc_sfc(1,i,j) = 0.0_rp
2477 
2478  do k = ks, ke
2479  qv(k,i,j) = ( env_rh + rndm(k,i,j) * random_rh ) * 1.e-2_rp * qsat(k,1,1)
2480  qc(k,i,j) = 0.0_rp
2481  enddo
2482  enddo
2483  enddo
2484 
2485  call random_get(rndm) ! make random
2486  do j = js, je
2487  do i = is, ie
2488  pres_sfc(1,i,j) = sfc_pres
2489  pott_sfc(1,i,j) = sfc_theta + rndm(ks-1,i,j) * random_theta
2490 
2491  do k = ks, ke
2492  pott(k,i,j) = env_theta + env_tlaps * grid_cz(k) + rndm(k,i,j) * random_theta
2493  enddo
2494  enddo
2495  enddo
2496 
2497  ! make density & pressure profile in moist condition
2498  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
2499  temp(:,:,:), & ! [OUT]
2500  pres(:,:,:), & ! [OUT]
2501  pott(:,:,:), & ! [IN]
2502  qv(:,:,:), & ! [IN]
2503  qc(:,:,:), & ! [IN]
2504  temp_sfc(:,:,:), & ! [OUT]
2505  pres_sfc(:,:,:), & ! [IN]
2506  pott_sfc(:,:,:), & ! [IN]
2507  qv_sfc(:,:,:), & ! [IN]
2508  qc_sfc(:,:,:) ) ! [IN]
2509 
2510  call comm_vars8( dens(:,:,:), 1 )
2511  call comm_wait ( dens(:,:,:), 1 )
2512 
2513  call random_get(rndm) ! make random
2514  do j = js, je
2515  do i = is, ie
2516  do k = ks, ke
2517  momx(k,i,j) = ( env_u + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u ) &
2518  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
2519  enddo
2520  enddo
2521  enddo
2522 
2523  call random_get(rndm) ! make random
2524  do j = js, je
2525  do i = is, ie
2526  do k = ks, ke
2527  momy(k,i,j) = ( env_v + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_v ) &
2528  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
2529  enddo
2530  enddo
2531  enddo
2532 
2533  do j = js, je
2534  do i = is, ie
2535  do k = ks, ke
2536  momz(k,i,j) = 0.0_rp
2537  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
2538 
2539  qtrc(k,i,j,i_qv) = qv(k,i,j)
2540  enddo
2541  enddo
2542  enddo
2543 
2544 #ifndef DRY
2545  if ( qa >= 2 ) then
2546  do iq = 2, qa
2547  do j = js, je
2548  do i = is, ie
2549  do k = ks, ke
2550  qtrc(k,i,j,iq) = 0.0_rp
2551  enddo
2552  enddo
2553  enddo
2554  enddo
2555  endif
2556 #endif
2557 
2558  if ( tracer_type == 'SUZUKI10' ) then
2559  write(*,*) 'xxx SBM cannot be used on turbulence. Check!'
2560  call prc_mpistop
2561  endif
2562 
2563  return
2564  end subroutine mkinit_turbulence
2565 
2566  !-----------------------------------------------------------------------------
2568  subroutine mkinit_cavityflow
2569  implicit none
2570 
2571  ! Nondimenstional numbers for a cavity flow problem
2572  real(RP) :: REYNOLDS_NUM = 4.d02
2573  real(RP) :: MACH_NUM = 3.d-2
2574  real(RP) :: TEMP0 = 300.d0
2575 
2576  namelist / param_mkinit_cavityflow / &
2577  reynolds_num, &
2578  mach_num, &
2579  temp0
2580 
2581  real(RP) :: DENS0
2582  real(RP) :: TEMP
2583  real(RP) :: Gam, Cs2
2584 
2585  real(RP), parameter :: PRES0 = 1000.e+2_rp
2586  real(RP), parameter :: Ulid = 10.0_rp
2587 
2588  integer :: k, i, j
2589  integer :: ierr
2590  !---------------------------------------------------------------------------
2591 
2592  if( io_l ) write(io_fid_log,*)
2593  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit CAVITYFLOW] / Categ[preprocess] / Origin[SCALE-RM]'
2594 
2595  !--- read namelist
2596  rewind(io_fid_conf)
2597  read(io_fid_conf,nml=param_mkinit_cavityflow,iostat=ierr)
2598 
2599  if( ierr < 0 ) then !--- missing
2600  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2601  elseif( ierr > 0 ) then !--- fatal error
2602  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_CAVITYFLOW. Check!'
2603  call prc_mpistop
2604  endif
2605  if( io_lnml ) write(io_fid_log,nml=param_mkinit_cavityflow)
2606 
2607  gam = cpdry / ( cpdry - rdry )
2608  cs2 = ( ulid / mach_num )**2
2609  temp = cs2 / ( gam * rdry )
2610  dens0 = pres0 / ( rdry * temp )
2611 
2612  if( io_l ) write(io_fid_log,*) "DENS = ", dens0
2613  if( io_l ) write(io_fid_log,*) "PRES = ", pres0
2614  if( io_l ) write(io_fid_log,*) "TEMP = ", rhot(10,10,4)/dens0, temp
2615  if( io_l ) write(io_fid_log,*) "Ulid = ", ulid
2616  if( io_l ) write(io_fid_log,*) "Cs = ", sqrt(cs2)
2617 
2618  do j = 1, ja
2619  do i = 1, ia
2620  do k = ks, ke
2621  dens(k,i,j) = dens0
2622  momz(k,i,j) = 0.0_rp
2623  momx(k,i,j) = 0.0_rp
2624  momy(k,i,j) = 0.0_rp
2625  pres(k,i,j) = pres0
2626  rhot(k,i,j) = p00/rdry * (p00/pres0)**((rdry - cpdry)/cpdry)
2627  qtrc(k,i,j,:) = 0.0_rp
2628  enddo
2629  enddo
2630  enddo
2631 
2632  momx(ke+1:ka,:,:) = dens0 * ulid
2633 
2634  return
2635  end subroutine mkinit_cavityflow
2636 
2637  !-----------------------------------------------------------------------------
2639  subroutine mkinit_mountainwave
2640  implicit none
2641 
2642  ! Surface state
2643  real(RP) :: SFC_THETA ! surface potential temperature [K]
2644  real(RP) :: SFC_PRES ! surface pressure [Pa]
2645  ! Environment state
2646  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
2647  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2648 
2649  real(RP) :: SCORER = 2.e-3_rp ! Scorer parameter (~=N/U) [1/m]
2650  real(RP) :: BBL_NC = 0.0_rp ! extremum of NC in bubble [kg/kg]
2651 
2652  namelist / param_mkinit_mountainwave / &
2653  sfc_theta, &
2654  sfc_pres, &
2655  env_u, &
2656  env_v, &
2657  scorer, &
2658  bbl_nc
2659 
2660  real(RP) :: Ustar2, N2
2661 
2662  integer :: ierr
2663  integer :: k, i, j
2664  !---------------------------------------------------------------------------
2665 
2666  if( io_l ) write(io_fid_log,*)
2667  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit MOUNTAINWAVE] / Categ[preprocess] / Origin[SCALE-RM]'
2668 
2669  sfc_theta = thetastd
2670  sfc_pres = pstd
2671 
2672  !--- read namelist
2673  rewind(io_fid_conf)
2674  read(io_fid_conf,nml=param_mkinit_mountainwave,iostat=ierr)
2675 
2676  if( ierr < 0 ) then !--- missing
2677  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2678  elseif( ierr > 0 ) then !--- fatal error
2679  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_MOUNTAINWAVE. Check!'
2680  call prc_mpistop
2681  endif
2682  if( io_lnml ) write(io_fid_log,nml=param_mkinit_mountainwave)
2683 
2684  ! calc in dry condition
2685  do j = js, je
2686  do i = is, ie
2687  pres_sfc(1,i,j) = sfc_pres
2688  pott_sfc(1,i,j) = sfc_theta
2689  qv_sfc(1,i,j) = 0.0_rp
2690  qc_sfc(1,i,j) = 0.0_rp
2691  enddo
2692  enddo
2693 
2694  do j = js, je
2695  do i = is, ie
2696  do k = ks, ke
2697  ustar2 = env_u * env_u + env_v * env_v
2698  n2 = ustar2 * (scorer*scorer)
2699 
2700  pott(k,i,j) = sfc_theta * exp( n2 / grav * real_cz(k,i,j) )
2701  qv(k,i,j) = 0.0_rp
2702  qc(k,i,j) = 0.0_rp
2703  enddo
2704  enddo
2705  enddo
2706 
2707  ! make density & pressure profile in dry condition
2708  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
2709  temp(:,:,:), & ! [OUT]
2710  pres(:,:,:), & ! [OUT]
2711  pott(:,:,:), & ! [IN]
2712  qv(:,:,:), & ! [IN]
2713  qc(:,:,:), & ! [IN]
2714  temp_sfc(:,:,:), & ! [OUT]
2715  pres_sfc(:,:,:), & ! [IN]
2716  pott_sfc(:,:,:), & ! [IN]
2717  qv_sfc(:,:,:), & ! [IN]
2718  qc_sfc(:,:,:) ) ! [IN]
2719 
2720  do j = js, je
2721  do i = is, ie
2722  do k = ks, ke
2723  dens(k,i,j) = dens(k,i,j)
2724  momz(k,i,j) = 0.0_rp
2725  momx(k,i,j) = env_u * dens(k,i,j)
2726  momy(k,i,j) = env_v * dens(k,i,j)
2727  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
2728 
2729  qtrc(k,i,j,:) = 0.0_rp
2730  enddo
2731  enddo
2732  enddo
2733 
2734  ! optional : add tracer bubble
2735  if ( bbl_nc > 0.0_rp ) then
2736  if ( i_nc > 0 ) then
2737  do j = js, je
2738  do i = is, ie
2739  do k = ks, ke
2740  qtrc(k,i,j,i_nc) = bbl_nc * bubble(k,i,j)
2741  enddo
2742  enddo
2743  enddo
2744  else
2745  write(*,*) 'xxx tracer I_NC is not defined. Check!'
2746  call prc_mpistop
2747  endif
2748  endif
2749 
2750  return
2751  end subroutine mkinit_mountainwave
2752 
2753  !-----------------------------------------------------------------------------
2755  subroutine mkinit_warmbubble
2756  implicit none
2757 
2758  ! Surface state
2759  real(RP) :: SFC_THETA ! surface potential temperature [K]
2760  real(RP) :: SFC_PRES ! surface pressure [Pa]
2761  real(RP) :: SFC_RH = 80.0_rp ! surface relative humidity [%]
2762  ! Environment state
2763  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
2764  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2765  real(RP) :: ENV_RH = 80.0_rp ! Relative Humidity of environment [%]
2766  real(RP) :: ENV_L1_ZTOP = 1.e3_rp ! top height of the layer1 (constant THETA) [m]
2767  real(RP) :: ENV_L2_ZTOP = 14.e3_rp ! top height of the layer2 (small THETA gradient) [m]
2768  real(RP) :: ENV_L2_TLAPS = 4.e-3_rp ! Lapse rate of THETA in the layer2 (small THETA gradient) [K/m]
2769  real(RP) :: ENV_L3_TLAPS = 3.e-2_rp ! Lapse rate of THETA in the layer3 (large THETA gradient) [K/m]
2770  ! Bubble
2771  real(RP) :: BBL_THETA = 1.0_rp ! extremum of temperature in bubble [K]
2772 
2773  namelist / param_mkinit_warmbubble / &
2774  sfc_theta, &
2775  sfc_pres, &
2776  env_u, &
2777  env_v, &
2778  env_rh, &
2779  env_l1_ztop, &
2780  env_l2_ztop, &
2781  env_l2_tlaps, &
2782  env_l3_tlaps, &
2783  bbl_theta
2784 
2785  integer :: ierr
2786  integer :: k, i, j
2787  !---------------------------------------------------------------------------
2788 
2789  if( io_l ) write(io_fid_log,*)
2790  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit WARMBUBBLE] / Categ[preprocess] / Origin[SCALE-RM]'
2791 
2792  sfc_theta = thetastd
2793  sfc_pres = pstd
2794 
2795  !--- read namelist
2796  rewind(io_fid_conf)
2797  read(io_fid_conf,nml=param_mkinit_warmbubble,iostat=ierr)
2798 
2799  if( ierr < 0 ) then !--- missing
2800  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2801  elseif( ierr > 0 ) then !--- fatal error
2802  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_WARMBUBBLE. Check!'
2803  call prc_mpistop
2804  endif
2805  if( io_lnml ) write(io_fid_log,nml=param_mkinit_warmbubble)
2806 
2807  ! calc in dry condition
2808  pres_sfc(1,1,1) = sfc_pres
2809  pott_sfc(1,1,1) = sfc_theta
2810  qv_sfc(1,1,1) = 0.0_rp
2811  qc_sfc(1,1,1) = 0.0_rp
2812 
2813  do k = ks, ke
2814  if ( grid_cz(k) <= env_l1_ztop ) then ! Layer 1
2815  pott(k,1,1) = sfc_theta
2816  elseif( grid_cz(k) < env_l2_ztop ) then ! Layer 2
2817  pott(k,1,1) = pott(k-1,1,1) + env_l2_tlaps * ( grid_cz(k)-grid_cz(k-1) )
2818  else ! Layer 3
2819  pott(k,1,1) = pott(k-1,1,1) + env_l3_tlaps * ( grid_cz(k)-grid_cz(k-1) )
2820  endif
2821  qv(k,1,1) = 0.0_rp
2822  qc(k,1,1) = 0.0_rp
2823  enddo
2824 
2825  ! make density & pressure profile in dry condition
2826  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2827  temp(:,1,1), & ! [OUT]
2828  pres(:,1,1), & ! [OUT]
2829  pott(:,1,1), & ! [IN]
2830  qv(:,1,1), & ! [IN]
2831  qc(:,1,1), & ! [IN]
2832  temp_sfc(1,1,1), & ! [OUT]
2833  pres_sfc(1,1,1), & ! [IN]
2834  pott_sfc(1,1,1), & ! [IN]
2835  qv_sfc(1,1,1), & ! [IN]
2836  qc_sfc(1,1,1) ) ! [IN]
2837 
2838  ! calc QV from RH
2839  call saturation_pres2qsat_all( qsat_sfc(1,1,1), temp_sfc(1,1,1), pres_sfc(1,1,1) )
2840  call saturation_pres2qsat_all( qsat(:,1,1), temp(:,1,1), pres(:,1,1) )
2841 
2842  qv_sfc(1,1,1) = sfc_rh * 1.e-2_rp * qsat_sfc(1,1,1)
2843  do k = ks, ke
2844  if ( grid_cz(k) <= env_l1_ztop ) then ! Layer 1
2845  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
2846  elseif( grid_cz(k) <= env_l2_ztop ) then ! Layer 2
2847  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
2848  else ! Layer 3
2849  qv(k,1,1) = 0.0_rp
2850  endif
2851  enddo
2852 
2853  ! make density & pressure profile in moist condition
2854  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
2855  temp(:,1,1), & ! [OUT]
2856  pres(:,1,1), & ! [OUT]
2857  pott(:,1,1), & ! [IN]
2858  qv(:,1,1), & ! [IN]
2859  qc(:,1,1), & ! [IN]
2860  temp_sfc(1,1,1), & ! [OUT]
2861  pres_sfc(1,1,1), & ! [IN]
2862  pott_sfc(1,1,1), & ! [IN]
2863  qv_sfc(1,1,1), & ! [IN]
2864  qc_sfc(1,1,1) ) ! [IN]
2865 
2866  do j = js, je
2867  do i = is, ie
2868  do k = ks, ke
2869  dens(k,i,j) = dens(k,1,1)
2870  momz(k,i,j) = 0.0_rp
2871  momx(k,i,j) = env_u * dens(k,i,j)
2872  momy(k,i,j) = env_v * dens(k,i,j)
2873 
2874  ! make warm bubble
2875  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
2876 
2877  qtrc(k,i,j,i_qv) = qv(k,1,1)
2878  enddo
2879  enddo
2880  enddo
2881 
2882  call flux_setup
2883 
2884  return
2885  end subroutine mkinit_warmbubble
2886 
2887  !-----------------------------------------------------------------------------
2889  subroutine mkinit_supercell
2890  implicit none
2891 
2892  real(RP) :: RHO(ka)
2893  real(RP) :: VELX(ka)
2894  real(RP) :: VELY(ka)
2895  real(RP) :: POTT(ka)
2896  real(RP) :: QV(ka)
2897 
2898  ! Bubble
2899  real(RP) :: BBL_THETA = 3.d0 ! extremum of temperature in bubble [K]
2900 
2901  namelist / param_mkinit_supercell / &
2902  bbl_theta
2903 
2904  integer :: ierr
2905  integer :: k, i, j
2906  !---------------------------------------------------------------------------
2907 
2908  if( io_l ) write(io_fid_log,*)
2909  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit SUPERCELL] / Categ[preprocess] / Origin[SCALE-RM]'
2910 
2911  !--- read namelist
2912  rewind(io_fid_conf)
2913  read(io_fid_conf,nml=param_mkinit_supercell,iostat=ierr)
2914 
2915  if( ierr < 0 ) then !--- missing
2916  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2917  elseif( ierr > 0 ) then !--- fatal error
2918  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_SUPERCELL. Check!'
2919  call prc_mpistop
2920  endif
2921  if( io_lnml ) write(io_fid_log,nml=param_mkinit_supercell)
2922 
2923  call read_sounding( rho, velx, vely, pott, qv ) ! (out)
2924 
2925  do j = js, je
2926  do i = is, ie
2927  do k = ks, ke
2928  dens(k,i,j) = rho(k)
2929  momz(k,i,j) = 0.0_rp
2930  momx(k,i,j) = rho(k) * velx(k)
2931  momy(k,i,j) = rho(k) * vely(k)
2932 
2933  ! make warm bubble
2934  rhot(k,i,j) = rho(k) * ( pott(k) + bbl_theta * bubble(k,i,j) )
2935 
2936  qtrc(k,i,j,i_qv) = qv(k)
2937  enddo
2938  enddo
2939  enddo
2940 
2941  call flux_setup
2942 
2943  return
2944  end subroutine mkinit_supercell
2945 
2946  !-----------------------------------------------------------------------------
2948  subroutine mkinit_squallline
2949  implicit none
2950 
2951  real(RP) :: RHO(ka)
2952  real(RP) :: VELX(ka)
2953  real(RP) :: VELY(ka)
2954  real(RP) :: POTT(ka)
2955  real(RP) :: QV(ka)
2956 
2957  real(RP) :: RANDOM_THETA = 0.01_rp
2958  real(RP) :: OFFSET_velx = 12.0_rp
2959  real(RP) :: OFFSET_vely = -2.0_rp
2960 
2961  namelist / param_mkinit_squallline / &
2962  random_theta, &
2963  offset_velx, &
2964  offset_vely
2965 
2966  integer :: ierr
2967  integer :: k, i, j
2968  !---------------------------------------------------------------------------
2969 
2970  if( io_l ) write(io_fid_log,*)
2971  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit SQUALLLINE] / Categ[preprocess] / Origin[SCALE-RM]'
2972 
2973  !--- read namelist
2974  rewind(io_fid_conf)
2975  read(io_fid_conf,nml=param_mkinit_squallline,iostat=ierr)
2976 
2977  if( ierr < 0 ) then !--- missing
2978  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
2979  elseif( ierr > 0 ) then !--- fatal error
2980  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_SQUALLLINE. Check!'
2981  call prc_mpistop
2982  endif
2983  if( io_lnml ) write(io_fid_log,nml=param_mkinit_squallline)
2984 
2985  call read_sounding( rho, velx, vely, pott, qv ) ! (out)
2986 
2987  call random_get(rndm) ! make random
2988  do j = js, je
2989  do i = is, ie
2990  do k = ks, ke
2991  dens(k,i,j) = rho(k)
2992  momz(k,i,j) = 0.0_rp
2993  momx(k,i,j) = ( velx(k) - offset_velx ) * rho(k)
2994  momy(k,i,j) = ( vely(k) - offset_vely ) * rho(k)
2995  rhot(k,i,j) = rho(k) * ( pott(k) + rndm(k,i,j) * random_theta )
2996 
2997  qtrc(k,i,j,i_qv) = qv(k)
2998  enddo
2999  enddo
3000  enddo
3001 
3002  call flux_setup
3003 
3004  call ocean_setup
3005 
3006  return
3007  end subroutine mkinit_squallline
3008 
3009  !-----------------------------------------------------------------------------
3011  subroutine mkinit_wk1982
3012  implicit none
3013 
3014  ! Surface state
3015  real(RP) :: SFC_THETA = 300.0_rp ! surface pot. temperature [K]
3016  real(RP) :: SFC_PRES ! surface pressure [Pa]
3017  ! Parameter in Weisman and Klemp (1982)
3018  real(RP) :: TR_Z = 12000.0_rp ! height of tropopause [m]
3019  real(RP) :: TR_THETA = 343.0_rp ! pot. temperature at tropopause [K]
3020  real(RP) :: TR_TEMP = 213.0_rp ! temperature at tropopause [K]
3021  real(RP) :: SHEAR_Z = 3000.0_rp ! center height of shear layer [m]
3022  real(RP) :: SHEAR_U = 15.0_rp ! velocity u over the shear layer [m/s]
3023  ! Bubble
3024  real(RP) :: BBL_THETA = 3.d0 ! extremum of temperature in bubble [K]
3025 
3026  namelist / param_mkinit_wk1982 / &
3027  sfc_theta, &
3028  sfc_pres, &
3029  tr_z, &
3030  tr_theta, &
3031  tr_temp, &
3032  shear_z, &
3033  shear_u, &
3034  bbl_theta
3035 
3036  real(RP) :: rh (ka,ia,ja)
3037  real(RP) :: rh_sfc(1 ,ia,ja)
3038 
3039  integer :: ierr
3040  integer :: k, i, j, iq
3041  !---------------------------------------------------------------------------
3042 
3043  if( io_l ) write(io_fid_log,*)
3044  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit WK1982] / Categ[preprocess] / Origin[SCALE-RM]'
3045 
3046  sfc_pres = pstd
3047 
3048  rewind(io_fid_conf)
3049  read(io_fid_conf,nml=param_mkinit_wk1982,iostat=ierr)
3050  if( ierr < 0 ) then !--- missing
3051  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
3052  elseif( ierr > 0 ) then !--- fatal error
3053  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_WK1982. Check!'
3054  call prc_mpistop
3055  endif
3056  if( io_lnml ) write(io_fid_log,nml=param_mkinit_wk1982)
3057 
3058  ! calc in dry condition
3059  do j = js, je
3060  do i = is, ie
3061  pres_sfc(1,i,j) = sfc_pres
3062  pott_sfc(1,i,j) = sfc_theta
3063  qv_sfc(1,i,j) = 0.0_rp
3064  qc_sfc(1,i,j) = 0.0_rp
3065 
3066  do k = ks, ke
3067  if ( real_cz(k,i,j) <= tr_z ) then ! below initial cloud top
3068  pott(k,i,j) = pott_sfc(1,i,j) &
3069  + ( tr_theta - pott_sfc(1,i,j) ) * ( real_cz(k,i,j) / tr_z )**1.25_rp
3070  else
3071  pott(k,i,j) = tr_theta * exp( grav * ( real_cz(k,i,j) - tr_z ) / cpdry / tr_temp )
3072  endif
3073 
3074  qv(k,i,j) = 0.0_rp
3075  qc(k,i,j) = 0.0_rp
3076  enddo
3077  enddo
3078  enddo
3079 
3080  ! make density & pressure profile in dry condition
3081  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3082  temp(:,:,:), & ! [OUT]
3083  pres(:,:,:), & ! [OUT]
3084  pott(:,:,:), & ! [IN]
3085  qv(:,:,:), & ! [IN]
3086  qc(:,:,:), & ! [IN]
3087  temp_sfc(:,:,:), & ! [OUT]
3088  pres_sfc(:,:,:), & ! [IN]
3089  pott_sfc(:,:,:), & ! [IN]
3090  qv_sfc(:,:,:), & ! [IN]
3091  qc_sfc(:,:,:) ) ! [IN]
3092 
3093  ! calc QV from RH
3094  do j = js, je
3095  do i = is, ie
3096  rh_sfc(1,i,j) = 1.0_rp - 0.75_rp * ( real_fz(ks-1,i,j) / tr_z )**1.25_rp
3097 
3098  do k = ks, ke
3099  if ( real_cz(k,i,j) <= tr_z ) then ! below initial cloud top
3100  rh(k,i,j) = 1.0_rp - 0.75_rp * ( real_cz(k,i,j) / tr_z )**1.25_rp
3101  else
3102  rh(k,i,j) = 0.25_rp
3103  endif
3104  enddo
3105  enddo
3106  enddo
3107 
3108  call saturation_pres2qsat_all( qsat_sfc(1,:,:), temp_sfc(1,:,:), pres_sfc(1,:,:) )
3109  call saturation_pres2qsat_all( qsat(:,:,:), temp(:,:,:), pres(:,:,:) )
3110 
3111  do j = js, je
3112  do i = is, ie
3113  qv_sfc(1,i,j) = rh_sfc(1,i,j) * qsat_sfc(1,i,j)
3114  do k = ks, ke
3115  qv(k,i,j) = rh(k,i,j) * qsat(k,i,j)
3116  enddo
3117  enddo
3118  enddo
3119 
3120  ! make density & pressure profile in moist condition
3121  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3122  temp(:,:,:), & ! [OUT]
3123  pres(:,:,:), & ! [OUT]
3124  pott(:,:,:), & ! [IN]
3125  qv(:,:,:), & ! [IN]
3126  qc(:,:,:), & ! [IN]
3127  temp_sfc(:,:,:), & ! [OUT]
3128  pres_sfc(:,:,:), & ! [IN]
3129  pott_sfc(:,:,:), & ! [IN]
3130  qv_sfc(:,:,:), & ! [IN]
3131  qc_sfc(:,:,:) ) ! [IN]
3132 
3133 do k = ks, ke
3134  if( io_l ) write(io_fid_log,*) k, real_cz(k,is,js), pres(k,is,js), pott(k,is,js), rh(k,is,js), qv(k,is,js)*1000
3135 enddo
3136 
3137  call comm_vars8( dens(:,:,:), 1 )
3138  call comm_wait ( dens(:,:,:), 1 )
3139 
3140  do j = js, je
3141  do i = is, ie
3142  do k = ks, ke
3143  momx(k,i,j) = shear_u * tanh( real_cz(k,i,j) / shear_z ) &
3144  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3145  enddo
3146  enddo
3147  enddo
3148 
3149  do j = js, je
3150  do i = is, ie
3151  do k = ks, ke
3152  momy(k,i,j) = 0.0_rp
3153  momz(k,i,j) = 0.0_rp
3154  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3155 
3156  ! make warm bubble
3157  rhot(k,i,j) = dens(k,i,j) * ( pott(k,i,j) + bbl_theta * bubble(k,i,j) )
3158 
3159  qtrc(k,i,j,i_qv) = qv(k,i,j)
3160  enddo
3161  enddo
3162  enddo
3163 
3164 #ifndef DRY
3165  if ( qa >= 2 ) then
3166  do iq = 2, qa
3167  do j = js, je
3168  do i = is, ie
3169  do k = ks, ke
3170  qtrc(k,i,j,iq) = 0.0_rp
3171  enddo
3172  enddo
3173  enddo
3174  enddo
3175  endif
3176 #endif
3177 
3178  call flux_setup
3179 
3180  return
3181  end subroutine mkinit_wk1982
3182 
3183  !-----------------------------------------------------------------------------
3185  subroutine mkinit_dycoms2_rf01
3186  implicit none
3187 
3188 #ifndef DRY
3189  real(RP) :: PERTURB_AMP = 0.0_rp
3190  integer :: RANDOM_LIMIT = 5
3191  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
3192  ! 1 -> petrurbation for pt
3193  ! 2 -> perturbation for u, v, w
3194  logical :: USE_LWSET = .false. ! use liq. water. static energy temp.?
3195 
3196  namelist / param_mkinit_rf01 / &
3197  perturb_amp, &
3198  random_limit, &
3199  random_flag, &
3200  use_lwset
3201 
3202  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3203 
3204  real(RP) :: qall ! QV+QC
3205  real(RP) :: fact
3206  real(RP) :: pi2
3207  real(RP) :: sint
3208  real(RP) :: GEOP_sw ! switch for geopotential energy correction
3209 
3210  integer :: ierr
3211  integer :: k, i, j, iq
3212  !---------------------------------------------------------------------------
3213 
3214  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
3215 
3216  if( io_l ) write(io_fid_log,*)
3217  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit DYCOMS2RF01] / Categ[preprocess] / Origin[SCALE-RM]'
3218 
3219  rewind(io_fid_conf)
3220  read(io_fid_conf,nml=param_mkinit_rf01,iostat=ierr)
3221  if( ierr < 0 ) then !--- missing
3222  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
3223  elseif( ierr > 0 ) then !--- fatal error
3224  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_RF01. Check!'
3225  call prc_mpistop
3226  endif
3227  if( io_lnml ) write(io_fid_log,nml=param_mkinit_rf01)
3228 
3229  if ( use_lwset ) then
3230  geop_sw = 1.0_rp
3231  else
3232  geop_sw = 0.0_rp
3233  endif
3234 
3235  ! calc in dry condition
3236  do j = js, je
3237  do i = is, ie
3238 
3239  pres_sfc(1,i,j) = 1017.8e2_rp ! [Pa]
3240  pott_sfc(1,i,j) = 289.0_rp ! [K]
3241  qv_sfc(1,i,j) = 0.0_rp
3242  qc_sfc(1,i,j) = 0.0_rp
3243 
3244  do k = ks, ke
3245  velx(k,i,j) = 7.0_rp
3246  vely(k,i,j) = -5.5_rp
3247  if ( grid_cz(k) < 820.0_rp ) then ! below initial cloud top
3248  potl(k,i,j) = 289.0_rp - grav / cpdry * grid_cz(k) * geop_sw
3249  elseif( grid_cz(k) <= 860.0_rp ) then
3250  sint = sin( pi2 * ( grid_cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
3251  potl(k,i,j) = ( 289.0_rp - grav / cpdry * grid_cz(k) * geop_sw ) * (0.5_rp-sint) &
3252  + ( 297.5_rp+sign(abs(grid_cz(k)-840.0_rp)**(1.0_rp/3.0_rp),grid_cz(k)-840.0_rp) &
3253  - grav / cpdry * grid_cz(k) * geop_sw ) * (0.5_rp+sint)
3254  else
3255  potl(k,i,j) = 297.5_rp + ( grid_cz(k)-840.0_rp )**(1.0_rp/3.0_rp) &
3256  - grav / cpdry * grid_cz(k) * geop_sw
3257  endif
3258 
3259  qv(k,i,j) = 0.0_rp
3260  qc(k,i,j) = 0.0_rp
3261  enddo
3262 
3263  enddo
3264  enddo
3265 
3266  ! make density & pressure profile in dry condition
3267  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3268  temp(:,:,:), & ! [OUT]
3269  pres(:,:,:), & ! [OUT]
3270  potl(:,:,:), & ! [IN]
3271  qv(:,:,:), & ! [IN]
3272  qc(:,:,:), & ! [IN]
3273  temp_sfc(:,:,:), & ! [OUT]
3274  pres_sfc(:,:,:), & ! [IN]
3275  pott_sfc(:,:,:), & ! [IN]
3276  qv_sfc(:,:,:), & ! [IN]
3277  qc_sfc(:,:,:) ) ! [IN]
3278 
3279  ! calc in moist condition
3280  do j = js, je
3281  do i = is, ie
3282  qv_sfc(1,i,j) = 9.0e-3_rp ! [kg/kg]
3283  qc_sfc(1,i,j) = 0.0_rp
3284 
3285  do k = ks, ke
3286  if ( grid_cz(k) < 820.0_rp ) then ! below initial cloud top
3287  qall = 9.0e-3_rp
3288  elseif( grid_cz(k) <= 860.0_rp ) then ! boundary
3289  sint = sin( pi2 * ( grid_cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
3290  qall = 9.0e-3_rp * (0.5_rp-sint) &
3291  + 1.5e-3_rp * (0.5_rp+sint)
3292  elseif( grid_cz(k) <= 5000.0_rp ) then
3293  qall = 1.5e-3_rp
3294  else
3295  qall = 0.0_rp
3296  endif
3297 
3298  if ( grid_cz(k) <= 600.0_rp ) then
3299  qc(k,i,j) = 0.0_rp
3300  elseif( grid_cz(k) < 820.0_rp ) then ! in the cloud
3301  fact = ( grid_cz(k)-600.0_rp ) / ( 840.0_rp-600.0_rp )
3302  qc(k,i,j) = 0.45e-3_rp * fact
3303  elseif( grid_cz(k) <= 860.0_rp ) then ! boundary
3304  sint = sin( pi2 * ( grid_cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
3305  fact = ( grid_cz(k)-600.0_rp ) / ( 840.0_rp-600.0_rp )
3306  qc(k,i,j) = 0.45e-3_rp * fact * (0.5_rp-sint)
3307  else
3308  qc(k,i,j) = 0.0_rp
3309  endif
3310 
3311  qv(k,i,j) = qall - qc(k,i,j)
3312  enddo
3313 
3314  enddo
3315  enddo
3316 
3317  do j = js, je
3318  do i = is, ie
3319  do k = ks, ke
3320  temp(k,i,j) = temp(k,i,j) + lhv / cpdry * qc(k,i,j)
3321  enddo
3322  enddo
3323  enddo
3324 
3325  ! make density & pressure profile in moist condition
3326  call hydrostatic_buildrho_bytemp( dens(:,:,:), & ! [OUT]
3327  pott(:,:,:), & ! [OUT]
3328  pres(:,:,:), & ! [OUT]
3329  temp(:,:,:), & ! [IN]
3330  qv(:,:,:), & ! [IN]
3331  qc(:,:,:), & ! [IN]
3332  pott_sfc(:,:,:), & ! [OUT]
3333  pres_sfc(:,:,:), & ! [IN]
3334  temp_sfc(:,:,:), & ! [IN]
3335  qv_sfc(:,:,:), & ! [IN]
3336  qc_sfc(:,:,:) ) ! [IN]
3337 
3338  do j = js, je
3339  do i = is, ie
3340  dens( 1:ks-1,i,j) = dens(ks,i,j)
3341  dens(ke+1:ka, i,j) = dens(ke,i,j)
3342  enddo
3343  enddo
3344 
3345  call comm_vars8( dens(:,:,:), 1 )
3346  call comm_wait ( dens(:,:,:), 1 )
3347 
3348  call random_get(rndm) ! make random
3349  do j = js, je
3350  do i = is, ie
3351  do k = ks, ke
3352  if ( random_flag == 2 .and. k <= random_limit ) then ! below initial cloud top
3353  momz(k,i,j) = ( 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3354  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3355  else
3356  momz(k,i,j) = 0.0_rp
3357  endif
3358  enddo
3359  enddo
3360  enddo
3361 
3362  call random_get(rndm) ! make random
3363  do j = js, je
3364  do i = is, ie
3365  do k = ks, ke
3366  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
3367  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3368  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3369  else
3370  momx(k,i,j) = velx(k,i,j) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3371  endif
3372  enddo
3373  enddo
3374  enddo
3375 
3376  call random_get(rndm) ! make random
3377  do j = js, je
3378  do i = is, ie
3379  do k = ks, ke
3380  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
3381  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3382  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3383  else
3384  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3385  endif
3386  enddo
3387  enddo
3388  enddo
3389 
3390  call random_get(rndm) ! make random
3391  do j = js, je
3392  do i = is, ie
3393  do k = ks, ke
3394  if ( random_flag == 1 .and. k <= random_limit ) then ! below initial cloud top
3395  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3396  * dens(k,i,j)
3397  else
3398  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3399  endif
3400  enddo
3401  enddo
3402  enddo
3403 
3404  if ( tracer_type == 'SUZUKI10' ) then
3405  do j = js, je
3406  do i = is, ie
3407  do k = ks, ke
3408  qtrc(k,i,j,i_qv) = qv(k,i,j) + qc(k,i,j) !--- Super saturated air at initial
3409  do iq = qqa+1, qa
3410  qtrc(k,i,j,iq) = qtrc(k,i,j,iq) / dens(k,i,j)
3411  enddo
3412  enddo
3413  enddo
3414  enddo
3415  else
3416  do j = js, je
3417  do i = is, ie
3418  do k = ks, ke
3419  qtrc(k,i,j,i_qv) = qv(k,i,j)
3420  qtrc(k,i,j,i_qc) = qc(k,i,j)
3421  enddo
3422  enddo
3423  enddo
3424 
3425  if ( i_nc > 0 ) then
3426  do j = js, je
3427  do i = is, ie
3428  do k = ks, ke
3429  if ( qc(k,i,j) > 0.0_rp ) then
3430  qtrc(k,i,j,i_nc) = 120.e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3431  endif
3432  enddo
3433  enddo
3434  enddo
3435  endif
3436  endif
3437 
3438 #endif
3439  if ( aetracer_type == 'KAJINO13' ) then
3440  call aerosol_setup
3441  endif
3442  return
3443  end subroutine mkinit_dycoms2_rf01
3444 
3445  !-----------------------------------------------------------------------------
3447  subroutine mkinit_dycoms2_rf02
3448  implicit none
3449 
3450 #ifndef DRY
3451  real(RP) :: PERTURB_AMP = 0.0_rp
3452  integer :: RANDOM_LIMIT = 5
3453  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
3454  ! 1 -> perturbation for PT
3455  ! 2 -> perturbation for u,v,w
3456 
3457  namelist / param_mkinit_rf02 / &
3458  perturb_amp, &
3459  random_limit, &
3460  random_flag
3461 
3462  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3463 
3464  real(RP) :: qall ! QV+QC
3465  real(RP) :: fact
3466  real(RP) :: pi2
3467  real(RP) :: sint
3468 
3469  integer :: ierr
3470  integer :: k, i, j, iq
3471  !---------------------------------------------------------------------------
3472 
3473  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
3474  if( io_l ) write(io_fid_log,*)
3475  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit DYCOMS2RF02] / Categ[preprocess] / Origin[SCALE-RM]'
3476 
3477  rewind(io_fid_conf)
3478  read(io_fid_conf,nml=param_mkinit_rf02,iostat=ierr)
3479  if( ierr < 0 ) then !--- missing
3480  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
3481  elseif( ierr > 0 ) then !--- fatal error
3482  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_RF02. Check!'
3483  call prc_mpistop
3484  endif
3485  if( io_lnml ) write(io_fid_log,nml=param_mkinit_rf02)
3486 
3487  ! calc in dry condition
3488  call random_get(rndm) ! make random
3489  do j = js, je
3490  do i = is, ie
3491 
3492  pres_sfc(1,i,j) = 1017.8e2_rp ! [Pa]
3493  pott_sfc(1,i,j) = 288.3_rp ! [K]
3494  qv_sfc(1,i,j) = 0.0_rp
3495  qc_sfc(1,i,j) = 0.0_rp
3496 
3497  do k = ks, ke
3498  velx(k,i,j) = 3.0_rp + 4.3 * grid_cz(k)*1.e-3_rp
3499  vely(k,i,j) = -9.0_rp + 5.6 * grid_cz(k)*1.e-3_rp
3500 
3501  if ( grid_cz(k) < 775.0_rp ) then ! below initial cloud top
3502  potl(k,i,j) = 288.3_rp ! [K]
3503  else if ( grid_cz(k) <= 815.0_rp ) then
3504  sint = sin( pi2 * (grid_cz(k) - 795.0_rp)/20.0_rp )
3505  potl(k,i,j) = 288.3_rp * (1.0_rp-sint)*0.5_rp &
3506  + ( 295.0_rp+sign(abs(grid_cz(k)-795.0_rp)**(1.0_rp/3.0_rp),grid_cz(k)-795.0_rp) ) &
3507  * (1.0_rp+sint)*0.5_rp
3508  else
3509  potl(k,i,j) = 295.0_rp + ( grid_cz(k)-795.0_rp )**(1.0_rp/3.0_rp)
3510  endif
3511 
3512  qv(k,i,j) = 0.0_rp
3513  qc(k,i,j) = 0.0_rp
3514  enddo
3515  enddo
3516  enddo
3517 
3518  ! make density & pressure profile in dry condition
3519  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3520  temp(:,:,:), & ! [OUT]
3521  pres(:,:,:), & ! [OUT]
3522  potl(:,:,:), & ! [IN]
3523  qv(:,:,:), & ! [IN]
3524  qc(:,:,:), & ! [IN]
3525  temp_sfc(:,:,:), & ! [OUT]
3526  pres_sfc(:,:,:), & ! [IN]
3527  pott_sfc(:,:,:), & ! [IN]
3528  qv_sfc(:,:,:), & ! [IN]
3529  qc_sfc(:,:,:) ) ! [IN]
3530 
3531  ! calc in moist condition
3532  do j = js, je
3533  do i = is, ie
3534  qv_sfc(1,i,j) = 9.45e-3_rp
3535  qc_sfc(1,i,j) = 0.0_rp
3536 
3537  do k = ks, ke
3538  if ( grid_cz(k) < 775.0_rp ) then ! below initial cloud top
3539  qall = 9.45e-3_rp ! [kg/kg]
3540  else if ( grid_cz(k) <= 815.0_rp ) then
3541  sint = sin( pi2 * (grid_cz(k) - 795.0_rp)/20.0_rp )
3542  qall = 9.45e-3_rp * (1.0_rp-sint)*0.5_rp + &
3543  ( 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-grid_cz(k))/500.0_rp ) ) ) * (1.0_rp+sint)*0.5_rp
3544  else
3545  qall = 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-grid_cz(k))/500.0_rp ) ) ! [kg/kg]
3546  endif
3547 
3548  if( grid_cz(k) < 400.0_rp ) then
3549  qc(k,i,j) = 0.0_rp
3550  elseif( grid_cz(k) < 775.0_rp ) then
3551  fact = ( grid_cz(k)-400.0_rp ) / ( 795.0_rp-400.0_rp )
3552  qc(k,i,j) = 0.65e-3_rp * fact
3553  elseif( grid_cz(k) <= 815.0_rp ) then
3554  sint = sin( pi2 * ( grid_cz(k)-795.0_rp )/20.0_rp )
3555  fact = ( grid_cz(k)-400.0_rp ) / ( 795.0_rp-400.0_rp )
3556  qc(k,i,j) = 0.65e-3_rp * fact * (1.0_rp-sint) * 0.5_rp
3557  else
3558  qc(k,i,j) = 0.0_rp
3559  endif
3560  qv(k,i,j) = qall - qc(k,i,j)
3561  enddo
3562 
3563  enddo
3564  enddo
3565 
3566  do j = js, je
3567  do i = is, ie
3568  do k = ks, ke
3569  temp(k,i,j) = temp(k,i,j) + lhv / cpdry * qc(k,i,j)
3570  enddo
3571  enddo
3572  enddo
3573 
3574  ! make density & pressure profile in moist condition
3575  call hydrostatic_buildrho_bytemp( dens(:,:,:), & ! [OUT]
3576  pott(:,:,:), & ! [OUT]
3577  pres(:,:,:), & ! [OUT]
3578  temp(:,:,:), & ! [IN]
3579  qv(:,:,:), & ! [IN]
3580  qc(:,:,:), & ! [IN]
3581  pott_sfc(:,:,:), & ! [OUT]
3582  pres_sfc(:,:,:), & ! [IN]
3583  temp_sfc(:,:,:), & ! [IN]
3584  qv_sfc(:,:,:), & ! [IN]
3585  qc_sfc(:,:,:) ) ! [IN]
3586 
3587  do j = js, je
3588  do i = is, ie
3589  dens( 1:ks-1,i,j) = dens(ks,i,j)
3590  dens(ke+1:ka, i,j) = dens(ke,i,j)
3591  enddo
3592  enddo
3593 
3594  call comm_vars8( dens(:,:,:), 1 )
3595  call comm_wait ( dens(:,:,:), 1 )
3596 
3597  call random_get(rndm) ! make random
3598  do j = js, je
3599  do i = is, ie
3600  do k = ks, ke
3601  if( random_flag == 2 .and. k <= random_limit ) then
3602  momz(k,i,j) = ( 0.0_rp + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3603  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3604  else
3605  momz(k,i,j) = 0.0_rp
3606  endif
3607  enddo
3608  enddo
3609  enddo
3610 
3611  call random_get(rndm) ! make random
3612  do j = js, je
3613  do i = is, ie
3614  do k = ks, ke
3615  if( random_flag == 2 .and. k <= random_limit ) then
3616  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3617  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3618  else
3619  momx(k,i,j) = ( velx(k,i,j) ) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3620  endif
3621  enddo
3622  enddo
3623  enddo
3624 
3625  call random_get(rndm) ! make random
3626  do j = js, je
3627  do i = is, ie
3628  do k = ks, ke
3629  if( random_flag == 2 .and. k <= random_limit ) then
3630  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3631  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3632  else
3633  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3634  endif
3635  enddo
3636  enddo
3637  enddo
3638 
3639  call random_get(rndm) ! make random
3640  do j = js, je
3641  do i = is, ie
3642  do k = ks, ke
3643  if( random_flag == 1 .and. k <= random_limit ) then
3644  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3645  * dens(k,i,j)
3646  else
3647  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3648  endif
3649  enddo
3650  enddo
3651  enddo
3652 
3653  if ( tracer_type == 'SUZUKI10' ) then
3654  do j = js, je
3655  do i = is, ie
3656  do k = ks, ke
3657  !--- Super saturated air at initial
3658  qtrc(k,i,j,i_qv) = qv(k,i,j) + qc(k,i,j)
3659  do iq = qqa+1, qa
3660  qtrc(k,i,j,iq) = qtrc(k,i,j,iq) / dens(k,i,j)
3661  enddo
3662  enddo
3663  enddo
3664  enddo
3665  else
3666  do j = js, je
3667  do i = is, ie
3668  do k = ks, ke
3669  qtrc(k,i,j,i_qv) = qv(k,i,j)
3670  qtrc(k,i,j,i_qc) = qc(k,i,j)
3671  enddo
3672  enddo
3673  enddo
3674 
3675  if ( i_nc > 0 ) then
3676  do j = js, je
3677  do i = is, ie
3678  do k = ks, ke
3679  if ( qc(k,i,j) > 0.0_rp ) then
3680  qtrc(k,i,j,i_nc) = 55.0e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3681  endif
3682  enddo
3683  enddo
3684  enddo
3685  endif
3686  endif
3687 
3688 #endif
3689  if ( aetracer_type == 'KAJINO13' ) then
3690  call aerosol_setup
3691  endif
3692  return
3693  end subroutine mkinit_dycoms2_rf02
3694 
3695  !-----------------------------------------------------------------------------
3697  subroutine mkinit_dycoms2_rf02_dns
3698  implicit none
3699 
3700 #ifndef DRY
3701  real(RP) :: ZB = 750.0_rp ! domain bottom
3702 ! real(RP) :: ZT = 900.0_RP ! domain top
3703  real(RP) :: CONST_U = 0.0_rp
3704  real(RP) :: CONST_V = 0.0_rp
3705  real(RP) :: PRES_ZB = 93060.0_rp
3706  real(RP) :: PERTURB_AMP = 0.0_rp
3707  integer :: RANDOM_LIMIT = 5
3708  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
3709  ! 1 -> perturbation for PT
3710  ! 2 -> perturbation for u,v,w
3711 
3712  namelist / param_mkinit_rf02_dns / &
3713  zb, const_u, const_v,pres_zb,&
3714  perturb_amp, &
3715  random_limit, &
3716  random_flag
3717 
3718  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3719 
3720  real(RP) :: qall ! QV+QC
3721  real(RP) :: fact
3722  real(RP) :: pi2
3723  real(RP) :: RovCP
3724 
3725  integer :: ierr
3726  integer :: k, i, j, iq
3727  !---------------------------------------------------------------------------
3728 
3729  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
3730 
3731  if( io_l ) write(io_fid_log,*)
3732  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit DYCOMS2RF02_DNS] / Categ[preprocess] / Origin[SCALE-RM]'
3733 
3734  rewind(io_fid_conf)
3735  read(io_fid_conf,nml=param_mkinit_rf02_dns,iostat=ierr)
3736  if( ierr < 0 ) then !--- missing
3737  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
3738  elseif( ierr > 0 ) then !--- fatal error
3739  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_RF02_DNS. Check!'
3740  call prc_mpistop
3741  endif
3742  if( io_lnml ) write(io_fid_log,nml=param_mkinit_rf02_dns)
3743 
3744  ! calc in dry condition
3745  call random_get(rndm) ! make random
3746  do j = js, je
3747  do i = is, ie
3748 
3749  pres_sfc(1,i,j) = pres_zb
3750 ! pott_sfc(1,i,j) = 288.3_RP ! [K]
3751 ! qv_sfc (1,i,j) = 9.45E-3_RP
3752 ! qc_sfc (1,i,j) = 0.0_RP
3753 
3754  do k = ks, ke
3755 
3756  velx(k,i,j) = const_u
3757  vely(k,i,j) = const_v
3758 
3759 ! if ( ZB+GRID_CZ(k) < 775.0_RP ) then ! below initial cloud top
3760  if ( zb+grid_cz(k) <= 795.0_rp ) then ! below initial cloud top
3761  potl(k,i,j) = 288.3_rp ! [K]
3762  qall = 9.45e-3_rp ! [kg/kg]
3763 ! necessary?
3764 ! else if ( GRID_CZ(k) <= 815.0_RP ) then
3765 ! sint = sin( pi2 * (GRID_CZ(k) - 795.0_RP)/20.0_RP )
3766 ! potl(k,i,j) = 288.3_RP * (1.0_RP-sint)*0.5_RP + &
3767 ! ( 295.0_RP+sign(abs(GRID_CZ(k)-795.0_RP)**(1.0_RP/3.0_RP),GRID_CZ(k)-795.0_RP) ) * (1.0_RP+sint)*0.5_RP
3768 ! qall = 9.45E-3_RP * (1.0_RP-sint)*0.5_RP + &
3769 ! ( 5.E-3_RP - 3.E-3_RP * ( 1.0_RP - exp( (795.0_RP-GRID_CZ(k))/500.0_RP ) ) ) * (1.0_RP+sint)*0.5_RP
3770  else
3771  potl(k,i,j) = 295.0_rp + ( zb+grid_cz(k)-795.0_rp )**(1.0_rp/3.0_rp)
3772  qall = 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-(zb+grid_cz(k)))/500.0_rp ) ) ! [kg/kg]
3773  endif
3774 
3775  if( zb+grid_cz(k) < 400.0_rp ) then
3776  qc(k,i,j) = 0.0_rp
3777  elseif( zb+grid_cz(k) <= 795.0_rp ) then
3778  fact = ( (zb+grid_cz(k))-400.0_rp ) / ( 795.0_rp-400.0_rp )
3779  qc(k,i,j) = 0.8e-3_rp * fact
3780  else
3781  qc(k,i,j) = 0.0_rp
3782  endif
3783  qv(k,i,j) = qall - qc(k,i,j)
3784 
3785 !if(i==is.and.j==js)write(*,*)'chkk',k,cz(k)+zb,qc(k,i,j),qv(k,i,j)
3786  enddo
3787  enddo
3788  enddo
3789 
3790 !write(*,*)'chk3',ks,ke
3791  ! extrapolation (temtative)
3792  pott_sfc(1,:,:) = potl(ks,:,:)-0.5*(potl(ks+1,:,:)-potl(ks,:,:))
3793  qv_sfc(1,:,:) = qv(ks,:,:)-0.5*(qv(ks+1,:,:)-qv(ks,:,:))
3794  qc_sfc(1,:,:) = qc(ks,:,:)-0.5*(qc(ks+1,:,:)-qc(ks,:,:))
3795 
3796  ! make density & pressure profile in moist condition
3797  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3798  temp(:,:,:), & ! [OUT]
3799  pres(:,:,:), & ! [OUT]
3800  potl(:,:,:), & ! [IN]
3801  qv(:,:,:), & ! [IN]
3802  qc(:,:,:), & ! [IN]
3803  temp_sfc(:,:,:), & ! [OUT]
3804  pres_sfc(:,:,:), & ! [IN]
3805  pott_sfc(:,:,:), & ! [IN]
3806  qv_sfc(:,:,:), & ! [IN]
3807  qc_sfc(:,:,:) ) ! [IN]
3808 
3809 !write(*,*)'chk4.1'
3810  rovcp = rdry / cpdry
3811  do j = js, je
3812  do i = is, ie
3813  do k = ks, ke
3814  pott(k,i,j) = potl(k,i,j) + lhv / cpdry * qc(k,i,j) * ( p00/pres(k,i,j) )**rovcp
3815  enddo
3816  enddo
3817  enddo
3818 
3819 !write(*,*)'chk5'
3820  ! make density & pressure profile in moist condition
3821  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
3822  temp(:,:,:), & ! [OUT]
3823  pres(:,:,:), & ! [OUT]
3824  pott(:,:,:), & ! [IN]
3825  qv(:,:,:), & ! [IN]
3826  qc(:,:,:), & ! [IN]
3827  temp_sfc(:,:,:), & ! [OUT]
3828  pres_sfc(:,:,:), & ! [IN]
3829  pott_sfc(:,:,:), & ! [IN]
3830  qv_sfc(:,:,:), & ! [IN]
3831  qc_sfc(:,:,:) ) ! [IN]
3832 
3833  do j = js, je
3834  do i = is, ie
3835  dens( 1:ks-1,i,j) = dens(ks,i,j)
3836  dens(ke+1:ka, i,j) = dens(ke,i,j)
3837  enddo
3838  enddo
3839 
3840  call comm_vars8( dens(:,:,:), 1 )
3841  call comm_wait ( dens(:,:,:), 1 )
3842 
3843 !write(*,*)'chk7'
3844  call random_get(rndm) ! make random
3845  do j = js, je
3846  do i = is, ie
3847  do k = ks, ke
3848  if( random_flag == 2 .and. k <= random_limit ) then
3849  momz(k,i,j) = ( 0.0_rp + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3850  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3851  else
3852  momz(k,i,j) = 0.0_rp
3853  endif
3854  enddo
3855  enddo
3856  enddo
3857 
3858 !write(*,*)'chk8'
3859  call random_get(rndm) ! make random
3860  do j = js, je
3861  do i = is, ie
3862  do k = ks, ke
3863  if( random_flag == 2 .and. k <= random_limit ) then
3864  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3865  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3866  else
3867  momx(k,i,j) = ( velx(k,i,j) ) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3868  endif
3869  enddo
3870  enddo
3871  enddo
3872 !write(*,*)'chk9'
3873 
3874  call random_get(rndm) ! make random
3875  do j = js, je
3876  do i = is, ie
3877  do k = ks, ke
3878  if( random_flag == 2 .and. k <= random_limit ) then
3879  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3880  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3881  else
3882  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3883  endif
3884  enddo
3885  enddo
3886  enddo
3887 !write(*,*)'chk10'
3888 
3889  call random_get(rndm) ! make random
3890  do j = js, je
3891  do i = is, ie
3892  do k = ks, ke
3893  if( random_flag == 1 .and. k <= random_limit ) then
3894  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3895  * dens(k,i,j)
3896  else
3897  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3898  endif
3899  enddo
3900  enddo
3901  enddo
3902 
3903 !write(*,*)'chk11'
3904  do iq = 1, qa
3905  do j = js, je
3906  do i = is, ie
3907  do k = ks, ke
3908  qtrc(k,i,j,iq) = 0.0_rp
3909  enddo
3910  enddo
3911  enddo
3912  enddo
3913 
3914 !write(*,*)'chk12'
3915  if ( tracer_type == 'SUZUKI10' ) then
3916  do j = js, je
3917  do i = is, ie
3918  do k = ks, ke
3919  !--- Super saturated air at initial
3920  qtrc(k,i,j,i_qv) = qv(k,i,j) + qc(k,i,j)
3921 
3922  !--- for aerosol
3923  do iq = qqa+1, qa
3924  qtrc(k,i,j,iq) = gan(iq-qqa) / dens(k,i,j)
3925  enddo
3926  enddo
3927  enddo
3928  enddo
3929  else
3930  do j = js, je
3931  do i = is, ie
3932  do k = ks, ke
3933  qtrc(k,i,j,i_qv) = qv(k,i,j)
3934  qtrc(k,i,j,i_qc) = qc(k,i,j)
3935  enddo
3936  enddo
3937  enddo
3938 
3939  if ( i_nc > 0 ) then
3940  do j = js, je
3941  do i = is, ie
3942  do k = ks, ke
3943  if ( qc(k,i,j) > 0.0_rp ) then
3944  qtrc(k,i,j,i_nc) = 55.0e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3945  endif
3946  enddo
3947  enddo
3948  enddo
3949  endif
3950  endif
3951 
3952 #endif
3953  return
3954  end subroutine mkinit_dycoms2_rf02_dns
3955 
3956  !-----------------------------------------------------------------------------
3958  subroutine mkinit_rico
3959  implicit none
3960 
3961 #ifndef DRY
3962  real(RP):: PERTURB_AMP_PT = 0.1_rp
3963  real(RP):: PERTURB_AMP_QV = 2.5e-5_rp
3964 
3965  namelist / param_mkinit_rico / &
3966  perturb_amp_pt, &
3967  perturb_amp_qv
3968 
3969  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3970  real(RP) :: qall ! QV+QC
3971  real(RP) :: fact
3972 
3973  integer :: ierr
3974  integer :: k, i, j, iq
3975  !---------------------------------------------------------------------------
3976 
3977  if( io_l ) write(io_fid_log,*)
3978  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit RICO] / Categ[preprocess] / Origin[SCALE-RM]'
3979 
3980  rewind(io_fid_conf)
3981  read(io_fid_conf,nml=param_mkinit_rico,iostat=ierr)
3982  if( ierr < 0 ) then !--- missing
3983  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
3984  elseif( ierr > 0 ) then !--- fatal error
3985  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_RICO. Check!'
3986  call prc_mpistop
3987  endif
3988  if( io_lnml ) write(io_fid_log,nml=param_mkinit_rico)
3989 
3990  ! calc in moist condition
3991  do j = js, je
3992  do i = is, ie
3993 
3994  pres_sfc(1,i,j) = 1015.4e2_rp ! [Pa]
3995  pott_sfc(1,i,j) = 297.9_rp
3996  qv_sfc(1,i,j) = 0.0_rp
3997  qc_sfc(1,i,j) = 0.0_rp
3998 
3999  do k = ks, ke
4000  !--- potential temperature
4001  if ( grid_cz(k) < 740.0_rp ) then ! below initial cloud top
4002  potl(k,i,j) = 297.9_rp
4003  else
4004  fact = ( grid_cz(k)-740.0_rp ) * ( 317.0_rp-297.9_rp ) / ( 4000.0_rp-740.0_rp )
4005  potl(k,i,j) = 297.9_rp + fact
4006  endif
4007 
4008  !--- horizontal wind velocity
4009  if ( grid_cz(k) <= 4000.0_rp ) then ! below initial cloud top
4010  fact = ( grid_cz(k)-0.0_rp ) * ( -1.9_rp+9.9_rp ) / ( 4000.0_rp-0.0_rp )
4011  velx(k,i,j) = -9.9_rp + fact
4012  vely(k,i,j) = -3.8_rp
4013  else
4014  velx(k,i,j) = -1.9_rp
4015  vely(k,i,j) = -3.8_rp
4016  endif
4017 
4018  qv(k,i,j) = 0.0_rp
4019  qc(k,i,j) = 0.0_rp
4020 
4021  enddo
4022 
4023  enddo
4024  enddo
4025 
4026  ! make density & pressure profile in moist condition
4027  call hydrostatic_buildrho( dens(:,:,:), & ! [OUT]
4028  temp(:,:,:), & ! [OUT]
4029  pres(:,:,:), & ! [OUT]
4030  potl(:,:,:), & ! [IN]
4031  qv(:,:,:), & ! [IN]
4032  qc(:,:,:), & ! [IN]
4033  temp_sfc(:,:,:), & ! [OUT]
4034  pres_sfc(:,:,:), & ! [IN]
4035  pott_sfc(:,:,:), & ! [IN]
4036  qv_sfc(:,:,:), & ! [IN]
4037  qc_sfc(:,:,:) ) ! [IN]
4038 
4039 
4040  do j = js, je
4041  do i = is, ie
4042  qv_sfc(1,i,j) = 16.0e-3_rp ! [kg/kg]
4043  qc_sfc(1,i,j) = 0.0_rp
4044 
4045  do k = ks, ke
4046  !--- mixing ratio of vapor
4047  if ( grid_cz(k) <= 740.0_rp ) then ! below initial cloud top
4048  fact = ( grid_cz(k)-0.0_rp ) * ( 13.8e-3_rp-16.0e-3_rp ) / ( 740.0_rp-0.0_rp )
4049  qall = 16.0e-3_rp + fact
4050  elseif ( grid_cz(k) <= 3260.0_rp ) then ! boundary
4051  fact = ( grid_cz(k)-740.0_rp ) * ( 2.4e-3_rp-13.8e-3_rp ) / ( 3260.0_rp-740.0_rp )
4052  qall = 13.8e-3_rp + fact
4053  elseif( grid_cz(k) <= 4000.0_rp ) then
4054  fact = ( grid_cz(k)-3260.0_rp ) * ( 1.8e-3_rp-2.4e-3_rp ) / ( 4000.0_rp-3260.0_rp )
4055  qall = 2.4e-3_rp + fact
4056  else
4057  qall = 0.0_rp
4058  endif
4059 
4060  qc(k,i,j) = 0.0_rp
4061  qv(k,i,j) = qall - qc(k,i,j)
4062  enddo
4063 
4064  enddo
4065  enddo
4066 
4067  do j = js, je
4068  do i = is, ie
4069  do k = ks, ke
4070  temp(k,i,j) = temp(k,i,j) + lhv / cpdry * qc(k,i,j)
4071  enddo
4072  enddo
4073  enddo
4074 
4075  ! make density & pressure profile in moist condition
4076  call hydrostatic_buildrho_bytemp( dens(:,:,:), & ! [OUT]
4077  pott(:,:,:), & ! [OUT]
4078  pres(:,:,:), & ! [OUT]
4079  temp(:,:,:), & ! [IN]
4080  qv(:,:,:), & ! [IN]
4081  qc(:,:,:), & ! [IN]
4082  pott_sfc(:,:,:), & ! [OUT]
4083  pres_sfc(:,:,:), & ! [IN]
4084  temp_sfc(:,:,:), & ! [IN]
4085  qv_sfc(:,:,:), & ! [IN]
4086  qc_sfc(:,:,:) ) ! [IN]
4087 
4088 
4089  do j = js, je
4090  do i = is, ie
4091  dens( 1:ks-1,i,j) = dens(ks,i,j)
4092  dens(ke+1:ka ,i,j) = dens(ke,i,j)
4093  enddo
4094  enddo
4095 
4096  call comm_vars8( dens(:,:,:), 1 )
4097  call comm_wait ( dens(:,:,:), 1 )
4098 
4099  do j = js, je
4100  do i = is, ie
4101  do k = ks, ke
4102  momz(k,i,j) = 0.0_rp
4103  enddo
4104  enddo
4105  enddo
4106 
4107  do j = js, je
4108  do i = is, ie
4109  do k = ks, ke
4110  momx(k,i,j) = velx(k,i,j) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
4111  enddo
4112  enddo
4113  enddo
4114 
4115  do j = js, je
4116  do i = is, ie
4117  do k = ks, ke
4118  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
4119  enddo
4120  enddo
4121  enddo
4122 
4123  call random_get(rndm) ! make random
4124  do j = js, je
4125  do i = is, ie
4126  do k = ks, ke
4127  rhot(k,i,j) = ( pott(k,i,j)+2.0_rp*( rndm(k,i,j)-0.5_rp )*perturb_amp_pt ) * dens(k,i,j)
4128  enddo
4129  enddo
4130  enddo
4131 
4132  call random_get(rndm) ! make random
4133 
4134  if ( tracer_type == 'SUZUKI10' ) then
4135  do j = js, je
4136  do i = is, ie
4137  do k = ks, ke
4138  !--- Super saturated air at initial
4139  qtrc(k,i,j,i_qv) = qv(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp_qv &
4140  + qc(k,i,j)
4141  do iq = qqa+1, qa
4142  qtrc(k,i,j,iq) = qtrc(k,i,j,iq) / dens(k,i,j)
4143  enddo
4144  enddo
4145  enddo
4146  enddo
4147  else
4148  do j = js, je
4149  do i = is, ie
4150  do k = ks, ke
4151  qtrc(k,i,j,i_qv) = qv(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp_qv
4152  qtrc(k,i,j,i_qc) = qc(k,i,j)
4153  enddo
4154  enddo
4155  enddo
4156 
4157  if ( i_nc > 0 ) then
4158  do j = js, je
4159  do i = is, ie
4160  do k = ks, ke
4161  if ( qc(k,i,j) > 0.0_rp ) then
4162  qtrc(k,i,j,i_nc) = 70.e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
4163  endif
4164  enddo
4165  enddo
4166  enddo
4167  endif
4168  endif
4169 
4170 #endif
4171  if ( aetracer_type == 'KAJINO13' ) then
4172  call aerosol_setup
4173  endif
4174  return
4175  end subroutine mkinit_rico
4176 
4177  !-----------------------------------------------------------------------------
4178  subroutine mkinit_interporation
4179  use gtool_file, only: &
4180  filegetshape, &
4181  fileread
4182  implicit none
4183 
4184  real(RP) :: dz(ka,ia,ja)
4185 
4186  real(RP) :: W(ka,ia,ja)
4187  real(RP) :: U(ka,ia,ja)
4188  real(RP) :: V(ka,ia,ja)
4189 
4190  real(RP) :: fact_cz0(ka)
4191  real(RP) :: fact_cz1(ka)
4192  real(RP) :: fact_fz0(ka)
4193  real(RP) :: fact_fz1(ka)
4194  real(RP) :: fact_cx0(ia)
4195  real(RP) :: fact_cx1(ia)
4196  real(RP) :: fact_fx0(ia)
4197  real(RP) :: fact_fx1(ia)
4198  real(RP) :: fact_cy0(ja)
4199  real(RP) :: fact_cy1(ja)
4200  real(RP) :: fact_fy0(ja)
4201  real(RP) :: fact_fy1(ja)
4202 
4203  integer :: idx_cz0(ka)
4204  integer :: idx_cz1(ka)
4205  integer :: idx_fz0(ka)
4206  integer :: idx_fz1(ka)
4207  integer :: idx_cx0(ia)
4208  integer :: idx_cx1(ia)
4209  integer :: idx_fx0(ia)
4210  integer :: idx_fx1(ia)
4211  integer :: idx_cy0(ja)
4212  integer :: idx_cy1(ja)
4213  integer :: idx_fy0(ja)
4214  integer :: idx_fy1(ja)
4215 
4216  real(RP), allocatable :: DENS_ORG(:,:,:)
4217  real(RP), allocatable :: MOMZ_ORG(:,:,:)
4218  real(RP), allocatable :: MOMX_ORG(:,:,:)
4219  real(RP), allocatable :: MOMY_ORG(:,:,:)
4220  real(RP), allocatable :: RHOT_ORG(:,:,:)
4221  real(RP), allocatable :: QTRC_ORG(:,:,:,:)
4222 
4223  real(RP), allocatable :: W_ORG(:,:,:)
4224  real(RP), allocatable :: U_ORG(:,:,:)
4225  real(RP), allocatable :: V_ORG(:,:,:)
4226  real(RP), allocatable :: POTT_ORG(:,:,:)
4227 
4228  real(RP), allocatable :: CZ_ORG(:)
4229  real(RP), allocatable :: FZ_ORG(:)
4230  real(RP), allocatable :: CX_ORG(:)
4231  real(RP), allocatable :: FX_ORG(:)
4232  real(RP), allocatable :: CY_ORG(:)
4233  real(RP), allocatable :: FY_ORG(:)
4234 
4235  integer :: dims(3)
4236 
4237  character(len=H_LONG) :: BASENAME_ORG = ''
4238 
4239  namelist / param_mkinit_interporation / &
4240  basename_org
4241 
4242  integer :: ierr
4243  integer :: k, i, j, iq
4244  !---------------------------------------------------------------------------
4245 
4246  if( io_l ) write(io_fid_log,*)
4247  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit INTERPORATION] / Categ[preprocess] / Origin[SCALE-RM]'
4248 
4249  !--- read namelist
4250  rewind(io_fid_conf)
4251  read(io_fid_conf,nml=param_mkinit_interporation,iostat=ierr)
4252 
4253  if( ierr < 0 ) then !--- missing
4254  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
4255  elseif( ierr > 0 ) then !--- fatal error
4256  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_INTERPORATION. Check!'
4257  call prc_mpistop
4258  endif
4259  if( io_lnml ) write(io_fid_log,nml=param_mkinit_interporation)
4260 
4261  call filegetshape( dims(:), &
4262  basename_org, "DENS", 1, single=.true. )
4263 
4264  allocate( dens_org(dims(1),dims(2),dims(3)) )
4265  allocate( momz_org(dims(1),dims(2),dims(3)) )
4266  allocate( momx_org(dims(1),dims(2),dims(3)) )
4267  allocate( momy_org(dims(1),dims(2),dims(3)) )
4268  allocate( rhot_org(dims(1),dims(2),dims(3)) )
4269  allocate( qtrc_org(dims(1),dims(2),dims(3),qa) )
4270 
4271  allocate( w_org(dims(1),dims(2),dims(3)) )
4272  allocate( u_org(dims(1),dims(2),dims(3)) )
4273  allocate( v_org(dims(1),dims(2),dims(3)) )
4274  allocate( pott_org(dims(1),dims(2),dims(3)) )
4275 
4276  allocate( cz_org(dims(1)) )
4277  allocate( fz_org(dims(1)) )
4278  allocate( cx_org(dims(2)) )
4279  allocate( fx_org(dims(2)) )
4280  allocate( cy_org(dims(3)) )
4281  allocate( fy_org(dims(3)) )
4282 
4283  call fileread( dens_org(:,:,:), &
4284  basename_org, "DENS", 1, 1, single=.true. )
4285  call fileread( momz_org(:,:,:), &
4286  basename_org, "MOMZ", 1, 1, single=.true. )
4287  call fileread( momx_org(:,:,:), &
4288  basename_org, "MOMX", 1, 1, single=.true. )
4289  call fileread( momy_org(:,:,:), &
4290  basename_org, "MOMY", 1, 1, single=.true. )
4291  call fileread( rhot_org(:,:,:), &
4292  basename_org, "RHOT", 1, 1, single=.true. )
4293  do iq = 1, qa
4294  call fileread( qtrc_org(:,:,:,iq), &
4295  basename_org, aq_name(iq), 1, 1, single=.true. )
4296  end do
4297 
4298  call fileread( cz_org(:), &
4299  basename_org, "z" , 1, 1, single=.true. )
4300  call fileread( cx_org(:), &
4301  basename_org, "x" , 1, 1, single=.true. )
4302  call fileread( cy_org(:), &
4303  basename_org, "y" , 1, 1, single=.true. )
4304  call fileread( fx_org(:), &
4305  basename_org, "xh", 1, 1, single=.true. )
4306  call fileread( fy_org(:), &
4307  basename_org, "yh", 1, 1, single=.true. )
4308 
4309  do k = ks, ke
4310  call interporation_fact( fact_cz0(k), fact_cz1(k), & ! (OUT)
4311  idx_cz0( k), idx_cz1(k), & ! (OUT)
4312  grid_cz(k), cz_org, dims(1), & ! (IN)
4313  .false. ) ! (IN)
4314  call interporation_fact( fact_fz0(k), fact_fz1(k), & ! (OUT)
4315  idx_fz0(k), idx_fz1(k), & ! (OUT)
4316  grid_fz(k), fz_org, dims(1), & ! (IN)
4317  .false. ) ! (IN)
4318  enddo
4319  do i = is, ie
4320  call interporation_fact( fact_cx0(i), fact_cx1(i), & ! (OUT)
4321  idx_cx0(i), idx_cx1(i), & ! (OUT)
4322  grid_cx(i), cx_org, dims(2), & ! (IN)
4323  .true. ) ! (IN)
4324  call interporation_fact( fact_fx0(i), fact_fx1(i), & ! (OUT)
4325  idx_fx0(i), idx_fx1(i), & ! (OUT)
4326  grid_fx(i), fx_org, dims(2), & ! (IN)
4327  .true. ) ! (IN)
4328  enddo
4329  do j = js, je
4330  call interporation_fact( fact_cy0(j), fact_cy1(j), & ! (OUT)
4331  idx_cy0(j), idx_cy1(j), & ! (OUT)
4332  grid_cy(j), cy_org, dims(3), & ! (IN)
4333  .true. ) ! (IN)
4334  call interporation_fact( fact_fy0(j), fact_fy1(j), & ! (OUT)
4335  idx_fy0(j), idx_fy1(j), & ! (OUT)
4336  grid_fy(j), fy_org, dims(3), & ! (IN)
4337  .true. ) ! (IN)
4338  enddo
4339 
4340 
4341  do j = 1, dims(3)
4342  do i = 1, dims(2)
4343  do k = 1, dims(1)-1
4344  w_org(k,i,j) = 2.0_rp * momz_org(k,i,j) / ( dens_org(k+1,i,j) + dens_org(k,i,j) )
4345  end do
4346  end do
4347  end do
4348  do j = 1, dims(3)
4349  do i = 1, dims(2)
4350  w_org(dims(1),i,j) = 0.0_rp
4351  end do
4352  end do
4353 
4354  do j = 1, dims(3)
4355  do i = 1, dims(2)-1
4356  do k = 1, dims(1)
4357  u_org(k,i,j) = 2.0_rp * momx_org(k,i,j) / ( dens_org(k,i+1,j) + dens_org(k,i,j) )
4358  end do
4359  end do
4360  end do
4361  do j = 1, dims(3)
4362  do k = 1, dims(1)
4363  u_org(k,dims(2),j) = 2.0_rp * momx_org(k,dims(2),j) / ( dens_org(k,1,j) + dens_org(k,dims(2),j) )
4364  end do
4365  end do
4366 
4367  do j = 1, dims(3)-1
4368  do i = 1, dims(2)
4369  do k = 1, dims(1)
4370  v_org(k,i,j) = 2.0_rp * momy_org(k,i,j) / ( dens_org(k,i,j+1) + dens_org(k,i,j) )
4371  end do
4372  end do
4373  end do
4374  do i = 1, dims(2)
4375  do k = 1, dims(1)
4376  v_org(k,i,dims(3)) = 2.0_rp * momy_org(k,i,dims(3)) / ( dens_org(k,i,1) + dens_org(k,i,dims(3)) )
4377  end do
4378  end do
4379 
4380  do j = 1, dims(3)
4381  do i = 1, dims(2)
4382  do k = 1, dims(1)
4383  pott_org(k,i,j) = rhot_org(k,i,j) / dens_org(k,i,j)
4384  end do
4385  end do
4386  end do
4387 
4388  do j = 1, dims(3)
4389  do i = 1, dims(2)
4390  do k = 1, dims(1)
4391  dens_org(k,i,j) = log( dens_org(k,i,j) )
4392  end do
4393  end do
4394  end do
4395 
4396  do j = js, je
4397  do i = is, ie
4398  do k = ks, ie
4399  dens(k,i,j) = exp( &
4400  fact_cz0(k)*fact_cx0(i)*fact_cy0(j)*dens_org(idx_cz0(k),idx_cx0(i),idx_cy0(j)) &
4401  + fact_cz1(k)*fact_cx0(i)*fact_cy0(j)*dens_org(idx_cz1(k),idx_cx0(i),idx_cy0(j)) &
4402  + fact_cz0(k)*fact_cx1(i)*fact_cy0(j)*dens_org(idx_cz0(k),idx_cx1(i),idx_cy0(j)) &
4403  + fact_cz1(k)*fact_cx1(i)*fact_cy0(j)*dens_org(idx_cz1(k),idx_cx1(i),idx_cy0(j)) &
4404  + fact_cz0(k)*fact_cx0(i)*fact_cy1(j)*dens_org(idx_cz0(k),idx_cx0(i),idx_cy1(j)) &
4405  + fact_cz1(k)*fact_cx0(i)*fact_cy1(j)*dens_org(idx_cz1(k),idx_cx0(i),idx_cy1(j)) &
4406  + fact_cz0(k)*fact_cx1(i)*fact_cy1(j)*dens_org(idx_cz0(k),idx_cx1(i),idx_cy1(j)) &
4407  + fact_cz1(k)*fact_cx1(i)*fact_cy1(j)*dens_org(idx_cz1(k),idx_cx1(i),idx_cy1(j)) &
4408  )
4409 
4410  w(k,i,j) = fact_fz0(k)*fact_cx0(i)*fact_cy0(j)*w_org(idx_fz0(k),idx_cx0(i),idx_cy0(j)) &
4411  + fact_fz1(k)*fact_cx0(i)*fact_cy0(j)*w_org(idx_fz1(k),idx_cx0(i),idx_cy0(j)) &
4412  + fact_fz0(k)*fact_cx1(i)*fact_cy0(j)*w_org(idx_fz0(k),idx_cx1(i),idx_cy0(j)) &
4413  + fact_fz1(k)*fact_cx1(i)*fact_cy0(j)*w_org(idx_fz1(k),idx_cx1(i),idx_cy0(j)) &
4414  + fact_fz0(k)*fact_cx0(i)*fact_cy1(j)*w_org(idx_fz0(k),idx_cx0(i),idx_cy1(j)) &
4415  + fact_fz1(k)*fact_cx0(i)*fact_cy1(j)*w_org(idx_fz1(k),idx_cx0(i),idx_cy1(j)) &
4416  + fact_fz0(k)*fact_cx1(i)*fact_cy1(j)*w_org(idx_fz0(k),idx_cx1(i),idx_cy1(j)) &
4417  + fact_fz1(k)*fact_cx1(i)*fact_cy1(j)*w_org(idx_fz1(k),idx_cx1(i),idx_cy1(j))
4418 
4419  u(k,i,j) = fact_cz0(k)*fact_fx0(i)*fact_cy0(j)*u_org(idx_cz0(k),idx_fx0(i),idx_cy0(j)) &
4420  + fact_cz1(k)*fact_fx0(i)*fact_cy0(j)*u_org(idx_cz1(k),idx_fx0(i),idx_cy0(j)) &
4421  + fact_cz0(k)*fact_fx1(i)*fact_cy0(j)*u_org(idx_cz0(k),idx_fx1(i),idx_cy0(j)) &
4422  + fact_cz1(k)*fact_fx1(i)*fact_cy0(j)*u_org(idx_cz1(k),idx_fx1(i),idx_cy0(j)) &
4423  + fact_cz0(k)*fact_fx0(i)*fact_cy1(j)*u_org(idx_cz0(k),idx_fx0(i),idx_cy1(j)) &
4424  + fact_cz1(k)*fact_fx0(i)*fact_cy1(j)*u_org(idx_cz1(k),idx_fx0(i),idx_cy1(j)) &
4425  + fact_cz0(k)*fact_fx1(i)*fact_cy1(j)*u_org(idx_cz0(k),idx_fx1(i),idx_cy1(j)) &
4426  + fact_cz1(k)*fact_fx1(i)*fact_cy1(j)*u_org(idx_cz1(k),idx_fx1(i),idx_cy1(j))
4427 
4428  v(k,i,j) = fact_cz0(k)*fact_cx0(i)*fact_fy0(j)*v_org(idx_cz0(k),idx_cx0(i),idx_fy0(j)) &
4429  + fact_cz1(k)*fact_cx0(i)*fact_fy0(j)*v_org(idx_cz1(k),idx_cx0(i),idx_fy0(j)) &
4430  + fact_cz0(k)*fact_cx1(i)*fact_fy0(j)*v_org(idx_cz0(k),idx_cx1(i),idx_fy0(j)) &
4431  + fact_cz1(k)*fact_cx1(i)*fact_fy0(j)*v_org(idx_cz1(k),idx_cx1(i),idx_fy0(j)) &
4432  + fact_cz0(k)*fact_cx0(i)*fact_fy1(j)*v_org(idx_cz0(k),idx_cx0(i),idx_fy1(j)) &
4433  + fact_cz1(k)*fact_cx0(i)*fact_fy1(j)*v_org(idx_cz1(k),idx_cx0(i),idx_fy1(j)) &
4434  + fact_cz0(k)*fact_cx1(i)*fact_fy1(j)*v_org(idx_cz0(k),idx_cx1(i),idx_fy1(j)) &
4435  + fact_cz1(k)*fact_cx1(i)*fact_fy1(j)*v_org(idx_cz1(k),idx_cx1(i),idx_fy1(j))
4436 
4437  pott(k,i,j) = fact_cz0(k)*fact_cx0(i)*fact_cy0(j)*pott_org(idx_cz0(k),idx_cx0(i),idx_cy0(j)) &
4438  + fact_cz1(k)*fact_cx0(i)*fact_cy0(j)*pott_org(idx_cz1(k),idx_cx0(i),idx_cy0(j)) &
4439  + fact_cz0(k)*fact_cx1(i)*fact_cy0(j)*pott_org(idx_cz0(k),idx_cx1(i),idx_cy0(j)) &
4440  + fact_cz1(k)*fact_cx1(i)*fact_cy0(j)*pott_org(idx_cz1(k),idx_cx1(i),idx_cy0(j)) &
4441  + fact_cz0(k)*fact_cx0(i)*fact_cy1(j)*pott_org(idx_cz0(k),idx_cx0(i),idx_cy1(j)) &
4442  + fact_cz1(k)*fact_cx0(i)*fact_cy1(j)*pott_org(idx_cz1(k),idx_cx0(i),idx_cy1(j)) &
4443  + fact_cz0(k)*fact_cx1(i)*fact_cy1(j)*pott_org(idx_cz0(k),idx_cx1(i),idx_cy1(j)) &
4444  + fact_cz1(k)*fact_cx1(i)*fact_cy1(j)*pott_org(idx_cz1(k),idx_cx1(i),idx_cy1(j))
4445 
4446  do iq = 1, qa
4447  qtrc(k,i,j,iq) = fact_cz0(k)*fact_cx0(i)*fact_cy0(j)*qtrc_org(idx_cz0(k),idx_cx0(i),idx_cy0(j),iq) &
4448  + fact_cz1(k)*fact_cx0(i)*fact_cy0(j)*qtrc_org(idx_cz1(k),idx_cx0(i),idx_cy0(j),iq) &
4449  + fact_cz0(k)*fact_cx1(i)*fact_cy0(j)*qtrc_org(idx_cz0(k),idx_cx1(i),idx_cy0(j),iq) &
4450  + fact_cz1(k)*fact_cx1(i)*fact_cy0(j)*qtrc_org(idx_cz1(k),idx_cx1(i),idx_cy0(j),iq) &
4451  + fact_cz0(k)*fact_cx0(i)*fact_cy1(j)*qtrc_org(idx_cz0(k),idx_cx0(i),idx_cy1(j),iq) &
4452  + fact_cz1(k)*fact_cx0(i)*fact_cy1(j)*qtrc_org(idx_cz1(k),idx_cx0(i),idx_cy1(j),iq) &
4453  + fact_cz0(k)*fact_cx1(i)*fact_cy1(j)*qtrc_org(idx_cz0(k),idx_cx1(i),idx_cy1(j),iq) &
4454  + fact_cz1(k)*fact_cx1(i)*fact_cy1(j)*qtrc_org(idx_cz1(k),idx_cx1(i),idx_cy1(j),iq)
4455  enddo
4456  enddo
4457  enddo
4458  enddo
4459 
4460  deallocate( dens_org )
4461  deallocate( momz_org )
4462  deallocate( momx_org )
4463  deallocate( momy_org )
4464  deallocate( rhot_org )
4465  deallocate( qtrc_org )
4466 
4467  deallocate( w_org )
4468  deallocate( u_org )
4469  deallocate( v_org )
4470  deallocate( pott_org )
4471 
4472  deallocate( cz_org )
4473  deallocate( fz_org )
4474  deallocate( cx_org )
4475  deallocate( fx_org )
4476  deallocate( cy_org )
4477  deallocate( fy_org )
4478 
4479  if ( i_qv > 0 ) then
4480  do j = js, je
4481  do i = is, ie
4482  do k = ks, ke
4483  qv(k,i,j) = qtrc(k,i,j,i_qv)
4484  end do
4485  end do
4486  end do
4487  else
4488  do j = js, je
4489  do i = is, ie
4490  do k = ks, ke
4491  qv(k,i,j) = 0.0_rp
4492  end do
4493  end do
4494  end do
4495  end if
4496 
4497  if ( i_qc > 0 ) then
4498  do j = js, je
4499  do i = is, ie
4500  do k = ks, ke
4501  qc(k,i,j) = qtrc(k,i,j,i_qc)
4502  end do
4503  end do
4504  end do
4505  else
4506  do j = js, je
4507  do i = is, ie
4508  do k = ks, ke
4509  qc(k,i,j) = 0.0_rp
4510  end do
4511  end do
4512  end do
4513  end if
4514 
4515  do j = js, je
4516  do i = is, ie
4517  do k = ks+1, ke
4518  dz(k,i,j) = real_cz(k,i,j) - real_cz(k-1,i,j) ! distance from cell center to cell center
4519  enddo
4520  enddo
4521  enddo
4522 
4523  ! make density & pressure profile in moist condition
4524  call hydrostatic_buildrho_atmos( dens(:,:,:), & ! [INOUT]
4525  temp(:,:,:), & ! [OUT]
4526  pres(:,:,:), & ! [OUT]
4527  pott(:,:,:), & ! [IN]
4528  qv(:,:,:), & ! [IN]
4529  qc(:,:,:), & ! [IN]
4530  dz(:,:,:) ) ! [IN]
4531 
4532  do j = js, je
4533  do i = is, ie
4534  do k = ks, ke
4535  momz(k,i,j) = 0.5_rp * w(k,i,j) * ( dens(k,i,j) + dens(k+1,i,j) )
4536  enddo
4537  enddo
4538  enddo
4539 
4540  do j = js, je
4541  do i = is, ie
4542  do k = ks, ke
4543  momx(k,i,j) = 0.5_rp * u(k,i,j) * ( dens(k,i+1,j) + dens(k,i,j) )
4544  enddo
4545  enddo
4546  enddo
4547 
4548  do j = js, je
4549  do i = is, ie
4550  do k = ks, ke
4551  momy(k,i,j) = 0.5_rp * v(k,i,j) * ( dens(k,i,j+1) + dens(k,i,j) )
4552  enddo
4553  enddo
4554  enddo
4555 
4556  do j = js, je
4557  do i = is, ie
4558  do k = ks, ke
4559  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
4560  enddo
4561  enddo
4562  enddo
4563 
4564  return
4565  end subroutine mkinit_interporation
4566 
4567  !-----------------------------------------------------------------------------
4568  subroutine interporation_fact( &
4569  fact0, fact1, &
4570  idx0, idx1, &
4571  x, x_org, nx, &
4572  loop )
4573  implicit none
4574 
4575  real(RP), intent(out) :: fact0
4576  real(RP), intent(out) :: fact1
4577  integer, intent(out) :: idx0
4578  integer, intent(out) :: idx1
4579  real(RP), intent(in) :: x
4580  integer, intent(in) :: nx
4581  real(RP), intent(in) :: x_org(nx)
4582  logical, intent(in) :: loop
4583 
4584  real(RP) :: xwork
4585  integer :: i
4586 
4587  if ( x < x_org(1) ) then
4588  if ( loop ) then
4589  xwork = x_org(1) - ( x_org(2) - x_org(1) )**2 / ( x_org(3) - x_org(2) )
4590  fact0 = ( x_org(1) - x ) / ( x_org(1) - xwork )
4591  fact1 = ( x - xwork ) / ( x_org(1) - xwork )
4592  idx0 = nx
4593  idx1 = 1
4594  else
4595  fact0 = ( x_org(2) - x ) / ( x_org(2) - x_org(1) )
4596  fact1 = ( x - x_org(1) ) / ( x_org(2) - x_org(1) )
4597  idx0 = 1
4598  idx1 = 2
4599  end if
4600  else if ( x > x_org(nx) ) then
4601  if ( loop ) then
4602  xwork = x_org(nx) + ( x_org(nx) - x_org(nx-1) )**2 / ( x_org(nx-1) - x_org(nx-2) )
4603  fact0 = ( xwork - x ) / ( xwork - x_org(nx) )
4604  fact1 = ( x - x_org(nx) ) / ( xwork - x_org(nx) )
4605  idx0 = nx
4606  idx1 = 1
4607  else
4608  fact0 = ( x_org(nx) - x ) / ( x_org(nx) - x_org(nx-1) )
4609  fact1 = ( x - x_org(nx-1) ) / ( x_org(nx) - x_org(nx-1) )
4610  idx0 = nx-1
4611  idx1 = nx
4612  end if
4613  else
4614  do i = 2, nx
4615  if ( x <= x_org(i) ) then
4616  fact0 = ( x_org(i) - x ) / ( x_org(i) - x_org(i-1) )
4617  fact1 = ( x - x_org(i-1) ) / ( x_org(i) - x_org(i-1) )
4618  idx0 = i-1
4619  idx1 = i
4620  exit
4621  end if
4622  end do
4623  end if
4624 
4625  return
4626  end subroutine interporation_fact
4627  !-----------------------------------------------------------------------------
4629  subroutine mkinit_oceancouple
4630  implicit none
4631 
4632  if( io_l ) write(io_fid_log,*)
4633  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit OceanCouple] / Categ[preprocess] / Origin[SCALE-RM]'
4634 
4635  call flux_setup
4636 
4637  call ocean_setup
4638 
4639  return
4640  end subroutine mkinit_oceancouple
4641 
4642  !-----------------------------------------------------------------------------
4644  subroutine mkinit_landcouple
4645  implicit none
4646 
4647  if( io_l ) write(io_fid_log,*)
4648  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit LandCouple] / Categ[preprocess] / Origin[SCALE-RM]'
4649 
4650  call flux_setup
4651 
4652  call land_setup
4653 
4654  return
4655  end subroutine mkinit_landcouple
4656 
4657  !-----------------------------------------------------------------------------
4659  subroutine mkinit_urbancouple
4660  implicit none
4661 
4662  if( io_l ) write(io_fid_log,*)
4663  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit UrbanCouple] / Categ[preprocess] / Origin[SCALE-RM]'
4664 
4665  call flux_setup
4666 
4667  call urban_setup
4668 
4669  return
4670  end subroutine mkinit_urbancouple
4671 
4672  !-----------------------------------------------------------------------------
4674  subroutine mkinit_seabreeze
4675  use scale_rm_process, only: &
4676  prc_num_x
4677  use scale_landuse, only: &
4680  implicit none
4681 
4682  real(RP) :: dist
4683 
4684  integer :: i, j
4685  !---------------------------------------------------------------------------
4686 
4687  if( io_l ) write(io_fid_log,*)
4688  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit SEABREEZE] / Categ[preprocess] / Origin[SCALE-RM]'
4689 
4690  call flux_setup
4691 
4692  call land_setup
4693 
4694  call ocean_setup
4695 
4696  ! quartor size of domain
4697  dist = ( grid_cxg(imax*prc_num_x) - grid_cxg(1) ) / 8.0_rp
4698 
4699  ! make landuse conditions
4700  do j = js, je
4701  do i = is, ie
4702  if ( grid_cx(i) >= dist * 3.0_rp &
4703  .AND. grid_cx(i) <= dist * 5.0_rp ) then
4704  landuse_frac_land(i,j) = 1.0_rp
4705  else
4706  landuse_frac_land(i,j) = 0.0_rp
4707  endif
4708  enddo
4709  enddo
4710 
4711  call landuse_calc_fact
4712 
4713  return
4714  end subroutine mkinit_seabreeze
4715 
4716  !-----------------------------------------------------------------------------
4718  subroutine mkinit_heatisland
4719  use scale_rm_process, only: &
4720  prc_num_x
4721  use scale_landuse, only: &
4725  implicit none
4726 
4727  real(RP) :: dist
4728 
4729  integer :: i, j
4730  !---------------------------------------------------------------------------
4731 
4732  if( io_l ) write(io_fid_log,*)
4733  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit HEATISLAND] / Categ[preprocess] / Origin[SCALE-RM]'
4734 
4735  call flux_setup
4736 
4737  call land_setup
4738 
4739  call urban_setup
4740 
4741  ! 1/9 size of domain
4742  dist = ( grid_cxg(imax*prc_num_x) - grid_cxg(1) ) / 9.0_rp
4743 
4744  ! make landuse conditions
4745  do j = js, je
4746  do i = is, ie
4747  if ( grid_cx(i) >= dist * 4.0_rp &
4748  .AND. grid_cx(i) < dist * 5.0_rp ) then
4749  landuse_frac_land(i,j) = 1.0_rp
4750  landuse_frac_urban(i,j) = 1.0_rp
4751  else
4752  landuse_frac_land(i,j) = 1.0_rp
4753  landuse_frac_urban(i,j) = 0.0_rp
4754  endif
4755  enddo
4756  enddo
4757 
4758  call landuse_calc_fact
4759 
4760  return
4761  end subroutine mkinit_heatisland
4762 
4763  !-----------------------------------------------------------------------------
4765  subroutine mkinit_grayzone
4766  implicit none
4767 
4768  real(RP) :: RHO(ka)
4769  real(RP) :: VELX(ka)
4770  real(RP) :: VELY(ka)
4771  real(RP) :: POTT(ka)
4772  real(RP) :: QV(ka)
4773 
4774  real(RP) :: PERTURB_AMP = 0.0_rp
4775  integer :: RANDOM_LIMIT = 0
4776  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
4777  ! 1 -> petrurbation for pt
4778  ! 2 -> perturbation for u, v, w
4779 
4780  namelist / param_mkinit_grayzone / &
4781  perturb_amp, &
4782  random_limit, &
4783  random_flag
4784 
4785  integer :: ierr
4786  integer :: k, i, j
4787  !---------------------------------------------------------------------------
4788 
4789  if( io_l ) write(io_fid_log,*)
4790  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit GRAYZONE] / Categ[preprocess] / Origin[SCALE-RM]'
4791 
4792  !--- read namelist
4793  rewind(io_fid_conf)
4794  read(io_fid_conf,nml=param_mkinit_grayzone,iostat=ierr)
4795 
4796  if( ierr < 0 ) then !--- missing
4797  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
4798  elseif( ierr > 0 ) then !--- fatal error
4799  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_GRAYZONE. Check!'
4800  call prc_mpistop
4801  endif
4802  if( io_l ) write(io_fid_log,nml=param_mkinit_grayzone)
4803 
4804  call read_sounding( rho, velx, vely, pott, qv ) ! (out)
4805 
4806 ! do j = JS, JE
4807 ! do i = IS, IE
4808  do j = 1, ja
4809  do i = 1, ia
4810  do k = ks, ke
4811  dens(k,i,j) = rho(k)
4812 ! MOMZ(k,i,j) = 0.0_RP
4813 ! MOMX(k,i,j) = RHO(k) * VELX(k)
4814 ! MOMY(k,i,j) = RHO(k) * VELY(k)
4815 
4816 ! RHOT(k,i,j) = RHO(k) * POTT(k)
4817 
4818  qtrc(k,i,j,i_qv) = qv(k)
4819  enddo
4820  enddo
4821  enddo
4822 
4823  do j = js, je
4824  do i = is, ie
4825  dens( 1:ks-1,i,j) = dens(ks,i,j)
4826  dens(ke+1:ka, i,j) = dens(ke,i,j)
4827  enddo
4828  enddo
4829 
4830  call random_get(rndm) ! make random
4831  do j = js, je
4832  do i = is, ie
4833  do k = ks, ke
4834  if ( random_flag == 2 .and. k <= random_limit ) then ! below initial cloud top
4835  momz(k,i,j) = ( 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4836  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
4837  else
4838  momz(k,i,j) = 0.0_rp
4839  endif
4840  enddo
4841  enddo
4842  enddo
4843 
4844  call random_get(rndm) ! make random
4845  do j = js, je
4846  do i = is, ie
4847  do k = ks, ke
4848  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
4849  momx(k,i,j) = ( velx(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4850  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
4851  else
4852  momx(k,i,j) = velx(k) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
4853  endif
4854  enddo
4855  enddo
4856  enddo
4857 
4858  call random_get(rndm) ! make random
4859  do j = js, je
4860  do i = is, ie
4861  do k = ks, ke
4862  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
4863  momy(k,i,j) = ( vely(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4864  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
4865  else
4866  momy(k,i,j) = vely(k) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
4867  endif
4868  enddo
4869  enddo
4870  enddo
4871 
4872  call random_get(rndm) ! make random
4873  do j = js, je
4874  do i = is, ie
4875  do k = ks, ke
4876  if ( random_flag == 1 .and. k <= random_limit ) then ! below initial cloud top
4877  rhot(k,i,j) = ( pott(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4878  * dens(k,i,j)
4879  else
4880  rhot(k,i,j) = pott(k) * dens(k,i,j)
4881  endif
4882  enddo
4883  enddo
4884  enddo
4885 
4886  return
4887  end subroutine mkinit_grayzone
4888  !-----------------------------------------------------------------------------
4890  subroutine mkinit_boxaero
4891 
4892  implicit none
4893 
4894  real(RP) :: init_dens = 1.12_rp ![kg/m3]
4895  real(RP) :: init_temp = 298.18_rp ![K]
4896  real(RP) :: init_pres = 1.e+5_rp ![Pa]
4897  real(RP) :: init_ssliq = 0.01_rp ![%]
4898 
4899  namelist / param_mkinit_boxaero / &
4900  init_dens, &
4901  init_temp, &
4902  init_pres, &
4903  init_ssliq
4904 
4905  real(RP) :: qsat
4906  integer :: i, j, k, ierr
4907 
4908  if ( aetracer_type /= 'KAJINO13' ) then
4909  if( io_l ) write(io_fid_log,*) '+++ For [Box model of aerosol],'
4910  if( io_l ) write(io_fid_log,*) '+++ AETRACER_TYPE should be KAJINO13. Stop!'
4911  call prc_mpistop
4912  endif
4913 
4914  if( io_l ) write(io_fid_log,*)
4915  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit Box model of aerosol] / Categ[preprocess] / Origin[SCALE-RM]'
4916 
4917  !--- read namelist
4918  rewind(io_fid_conf)
4919  read(io_fid_conf,nml=param_mkinit_boxaero,iostat=ierr)
4920  if( ierr < 0 ) then !--- missing
4921  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
4922  elseif( ierr > 0 ) then !--- fatal error
4923  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_BOXAERO. Check!'
4924  call prc_mpistop
4925  endif
4926  if( io_lnml ) write(io_fid_log,nml=param_mkinit_boxaero)
4927 
4928  qtrc(:,:,:,:) = 0.0_rp
4929  do k = ks, ke
4930  do i = is, ie
4931  do j = js, je
4932  dens(k,i,j) = init_dens
4933  momx(k,i,j) = 0.0_rp
4934  momy(k,i,j) = 0.0_rp
4935  momz(k,i,j) = 0.0_rp
4936  pott(k,i,j) = init_temp * ( p00/init_pres )**( rdry/cpdry )
4937  rhot(k,i,j) = dens(k,i,j) * pott(k,i,j)
4938  call saturation_pres2qsat_all( qsat,init_temp,init_pres )
4939  qtrc(k,i,j,i_qv) = ( init_ssliq + 1.0_rp )*qsat
4940  enddo
4941  enddo
4942  enddo
4943 
4944  if( aetracer_type == 'KAJINO13' ) then
4945  call aerosol_setup
4946  endif
4947 
4948  end subroutine mkinit_boxaero
4949  !-----------------------------------------------------------------------------
4951  subroutine mkinit_warmbubbleaero
4952  implicit none
4953 
4954  ! Surface state
4955  real(RP) :: SFC_THETA ! surface potential temperature [K]
4956  real(RP) :: SFC_PRES ! surface pressure [Pa]
4957  real(RP) :: SFC_RH = 80.0_rp ! surface relative humidity [%]
4958  ! Environment state
4959  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
4960  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
4961  real(RP) :: ENV_RH = 80.0_rp ! Relative Humidity of environment [%]
4962  real(RP) :: ENV_L1_ZTOP = 1.e3_rp ! top height of the layer1 (constant THETA) [m]
4963  real(RP) :: ENV_L2_ZTOP = 14.e3_rp ! top height of the layer2 (small THETA gradient) [m]
4964  real(RP) :: ENV_L2_TLAPS = 4.e-3_rp ! Lapse rate of THETA in the layer2 (small THETA gradient) [K/m]
4965  real(RP) :: ENV_L3_TLAPS = 3.e-2_rp ! Lapse rate of THETA in the layer3 (large THETA gradient) [K/m]
4966  ! Bubble
4967  real(RP) :: BBL_THETA = 1.0_rp ! extremum of temperature in bubble [K]
4968 
4969  namelist / param_mkinit_warmbubble / &
4970  sfc_theta, &
4971  sfc_pres, &
4972  env_u, &
4973  env_v, &
4974  env_rh, &
4975  env_l1_ztop, &
4976  env_l2_ztop, &
4977  env_l2_tlaps, &
4978  env_l3_tlaps, &
4979  bbl_theta
4980 
4981  integer :: ierr
4982  integer :: k, i, j
4983  !---------------------------------------------------------------------------
4984 
4985  if( io_l ) write(io_fid_log,*)
4986  if( io_l ) write(io_fid_log,*) '++++++ Module[mkinit WARMBUBBLEAERO] / Categ[preprocess] / Origin[SCALE-RM]'
4987 
4988  sfc_theta = thetastd
4989  sfc_pres = pstd
4990 
4991  !--- read namelist
4992  rewind(io_fid_conf)
4993  read(io_fid_conf,nml=param_mkinit_warmbubble,iostat=ierr)
4994 
4995  if( ierr < 0 ) then !--- missing
4996  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
4997  elseif( ierr > 0 ) then !--- fatal error
4998  write(*,*) 'xxx Not appropriate names in namelist PARAM_MKINIT_WARMBUBBLE. Check!'
4999  call prc_mpistop
5000  endif
5001  if( io_lnml ) write(io_fid_log,nml=param_mkinit_warmbubble)
5002 
5003  ! calc in dry condition
5004  pres_sfc(1,1,1) = sfc_pres
5005  pott_sfc(1,1,1) = sfc_theta
5006  qv_sfc(1,1,1) = 0.0_rp
5007  qc_sfc(1,1,1) = 0.0_rp
5008 
5009  do k = ks, ke
5010  if ( grid_cz(k) <= env_l1_ztop ) then ! Layer 1
5011  pott(k,1,1) = sfc_theta
5012  elseif( grid_cz(k) < env_l2_ztop ) then ! Layer 2
5013  pott(k,1,1) = pott(k-1,1,1) + env_l2_tlaps * ( grid_cz(k)-grid_cz(k-1) )
5014  else ! Layer 3
5015  pott(k,1,1) = pott(k-1,1,1) + env_l3_tlaps * ( grid_cz(k)-grid_cz(k-1) )
5016  endif
5017  qv(k,1,1) = 0.0_rp
5018  qc(k,1,1) = 0.0_rp
5019  enddo
5020 
5021  ! make density & pressure profile in dry condition
5022  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
5023  temp(:,1,1), & ! [OUT]
5024  pres(:,1,1), & ! [OUT]
5025  pott(:,1,1), & ! [IN]
5026  qv(:,1,1), & ! [IN]
5027  qc(:,1,1), & ! [IN]
5028  temp_sfc(1,1,1), & ! [OUT]
5029  pres_sfc(1,1,1), & ! [IN]
5030  pott_sfc(1,1,1), & ! [IN]
5031  qv_sfc(1,1,1), & ! [IN]
5032  qc_sfc(1,1,1) ) ! [IN]
5033 
5034  ! calc QV from RH
5035  call saturation_pres2qsat_all( qsat_sfc(1,1,1), temp_sfc(1,1,1), pres_sfc(1,1,1) )
5036  call saturation_pres2qsat_all( qsat(:,1,1), temp(:,1,1), pres(:,1,1) )
5037 
5038  qv_sfc(1,1,1) = sfc_rh * 1.e-2_rp * qsat_sfc(1,1,1)
5039  do k = ks, ke
5040  if ( grid_cz(k) <= env_l1_ztop ) then ! Layer 1
5041  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
5042  elseif( grid_cz(k) <= env_l2_ztop ) then ! Layer 2
5043  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
5044  else ! Layer 3
5045  qv(k,1,1) = 0.0_rp
5046  endif
5047  enddo
5048 
5049  ! make density & pressure profile in moist condition
5050  call hydrostatic_buildrho( dens(:,1,1), & ! [OUT]
5051  temp(:,1,1), & ! [OUT]
5052  pres(:,1,1), & ! [OUT]
5053  pott(:,1,1), & ! [IN]
5054  qv(:,1,1), & ! [IN]
5055  qc(:,1,1), & ! [IN]
5056  temp_sfc(1,1,1), & ! [OUT]
5057  pres_sfc(1,1,1), & ! [IN]
5058  pott_sfc(1,1,1), & ! [IN]
5059  qv_sfc(1,1,1), & ! [IN]
5060  qc_sfc(1,1,1) ) ! [IN]
5061 
5062  do j = js, je
5063  do i = is, ie
5064  do k = ks, ke
5065  dens(k,i,j) = dens(k,1,1)
5066  momz(k,i,j) = 0.0_rp
5067  momx(k,i,j) = env_u * dens(k,i,j)
5068  momy(k,i,j) = env_v * dens(k,i,j)
5069 
5070  ! make warm bubble
5071  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
5072 
5073  qtrc(k,i,j,i_qv) = qv(k,1,1)
5074  enddo
5075  enddo
5076  enddo
5077 
5078  call flux_setup
5079 
5080  if( aetracer_type == 'KAJINO13' ) then
5081  call aerosol_setup
5082  endif
5083 
5084  return
5085  end subroutine mkinit_warmbubbleaero
5086 
5087  !-----------------------------------------------------------------------------
5089  subroutine mkinit_real
5090  use mod_realinput, only: &
5091  realinput_atmos, &
5093  implicit none
5094 
5095  call realinput_atmos( flg_intrp )
5096 
5097  call realinput_surface
5098 
5099  call flux_setup
5100 
5101  return
5102  end subroutine mkinit_real
5103 
5104 end module mod_mkinit
5105 !-------------------------------------------------------------------------------
integer, public imax
of computational cells: x
integer, public prc_num_x
x length of 2D processor topology
integer, public is
start point of inner domain: x, local
module GTOOL_FILE
Definition: gtool_file.f90:17
real(rp), dimension(:,:,:), allocatable, target, public momz
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:59
integer, public je
end point of inner domain: y, local
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
module ATMOSPHERE / Saturation adjustment
real(rp), dimension(:,:,:), allocatable, public ocean_sfc_albedo
ocean surface albedo [0-1]
integer, public const_i_lw
long-wave radiation index
Definition: scale_const.F90:98
subroutine, public prc_mpistop
Abort MPI.
real(rp), dimension(:,:), allocatable, public urban_qc
logical, public land_restart_output
output restart file?
integer, parameter, public i_oceancouple
Definition: mod_mkinit.f90:126
real(rp), dimension(:,:,:), allocatable, target, public rhot
real(rp), dimension(:), allocatable, public grid_cxg
center coordinate [m]: x, global
character(len=h_short), public aetracer_type
integer, public qqe
module Atmosphere / Physics Cloud Microphysics
module URBAN driver
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0e
ocean surface roughness length for vapor [m]
real(rp), dimension(:), allocatable, public grid_cz
center coordinate [m]: z, local=global
real(rp), dimension(:,:), allocatable, public ocean_temp
temperature at uppermost ocean layer [K]
integer, parameter, public i_boxaero
Definition: mod_mkinit.f90:139
subroutine, public ocean_vars_restart_write
Write ocean restart.
subroutine land_setup
Land setup.
module ATMOSPHERIC Variables
real(rp), dimension(:,:,:), allocatable, target, public momx
subroutine, public atmos_vars_restart_write
Write restart of atmospheric variables.
module STDIO
Definition: scale_stdio.F90:12
subroutine ocean_setup
Ocean setup.
subroutine, public landuse_calc_fact
integer, public qaes
integer, parameter, public i_dycoms2_rf02_dns
Definition: mod_mkinit.f90:134
integer, public ke
end point of inner domain: z, local
subroutine, public filegetshape(dims, basename, varname, myrank, single)
Definition: gtool_file.f90:970
integer, public qa
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
integer, parameter, public i_turbulence
Definition: mod_mkinit.f90:112
subroutine, public mkinit
Driver.
Definition: mod_mkinit.f90:347
integer, parameter, public i_urbancouple
Definition: mod_mkinit.f90:127
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
integer, parameter, public i_gravitywave
Definition: mod_mkinit.f90:110
real(rp) function faero(f0, r0, x, alpha, rhoa)
real(rp), dimension(:,:), allocatable, public urban_tb
integer, parameter, public i_grayzone
Definition: mod_mkinit.f90:138
real(rp), dimension(:,:,:), allocatable, target, public dens
subroutine interporation_fact(fact0, fact1, idx0, idx1, x, x_org, nx, loop)
real(rp), dimension(:,:), allocatable, public landuse_frac_urban
urban fraction
module URBAN Variables
real(rp), dimension(:,:), allocatable, public urban_raing
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_lw_up
integer, parameter, public i_landcouple
Definition: mod_mkinit.f90:125
real(rp), dimension(:,:), allocatable, public atmos_phy_mp_sflx_rain
character(len=h_short), public tracer_type
module TRACER / suzuki10
real(rp), dimension(:,:), allocatable, public urban_uc
real(rp), dimension(:), allocatable, public grid_fx
face coordinate [m]: x, local
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_sw_up
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_lw_dn
logical, public atmos_restart_output
output restart file?
subroutine, public urban_surface_set(countup)
Set surface boundary to other model.
module Atmosphere / Physics Radiation
subroutine, public urban_vars_restart_write
Write urban restart.
module grid index
real(rp), dimension(:,:), allocatable, public land_sfc_temp
land surface skin temperature [K]
real(rp), public const_mvap
mass weight (water vapor) [g/mol]
Definition: scale_const.F90:64
subroutine, public atmos_surface_get
Get surface boundary condition.
module TRACER
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_sw_up
real(rp), dimension(:,:), allocatable, public urban_tr
module REAL input
module ATMOSPHERE / Typical vertical profile
integer, parameter, public i_seabreeze
Definition: mod_mkinit.f90:131
integer, public ia
of x whole cells (local, with HALO)
real(rp), dimension(:,:,:), allocatable, public urban_tgl
integer function, public io_get_available_fid()
search & get available file ID
subroutine rect_setup
Bubble.
Definition: mod_mkinit.f90:599
integer, dimension(:), allocatable, public nkap
module GRID (real space)
module LANDUSE
integer, public ka
of z whole cells (local, with HALO)
integer, public ae_ctg
integer, public i_qv
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_sw_dn
module OCEAN driver
integer, parameter, public dp
Definition: dc_types.f90:27
integer, parameter, public i_warmbubble
Definition: mod_mkinit.f90:115
subroutine, public ocean_surface_set(countup)
Put surface boundary to other model.
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:93
integer, public qqa
integer, public gas_ctg
subroutine read_sounding(DENS, VELX, VELY, POTT, QV)
Read sounding data from file.
real(rp), dimension(:), allocatable, public grid_fz
face coordinate [m]: z, local=global
module ATMOSPHERE driver
integer, parameter, public i_mountainwave
Definition: mod_mkinit.f90:113
subroutine, public landuse_write
Write landuse data.
character(len=h_short), dimension(:), allocatable, public aq_name
real(rp), dimension(:,:,:), allocatable, public land_temp
temperature of each soil layer [K]
module COMMUNICATION
Definition: scale_comm.F90:23
module ATMOSPHERE / Hydrostatic barance
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
integer, public js
start point of inner domain: y, local
integer, public i_nc
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_lw_dn
real(rp), dimension(:,:), allocatable, public urban_roff
module PROCESS
real(rp), dimension(:,:,:), allocatable, public atmos_phy_ae_ccn
subroutine, public realinput_atmos(flg_intrp)
integer, parameter, public i_real
Definition: mod_mkinit.f90:136
real(rp), dimension(:,:,:), allocatable, public land_sfc_albedo
land surface albedo [0-1]
integer, dimension(:), allocatable, public nsiz
integer, parameter, public i_heatisland
Definition: mod_mkinit.f90:132
subroutine, public log(type, message)
Definition: dc_log.f90:133
real(rp), public const_lhv
latent heat of vaporizaion for use
Definition: scale_const.F90:77
integer, parameter, public i_triplecouple
Definition: mod_mkinit.f90:128
module LAND Variables
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
module CONSTANT
Definition: scale_const.F90:14
real(rp), dimension(:,:,:), allocatable, public urban_sfc_albedo
real(rp), dimension(:), allocatable, public aq_cv
CV for each hydrometeors [J/kg/K].
integer, parameter, public i_planestate
Definition: mod_mkinit.f90:105
integer, parameter, public i_khwave
Definition: mod_mkinit.f90:111
integer, public ks
start point of inner domain: z, local
real(rp), dimension(:), allocatable, public grid_cx
center coordinate [m]: x, local
integer, parameter, public i_squallline
Definition: mod_mkinit.f90:117
real(rp), dimension(:,:), allocatable, public urban_tc
real(rp), dimension(:,:,:), allocatable, target, public momy
logical, public urban_restart_output
output restart file?
module GRID (cartesian)
real(rp), dimension(:,:), allocatable, public urban_rainr
subroutine, public random_get(var)
Get random number.
real(dp), parameter, public const_undef8
undefined value (REAL8)
Definition: scale_const.F90:42
real(rp), dimension(:,:,:), allocatable, public land_water
moisture of each soil layer [m3/m3]
real(rp), dimension(:,:,:,:), allocatable, public atmos_phy_rd_sflx_downall
integer, public mkinit_type
Definition: mod_mkinit.f90:102
module RM PROCESS
real(rp), dimension(:,:), allocatable, public atmos_phy_mp_sflx_snow
integer, parameter, public i_rico
Definition: mod_mkinit.f90:121
module profiler
Definition: scale_prof.F90:10
integer, public ie
end point of inner domain: x, local
integer, parameter, public i_ignore
Definition: mod_mkinit.f90:103
integer, public ic_mix
module ATMOSPHERE / Thermodynamics
integer, parameter, public i_lambwave
Definition: mod_mkinit.f90:109
integer, parameter, public i_supercell
Definition: mod_mkinit.f90:116
subroutine aerosol_activation_kajino13(c_kappa, super, temp_k, ia_m0, ia_m2, ia_m3, n_atr, n_siz_max, n_kap_max, n_ctg, n_siz, n_kap, d_ct, aerosol_procs, aerosol_activ)
Definition: mod_mkinit.f90:969
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0m
ocean surface roughness length for momentum [m]
integer, public qqs
logical, public io_lnml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:60
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_sw_dn
integer, public const_i_sw
short-wave radiation index
Definition: scale_const.F90:99
module PRECISION
subroutine diag_ds(m0, m2, m3, dg, sg, dm2)
integer, public n_atr
integer, parameter, public i_dycoms2_rf01
Definition: mod_mkinit.f90:119
integer, parameter, public i_interporation
Definition: mod_mkinit.f90:123
real(rp), public const_pi
pi
Definition: scale_const.F90:34
real(rp), dimension(:,:), allocatable, public urban_tg
real(rp), dimension(:,:), allocatable, public ocean_sfc_temp
ocean surface skin temperature [K]
real(rp), dimension(:,:), allocatable, public urban_sfc_temp
real(rp), dimension(:,:,:), allocatable, public urban_trl
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
integer, parameter, public i_warmbubbleaero
Definition: mod_mkinit.f90:140
integer, parameter, public i_wk1982
Definition: mod_mkinit.f90:118
module INITIAL
Definition: mod_mkinit.f90:31
integer, parameter, public i_bubblecouple
Definition: mod_mkinit.f90:129
module RANDOM
integer, parameter, public i_tracerbubble
Definition: mod_mkinit.f90:106
logical, public ocean_restart_output
output restart file?
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
subroutine, public land_vars_restart_write
Write land restart.
integer, parameter, public i_coldbubble
Definition: mod_mkinit.f90:107
subroutine, public realinput_surface
integer, parameter, public i_cavityflow
Definition: mod_mkinit.f90:142
integer, parameter, public rp
real(rp), dimension(:,:,:), allocatable, public urban_tbl
module ATMOSPHERE / Physics Aerosol Microphysics
integer, public qaee
module LAND driver
real(rp), dimension(:,:), allocatable, public landuse_frac_land
land fraction
subroutine, public mkinit_setup
Setup.
Definition: mod_mkinit.f90:221
integer, parameter, public i_dycoms2_rf02
Definition: mod_mkinit.f90:120
subroutine, public land_surface_set(countup)
Put surface boundary to other model.
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0h
ocean surface roughness length for heat [m]
module OCEAN Variables
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:92
real(rp), dimension(:), allocatable, public grid_cy
center coordinate [m]: y, local
real(rp), dimension(:,:), allocatable, public urban_rainb
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_lw_up
subroutine urban_setup
Urban setup.
subroutine flux_setup
flux setup
integer, public i_qc
real(rp), dimension(:,:,:,:), allocatable, target, public qtrc
integer, public ja
of y whole cells (local, with HALO)
real(rp), dimension(:), allocatable, public grid_fy
face coordinate [m]: y, local