SCALE-RM
mod_mkinit.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module mod_mkinit
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
21  use scale_tracer
23 
24  use scale_prc, only: &
25  prc_abort
26  use scale_const, only: &
27  pi => const_pi, &
28  grav => const_grav, &
29  pstd => const_pstd, &
30  rdry => const_rdry, &
31  cpdry => const_cpdry, &
32  p00 => const_pre00, &
33  epsvap => const_epsvap
34  use scale_random, only: &
36  use scale_comm_cartesc, only: &
37  comm_vars8, &
38  comm_wait
39  use scale_atmos_grid_cartesc, only: &
40  cz => atmos_grid_cartesc_cz, &
41  cx => atmos_grid_cartesc_cx, &
42  cy => atmos_grid_cartesc_cy, &
43  fz => atmos_grid_cartesc_fz, &
44  fx => atmos_grid_cartesc_fx, &
45  fy => atmos_grid_cartesc_fy, &
46  cxg => atmos_grid_cartesc_cxg, &
47  fxg => atmos_grid_cartesc_fxg, &
50  real_cz => atmos_grid_cartesc_real_cz, &
51  real_fz => atmos_grid_cartesc_real_fz, &
53  use scale_atmos_profile, only: &
54  profile_isa => atmos_profile_isa
55  use scale_atmos_hydrometeor, only: &
56  hydrometeor_lhv => atmos_hydrometeor_lhv
57  use scale_atmos_hydrostatic, only: &
58  hydrostatic_buildrho => atmos_hydrostatic_buildrho, &
59  hydrostatic_buildrho_atmos => atmos_hydrostatic_buildrho_atmos, &
60  hydrostatic_buildrho_bytemp => atmos_hydrostatic_buildrho_bytemp
61  use scale_atmos_saturation, only: &
62  saturation_psat_all => atmos_saturation_psat_all, &
63  saturation_pres2qsat_all => atmos_saturation_pres2qsat_all, &
64  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq
65  use mod_atmos_vars, only: &
66  dens, &
67  momx, &
68  momy, &
69  momz, &
70  rhot, &
71  qtrc
72  use mod_atmos_phy_ae_vars, only: &
73  ccn => atmos_phy_ae_ccn
74  !-----------------------------------------------------------------------------
75  implicit none
76  private
77  !-----------------------------------------------------------------------------
78  !
79  !++ Public procedure
80  !
81  public :: mkinit_setup
82  public :: mkinit
83 
84  !-----------------------------------------------------------------------------
85  !
86  !++ Public parameters & variables
87  !
88  integer, public :: mkinit_type = -1
89  integer, public, parameter :: i_ignore = 0
90 
91  integer, public, parameter :: i_planestate = 1
92  integer, public, parameter :: i_tracerbubble = 2
93  integer, public, parameter :: i_coldbubble = 3
94 
95  integer, public, parameter :: i_lambwave = 4
96  integer, public, parameter :: i_gravitywave = 5
97  integer, public, parameter :: i_khwave = 6
98  integer, public, parameter :: i_turbulence = 7
99  integer, public, parameter :: i_mountainwave = 8
100 
101  integer, public, parameter :: i_warmbubble = 9
102  integer, public, parameter :: i_supercell = 10
103  integer, public, parameter :: i_squallline = 11
104  integer, public, parameter :: i_wk1982 = 12
105  integer, public, parameter :: i_dycoms2_rf01 = 13
106  integer, public, parameter :: i_dycoms2_rf02 = 14
107  integer, public, parameter :: i_rico = 15
108 
109  integer, public, parameter :: i_interporation = 16
110 
111  integer, public, parameter :: i_landcouple = 17
112  integer, public, parameter :: i_oceancouple = 18
113  integer, public, parameter :: i_urbancouple = 19
114  integer, public, parameter :: i_triplecouple = 20
115  integer, public, parameter :: i_bubblecouple = 21
116 
117  integer, public, parameter :: i_seabreeze = 22
118  integer, public, parameter :: i_heatisland = 23
119 
120  integer, public, parameter :: i_dycoms2_rf02_dns = 24
121 
122  integer, public, parameter :: i_real = 25
123 
124  integer, public, parameter :: i_grayzone = 26
125  integer, public, parameter :: i_boxaero = 27
126  integer, public, parameter :: i_warmbubbleaero = 28
127 
128  integer, public, parameter :: i_cavityflow = 29
129  integer, public, parameter :: i_barocwave = 30
130  integer, public, parameter :: i_bomex = 31
131 
132  !-----------------------------------------------------------------------------
133  !
134  !++ Private procedure
135  !
136  private :: bubble_setup
137  private :: sbmaero_setup
138  private :: aerosol_setup
139 
140  private :: mkinit_planestate
141  private :: mkinit_tracerbubble
142  private :: mkinit_coldbubble
143  private :: mkinit_lambwave
144  private :: mkinit_gravitywave
145  private :: mkinit_khwave
146  private :: mkinit_turbulence
147  private :: mkinit_cavityflow
148  private :: mkinit_mountainwave
149  private :: mkinit_barocwave
150 
151  private :: mkinit_warmbubble
152  private :: mkinit_supercell
153  private :: mkinit_squallline
154  private :: mkinit_wk1982
155  private :: mkinit_dycoms2_rf01
156  private :: mkinit_dycoms2_rf02
157  private :: mkinit_rico
158  private :: mkinit_bomex
159 
160  private :: mkinit_landcouple
161  private :: mkinit_oceancouple
162  private :: mkinit_urbancouple
163  private :: mkinit_seabreeze
164  private :: mkinit_heatisland
165 
166  private :: mkinit_dycoms2_rf02_dns
167 
168  private :: mkinit_real
169 
170  private :: mkinit_grayzone
171 
172  private :: mkinit_boxaero
173  private :: mkinit_warmbubbleaero
174 
175  !-----------------------------------------------------------------------------
176  !
177  !++ Private parameters & variables
178  !
179  real(RP), private, parameter :: thetastd = 300.0_rp ! [K]
180 
181  real(RP), private, allocatable :: pres (:,:,:) ! pressure [Pa]
182  real(RP), private, allocatable :: temp (:,:,:) ! temperature [K]
183  real(RP), private, allocatable :: pott (:,:,:) ! potential temperature [K]
184  real(RP), private, allocatable :: qdry (:,:,:) ! dry air mass ratio [kg/kg]
185  real(RP), private, allocatable :: qsat (:,:,:) ! satulated water vapor [kg/kg]
186  real(RP), private, allocatable :: qv (:,:,:) ! water vapor [kg/kg]
187  real(RP), private, allocatable :: qc (:,:,:) ! cloud water [kg/kg]
188  real(RP), private, allocatable :: nc (:,:,:) ! cloud water number density [1/kg]
189  real(RP), private, allocatable :: velx (:,:,:) ! velocity u [m/s]
190  real(RP), private, allocatable :: vely (:,:,:) ! velocity v [m/s]
191 
192  real(RP), private, allocatable :: pres_sfc(:,:) ! surface pressure [Pa]
193  real(RP), private, allocatable :: temp_sfc(:,:) ! surface temperature [K]
194  real(RP), private, allocatable :: pott_sfc(:,:) ! surface potential temperature [K]
195  real(RP), private, allocatable :: psat_sfc(:,:) ! surface satulated water pressure [Pa]
196  real(RP), private, allocatable :: qsat_sfc(:,:) ! surface satulated water vapor [kg/kg]
197  real(RP), private, allocatable :: qv_sfc (:,:) ! surface water vapor [kg/kg]
198  real(RP), private, allocatable :: qc_sfc (:,:) ! surface cloud water [kg/kg]
199 
200  real(RP), private, allocatable :: rndm (:,:,:) ! random number (0-1)
201  real(RP), private, allocatable, target :: bubble (:,:,:) ! bubble factor (0-1)
202  real(RP), private, allocatable, target :: rect (:,:,:) ! rectangle factor (0-1)
203  real(RP), private, allocatable :: gan (:) ! gamma factor (0-1)
204 
205  !-----------------------------------------------------------------------------
206 contains
207  !-----------------------------------------------------------------------------
209  subroutine mkinit_setup
210  implicit none
211 
212  character(len=H_SHORT) :: MKINIT_initname = 'NONE'
213 
214  namelist / param_mkinit / &
215  mkinit_initname
216 
217  integer :: ierr
218  !---------------------------------------------------------------------------
219 
220  log_newline
221  log_info("MKINIT_setup",*) 'Setup'
222 
223  !--- read namelist
224  rewind(io_fid_conf)
225  read(io_fid_conf,nml=param_mkinit,iostat=ierr)
226  if( ierr < 0 ) then !--- missing
227  log_info("MKINIT_setup",*) 'Not found namelist. Default used.'
228  elseif( ierr > 0 ) then !--- fatal error
229  log_error("MKINIT_setup",*) 'Not appropriate names in namelist PARAM_MKINIT. Check!'
230  call prc_abort
231  endif
232  log_nml(param_mkinit)
233 
234  allocate( pres(ka,ia,ja) )
235  allocate( temp(ka,ia,ja) )
236  allocate( pott(ka,ia,ja) )
237  allocate( qdry(ka,ia,ja) )
238  allocate( qsat(ka,ia,ja) )
239  allocate( qv(ka,ia,ja) )
240  allocate( qc(ka,ia,ja) )
241  allocate( nc(ka,ia,ja) )
242  allocate( velx(ka,ia,ja) )
243  allocate( vely(ka,ia,ja) )
244 
245  allocate( pres_sfc(ia,ja) )
246  allocate( temp_sfc(ia,ja) )
247  allocate( pott_sfc(ia,ja) )
248  allocate( psat_sfc(ia,ja) )
249  allocate( qsat_sfc(ia,ja) )
250  allocate( qv_sfc(ia,ja) )
251  allocate( qc_sfc(ia,ja) )
252 
253  allocate( rndm(ka,ia,ja) )
254  allocate( bubble(ka,ia,ja) )
255  allocate( rect(ka,ia,ja) )
256 
257 
258  select case(trim(mkinit_initname))
259  case('NONE')
261  case('PLANESTATE')
263  case('TRACERBUBBLE')
265  case('COLDBUBBLE')
267  call bubble_setup
268  case('LAMBWAVE')
270  call bubble_setup
271  case('GRAVITYWAVE')
273  call bubble_setup
274  case('KHWAVE')
276  case('TURBULENCE')
278  case('MOUNTAINWAVE')
280  call bubble_setup
281  case('WARMBUBBLE')
283  call bubble_setup
284  case('SUPERCELL')
286  call bubble_setup
287  case('SQUALLLINE')
289  case('WK1982')
291  call bubble_setup
292  case('DYCOMS2_RF01')
294  case('DYCOMS2_RF02')
296  case('RICO')
298  case('BOMEX')
300  case('INTERPORATION')
302  case('LANDCOUPLE')
304  case('OCEANCOUPLE')
306  case('URBANCOUPLE')
308  case('TRIPLECOUPLE')
310  case('BUBBLECOUPLE')
312  call bubble_setup
313  case('SEABREEZE')
315  case('HEATISLAND')
317  case('DYCOMS2_RF02_DNS')
319  case('REAL')
321  case('GRAYZONE')
323  case('BOXAERO')
325  case('WARMBUBBLEAERO')
327  call bubble_setup
328  case('CAVITYFLOW')
330  case('BAROCWAVE')
332  case default
333  log_error("MKINIT_setup",*) 'Unsupported TYPE:', trim(mkinit_initname)
334  call prc_abort
335  endselect
336 
337  return
338  end subroutine mkinit_setup
339 
340  !-----------------------------------------------------------------------------
342  subroutine mkinit( output )
343  use scale_const, only: &
345  use scale_atmos_hydrometeor, only: &
346  n_hyd, &
347  i_hc
348  use mod_atmos_phy_mp_vars, only: &
349  qa_mp, &
350  qs_mp, &
351  qe_mp
352  use mod_atmos_admin, only: &
354  use mod_atmos_phy_mp_driver, only: &
356  implicit none
357  logical, intent(out) :: output
358 
359  real(RP) :: QHYD(ka,ia,ja,n_hyd)
360  real(RP) :: QNUM(ka,ia,ja,n_hyd)
361 
362  logical :: convert_qtrc
363  !---------------------------------------------------------------------------
364 
365  if ( mkinit_type == i_ignore ) then
366  log_newline
367  log_progress(*) 'skip making initial data'
368  output = .false.
369  else
370  log_newline
371  log_progress(*) 'start making initial data'
372 
373  !--- Initialize variables
374  pres(:,:,:) = const_undef8
375  temp(:,:,:) = const_undef8
376  pott(:,:,:) = const_undef8
377  qsat(:,:,:) = const_undef8
378  velx(:,:,:) = const_undef8
379  vely(:,:,:) = const_undef8
380 
381  rndm(:,:,:) = const_undef8
382 
383  pres_sfc(:,:) = const_undef8
384  temp_sfc(:,:) = const_undef8
385  pott_sfc(:,:) = const_undef8
386  psat_sfc(:,:) = const_undef8
387  qsat_sfc(:,:) = const_undef8
388 
389  qv(:,:,:) = 0.0_rp
390  qc(:,:,:) = 0.0_rp
391  nc(:,:,:) = 0.0_rp
392  qv_sfc(:,:) = 0.0_rp
393  qc_sfc(:,:) = 0.0_rp
394 
395 !OCL XFILL
396  qtrc(:,:,:,:) = 0.0_rp
397 !OCL XFILL
398  qhyd(:,:,:,:) = 0.0_rp
399 !OCL XFILL
400  qnum(:,:,:,:) = 0.0_rp
401 
402  call prof_rapstart('_MkInit_main',3)
403 
404  convert_qtrc = .true.
405 
406  select case(mkinit_type)
407  case(i_planestate)
408  call mkinit_planestate
409  case(i_tracerbubble)
410  call mkinit_tracerbubble
411  case(i_coldbubble)
412  call mkinit_coldbubble
413  case(i_lambwave)
414  call mkinit_lambwave
415  case(i_gravitywave)
416  call mkinit_gravitywave
417  case(i_khwave)
418  call mkinit_khwave
419  case(i_turbulence)
420  call mkinit_turbulence
421  case(i_mountainwave)
422  call mkinit_mountainwave
423  case(i_warmbubble)
424  call mkinit_warmbubble
425  case(i_supercell)
426  call mkinit_supercell
427  case(i_squallline)
428  call mkinit_squallline
429  case(i_wk1982)
430  call mkinit_wk1982
431  case(i_dycoms2_rf01)
432  call mkinit_dycoms2_rf01
433  case(i_dycoms2_rf02)
434  call mkinit_dycoms2_rf02
435  case(i_rico)
436  call mkinit_rico
437  case(i_bomex)
438  call mkinit_bomex
439  case(i_oceancouple)
440  call mkinit_planestate
441  call mkinit_oceancouple
442  case(i_landcouple)
443  call mkinit_planestate
444  call mkinit_landcouple
445  case(i_urbancouple)
446  call mkinit_planestate
447  call mkinit_urbancouple
448  case(i_triplecouple)
449  call mkinit_planestate
450  call mkinit_oceancouple
451  call mkinit_landcouple
452  call mkinit_urbancouple
453  case(i_bubblecouple)
454  call mkinit_planestate
455  call mkinit_warmbubble
456  call mkinit_oceancouple
457  call mkinit_landcouple
458  call mkinit_urbancouple
459  case(i_seabreeze)
460  call mkinit_planestate
461  call mkinit_seabreeze
462  case(i_heatisland)
463  call mkinit_planestate
464  call mkinit_heatisland
465  case(i_dycoms2_rf02_dns)
466  call mkinit_dycoms2_rf02_dns
467  case(i_real)
468  call mkinit_real
469  convert_qtrc = .false.
470  case(i_grayzone)
471  call mkinit_grayzone
472  case(i_boxaero)
473  call mkinit_boxaero
474  case(i_warmbubbleaero)
475  call mkinit_warmbubbleaero
476  case(i_cavityflow)
477  call mkinit_cavityflow
478  case(i_barocwave)
479  call mkinit_barocwave
480  case default
481  log_error("MKINIT",*) 'Unsupported TYPE:', mkinit_type
482  call prc_abort
483  endselect
484 
485  call tke_setup
486 
487  call aerosol_setup
488 
489  call sbmaero_setup( convert_qtrc ) ! [INOUT]
490 
491  if ( qa_mp > 0 .AND. convert_qtrc ) then
492 !OCL XFILL
493  qhyd(:,:,:,i_hc) = qc(:,:,:)
494 !OCL XFILL
495  qnum(:,:,:,i_hc) = nc(:,:,:)
496  call atmos_phy_mp_driver_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
497  qv(:,:,:), qhyd(:,:,:,:), & ! [IN]
498  qtrc(:,:,:,qs_mp:qe_mp), & ! [OUT]
499  qnum=qnum(:,:,:,:) ) ! [IN]
500  end if
501 
502  call prof_rapend ('_MkInit_main',3)
503 
504  log_progress(*) 'end making initial data'
505 
506  output = .true.
507 
508  endif
509 
510  return
511  end subroutine mkinit
512 
513  !-----------------------------------------------------------------------------
515  subroutine bubble_setup
516  use scale_const, only: &
518  implicit none
519 
520  ! Bubble
521  logical :: BBL_eachnode = .false. ! Arrange bubble at each node? [kg/kg]
522  real(RP) :: BBL_CZ = 2.e3_rp ! center location [m]: z
523  real(RP) :: BBL_CX = 2.e3_rp ! center location [m]: x
524  real(RP) :: BBL_CY = 2.e3_rp ! center location [m]: y
525  real(RP) :: BBL_RZ = 0.0_rp ! bubble radius [m]: z
526  real(RP) :: BBL_RX = 0.0_rp ! bubble radius [m]: x
527  real(RP) :: BBL_RY = 0.0_rp ! bubble radius [m]: y
528 
529  namelist / param_bubble / &
530  bbl_eachnode, &
531  bbl_cz, &
532  bbl_cx, &
533  bbl_cy, &
534  bbl_rz, &
535  bbl_rx, &
536  bbl_ry
537 
538  real(RP) :: CZ_offset
539  real(RP) :: CX_offset
540  real(RP) :: CY_offset
541  real(RP) :: distx, disty, distz
542 
543  real(RP) :: Domain_RX, Domain_RY
544 
545  integer :: ierr
546  integer :: k, i, j
547  !---------------------------------------------------------------------------
548 
549  log_newline
550  log_info("BUBBLE_setup",*) 'Setup'
551 
552  !--- read namelist
553  rewind(io_fid_conf)
554  read(io_fid_conf,nml=param_bubble,iostat=ierr)
555  if( ierr < 0 ) then !--- missing
556  log_info("BUBBLE_setup",*) 'Not found namelist. Default used.'
557  elseif( ierr > 0 ) then !--- fatal error
558  log_error("BUBBLE_setup",*) 'Not appropriate names in namelist PARAM_BUBBLE. Check!'
559  call prc_abort
560  endif
561  log_nml(param_bubble)
562 
563  if ( abs(bbl_rz*bbl_rx*bbl_ry) <= 0.0_rp ) then
564  log_info("BUBBLE_setup",*) 'no bubble'
565  bubble(:,:,:) = 0.0_rp
566  else
567 
568  bubble(:,:,:) = const_undef8
569 
570  if ( bbl_eachnode ) then
571  cz_offset = cz(ks)
572  cx_offset = cx(is)
573  cy_offset = cy(js)
574  domain_rx = fx(ie) - fx(is-1)
575  domain_ry = fy(je) - fy(js-1)
576  else
577  cz_offset = 0.0_rp
578  cx_offset = 0.0_rp
579  cy_offset = 0.0_rp
580  domain_rx = fxg(iag-ihalo) - fxg(ihalo)
581  domain_ry = fyg(jag-jhalo) - fyg(jhalo)
582  endif
583 
584  ! make bubble coefficient
585  do j = 1, ja
586  do i = 1, ia
587  do k = ks, ke
588 
589  distz = ( (cz(k)-cz_offset-bbl_cz)/bbl_rz )**2
590 
591  distx = min( ( (cx(i)-cx_offset-bbl_cx )/bbl_rx )**2, &
592  ( (cx(i)-cx_offset-bbl_cx-domain_rx)/bbl_rx )**2, &
593  ( (cx(i)-cx_offset-bbl_cx+domain_rx)/bbl_rx )**2 )
594 
595  disty = min( ( (cy(j)-cy_offset-bbl_cy )/bbl_ry )**2, &
596  ( (cy(j)-cy_offset-bbl_cy-domain_ry)/bbl_ry )**2, &
597  ( (cy(j)-cy_offset-bbl_cy+domain_ry)/bbl_ry )**2 )
598 
599  bubble(k,i,j) = cos( 0.5_rp*pi*sqrt( min(distz+distx+disty,1.0_rp) ) )**2
600 
601  enddo
602  enddo
603  enddo
604  endif
605 
606  return
607  end subroutine bubble_setup
608 
609  !-----------------------------------------------------------------------------
611  subroutine rect_setup
612  use scale_const, only: &
614  implicit none
615 
616  ! Bubble
617  logical :: RCT_eachnode = .false. ! Arrange rectangle at each node? [kg/kg]
618  real(RP) :: RCT_CZ = 2.e3_rp ! center location [m]: z
619  real(RP) :: RCT_CX = 2.e3_rp ! center location [m]: x
620  real(RP) :: RCT_CY = 2.e3_rp ! center location [m]: y
621  real(RP) :: RCT_RZ = 2.e3_rp ! rectangle z width [m]: z
622  real(RP) :: RCT_RX = 2.e3_rp ! rectangle x width [m]: x
623  real(RP) :: RCT_RY = 2.e3_rp ! rectangle y width [m]: y
624 
625  namelist / param_rect / &
626  rct_eachnode, &
627  rct_cz, &
628  rct_cx, &
629  rct_cy, &
630  rct_rz, &
631  rct_rx, &
632  rct_ry
633 
634  real(RP) :: CZ_offset
635  real(RP) :: CX_offset
636  real(RP) :: CY_offset
637  real(RP) :: dist
638 
639  integer :: ierr
640  integer :: k, i, j
641  !---------------------------------------------------------------------------
642 
643  log_newline
644  log_info("RECT_setup",*) 'Setup'
645 
646  !--- read namelist
647  rewind(io_fid_conf)
648  read(io_fid_conf,nml=param_rect,iostat=ierr)
649  if( ierr < 0 ) then !--- missing
650  log_error("RECT_setup",*) 'Not found namelist. Check!'
651  call prc_abort
652  elseif( ierr > 0 ) then !--- fatal error
653  log_error("RECT_setup",*) 'Not appropriate names in namelist PARAM_RECT. Check!'
654  call prc_abort
655  endif
656  log_nml(param_rect)
657 
658  rect(:,:,:) = const_undef8
659 
660  if ( rct_eachnode ) then
661  cz_offset = cz(ks)
662  cx_offset = cx(is)
663  cy_offset = cy(js)
664  else
665  cz_offset = 0.0_rp
666  cx_offset = 0.0_rp
667  cy_offset = 0.0_rp
668  endif
669 
670  do j = 1, ja
671  do i = 1, ia
672  do k = ks, ke
673 
674  ! make tracer rectangle
675  dist = 2.0_rp * max( &
676  abs(cz(k) - cz_offset - rct_cz)/rct_rz, &
677  abs(cx(i) - cx_offset - rct_cx)/rct_rx, &
678  abs(cy(j) - cy_offset - rct_cy)/rct_ry &
679  & )
680  if ( dist <= 1.0_rp ) then
681  rect(k,i,j) = 1.0_rp
682  else
683  rect(k,i,j) = 0.0_rp
684  end if
685  enddo
686  enddo
687  enddo
688 
689  return
690  end subroutine rect_setup
691 
692  !-----------------------------------------------------------------------------
694  subroutine aerosol_setup
695  use mod_atmos_admin, only: &
697  use scale_atmos_phy_ae_kajino13, only: &
699  use mod_atmos_phy_ae_vars, only: &
700  qa_ae, &
701  qs_ae, &
702  qe_ae
703  implicit none
704 
705  real(RP), parameter :: d_min_def = 1.e-9_rp ! default lower bound of 1st size bin
706  real(RP), parameter :: d_max_def = 1.e-5_rp ! upper bound of last size bin
707  integer, parameter :: n_kap_def = 1 ! number of kappa bins
708  real(RP), parameter :: k_min_def = 0.e0_rp ! lower bound of 1st kappa bin
709  real(RP), parameter :: k_max_def = 1.e0_rp ! upper bound of last kappa bin
710 
711  real(RP) :: m0_init = 0.0_rp ! initial total num. conc. of modes (Atk,Acm,Cor) [#/m3]
712  real(RP) :: dg_init = 80.e-9_rp ! initial number equivalen diameters of modes [m]
713  real(RP) :: sg_init = 1.6_rp ! initial standard deviation [-]
714 
715  real(RP) :: d_min_inp(3) = d_min_def
716  real(RP) :: d_max_inp(3) = d_max_def
717  real(RP) :: k_min_inp(3) = k_min_def
718  real(RP) :: k_max_inp(3) = k_max_def
719  integer :: n_kap_inp(3) = n_kap_def
720 
721  namelist / param_aero / &
722  m0_init, &
723  dg_init, &
724  sg_init, &
725  d_min_inp, &
726  d_max_inp, &
727  k_min_inp, &
728  k_max_inp, &
729  n_kap_inp
730 
731  integer :: ierr
732  !---------------------------------------------------------------------------
733 
734  if ( atmos_phy_ae_type /= 'KAJINO13' ) return
735 
736  log_newline
737  log_info("AEROSOL_setup",*) 'Setup'
738 
739  !--- read namelist
740  rewind(io_fid_conf)
741  read(io_fid_conf,nml=param_aero,iostat=ierr)
742  if( ierr < 0 ) then !--- missing
743  log_info("AEROSOL_setup",*) 'Not found namelist. Default used!'
744  elseif( ierr > 0 ) then !--- fatal error
745  log_error("AEROSOL_setup",*) 'Not appropriate names in namelist PARAM_AERO. Check!'
746  call prc_abort
747  endif
748  log_nml(param_aero)
749 
750  qdry(:,:,:) = 1.0_rp - qv(:,:,:) - qc(:,:,:)
751  call atmos_phy_ae_kajino13_mkinit( ka, ks, ke, ia, is, ie, ja, js, je, & ! (in)
752  qa_ae, & ! (in)
753  dens(:,:,:), & ! (in)
754  temp(:,:,:), & ! (in)
755  pres(:,:,:), & ! (in)
756  qdry(:,:,:), & ! (in)
757  qv(:,:,:), & ! (in)
758  m0_init, & ! (in)
759  dg_init, & ! (in)
760  sg_init, & ! (in)
761  d_min_inp(:), & ! (in)
762  d_max_inp(:), & ! (in)
763  k_min_inp(:), & ! (in)
764  k_max_inp(:), & ! (in)
765  n_kap_inp(:), & ! (in)
766  qtrc(:,:,:,qs_ae:qe_ae), & ! (out)
767  ccn(:,:,:) ) ! (out)
768 
769  return
770  end subroutine aerosol_setup
771 
772  !-----------------------------------------------------------------------------
774  subroutine sbmaero_setup( convert_qtrc )
775  use scale_atmos_hydrometeor, only: &
776  i_qv, &
777  qhs, &
778  qhe
779  use mod_atmos_admin, only: &
781  use scale_atmos_phy_mp_suzuki10, only: &
782  nccn, nbin
783  implicit none
784 
785  logical, intent(inout) :: convert_qtrc
786 
787  real(RP), allocatable :: xabnd(:), xactr(:)
788 
789  integer :: iq, i, j, k
790  !---------------------------------------------------------------------------
791 
792  if ( atmos_phy_mp_type /= 'SUZUKI10' ) return
793 
794  !--- Super saturated air at initial
795  do j = jsb, jeb
796  do i = isb, ieb
797  do k = ks, ke
798  qtrc(k,i,j,i_qv) = qv(k,i,j) + qc(k,i,j)
799  end do
800  end do
801  end do
802 
803  !-- Aerosol distribution
804  if ( nccn /= 0 ) then
805  do iq = 1, nccn
806  do j = jsb, jeb
807  do i = isb, ieb
808  do k = ks, ke
809  qtrc(k,i,j,qhe+iq) = gan(iq) / dens(k,i,j) ! [note] gan is never set.
810  enddo
811  enddo
812  enddo
813  enddo
814 
815  deallocate( xactr )
816  deallocate( xabnd )
817  endif
818 
819  convert_qtrc = .false.
820 
821  return
822  end subroutine sbmaero_setup
823 
824  !-----------------------------------------------------------------------------
825  function faero( f0,r0,x,alpha,rhoa )
826  use scale_const, only: &
827  pi => const_pi
828  implicit none
829 
830  real(RP), intent(in) :: x, f0, r0, alpha, rhoa
831  real(RP) :: faero
832  real(RP) :: rad
833  !---------------------------------------------------------------------------
834 
835  rad = ( exp(x) * 3.0_rp / 4.0_rp / pi / rhoa )**(1.0_rp/3.0_rp)
836 
837  faero = f0 * (rad/r0)**(-alpha)
838 
839  return
840  end function faero
841 
842  !-----------------------------------------------------------------------------
844  subroutine flux_setup
846  sflx_rain => atmos_phy_mp_sflx_rain, &
847  sflx_snow => atmos_phy_mp_sflx_snow
848  use mod_atmos_phy_rd_vars, only: &
849  sflx_lw_up => atmos_phy_rd_sflx_lw_up, &
850  sflx_lw_dn => atmos_phy_rd_sflx_lw_dn, &
851  sflx_sw_up => atmos_phy_rd_sflx_sw_up, &
852  sflx_sw_dn => atmos_phy_rd_sflx_sw_dn, &
853  toaflx_lw_up => atmos_phy_rd_toaflx_lw_up, &
854  toaflx_lw_dn => atmos_phy_rd_toaflx_lw_dn, &
855  toaflx_sw_up => atmos_phy_rd_toaflx_sw_up, &
856  toaflx_sw_dn => atmos_phy_rd_toaflx_sw_dn, &
857  sflx_rad_dn => atmos_phy_rd_sflx_down
858  implicit none
859 
860  ! Flux from Atmosphere
861  real(RP) :: FLX_rain = 0.0_rp ! surface rain flux [kg/m2/s]
862  real(RP) :: FLX_snow = 0.0_rp ! surface snow flux [kg/m2/s]
863  real(RP) :: FLX_IR_dn = 0.0_rp ! surface downwad radiation flux [J/m2/s]
864  real(RP) :: FLX_NIR_dn = 0.0_rp ! surface downwad radiation flux [J/m2/s]
865  real(RP) :: FLX_VIS_dn = 0.0_rp ! surface downwad radiation flux [J/m2/s]
866 
867  namelist / param_mkinit_flux / &
868  flx_rain, &
869  flx_snow, &
870  flx_ir_dn, &
871  flx_nir_dn, &
872  flx_vis_dn
873 
874  integer :: i, j
875  integer :: ierr
876  !---------------------------------------------------------------------------
877 
878  !--- read namelist
879  rewind(io_fid_conf)
880  read(io_fid_conf,nml=param_mkinit_flux,iostat=ierr)
881  if( ierr < 0 ) then !--- missing
882  log_info("flux_setup",*) 'Not found namelist. Default used.'
883  elseif( ierr > 0 ) then !--- fatal error
884  log_error("flux_setup",*) 'Not appropriate names in namelist PARAM_MKINIT_FLUX. Check!'
885  call prc_abort
886  endif
887  log_nml(param_mkinit_flux)
888 
889  do j = jsb, jeb
890  do i = isb, ieb
891  sflx_rain(i,j) = flx_rain
892  sflx_snow(i,j) = flx_snow
893 
894  sflx_lw_up(i,j) = 0.0_rp
895  sflx_lw_dn(i,j) = flx_ir_dn
896  sflx_sw_up(i,j) = 0.0_rp
897  sflx_sw_dn(i,j) = flx_nir_dn + flx_vis_dn
898 
899  toaflx_lw_up(i,j) = 0.0_rp
900  toaflx_lw_dn(i,j) = 0.0_rp
901  toaflx_sw_up(i,j) = 0.0_rp
902  toaflx_sw_dn(i,j) = 0.0_rp
903 
904  sflx_rad_dn(i,j,i_r_direct ,i_r_ir) = 0.0_rp
905  sflx_rad_dn(i,j,i_r_diffuse,i_r_ir) = flx_ir_dn
906  sflx_rad_dn(i,j,i_r_direct ,i_r_nir) = flx_nir_dn
907  sflx_rad_dn(i,j,i_r_diffuse,i_r_nir) = 0.0_rp
908  sflx_rad_dn(i,j,i_r_direct ,i_r_vis) = flx_vis_dn
909  sflx_rad_dn(i,j,i_r_diffuse,i_r_vis) = 0.0_rp
910  enddo
911  enddo
912 
913  return
914  end subroutine flux_setup
915 
916  !-----------------------------------------------------------------------------
918  subroutine land_setup
919  use mod_land_vars, only: &
920  land_temp, &
921  land_water, &
922  land_sfc_temp, &
924  implicit none
925 
926  real(RP) :: LND_TEMP ! land soil temperature [K]
927  real(RP) :: LND_WATER = 0.15_rp ! land soil moisture [m3/m3]
928  real(RP) :: SFC_TEMP ! land skin temperature [K]
929  real(RP) :: SFC_albedo_LW = 0.01_rp ! land surface albedo for LW (0-1)
930  real(RP) :: SFC_albedo_SW = 0.20_rp ! land surface albedo for SW (0-1)
931 
932  namelist / param_mkinit_land / &
933  lnd_temp, &
934  lnd_water, &
935  sfc_temp, &
936  sfc_albedo_lw, &
937  sfc_albedo_sw
938 
939  integer :: ierr
940  !---------------------------------------------------------------------------
941 
942  lnd_temp = thetastd
943  sfc_temp = thetastd
944 
945  !--- read namelist
946  rewind(io_fid_conf)
947  read(io_fid_conf,nml=param_mkinit_land,iostat=ierr)
948  if( ierr < 0 ) then !--- missing
949  log_info("land_setup",*) 'Not found namelist. Default used.'
950  elseif( ierr > 0 ) then !--- fatal error
951  log_error("land_setup",*) 'Not appropriate names in namelist PARAM_MKINIT_LAND. Check!'
952  call prc_abort
953  endif
954  log_nml(param_mkinit_land)
955 
956  land_temp(:,:,:) = lnd_temp
957  land_water(:,:,:) = lnd_water
958 
959  land_sfc_temp(:,:) = sfc_temp
960  land_sfc_albedo(:,:,:,i_r_ir) = sfc_albedo_lw
961  land_sfc_albedo(:,:,:,i_r_nir) = sfc_albedo_sw
962  land_sfc_albedo(:,:,:,i_r_vis) = sfc_albedo_sw
963 
964  return
965  end subroutine land_setup
966 
967  !-----------------------------------------------------------------------------
969  subroutine ocean_setup
970  use mod_ocean_vars, only: &
971  ocean_temp, &
972  ocean_salt, &
973  ocean_uvel, &
974  ocean_vvel, &
975  ocean_ocn_z0m, &
976  ocean_ice_temp, &
977  ocean_ice_mass, &
978  ocean_sfc_temp, &
980  ocean_sfc_z0m, &
981  ocean_sfc_z0h, &
983  implicit none
984 
985  real(RP) :: OCN_TEMP ! ocean temperature [K]
986  real(RP) :: OCN_SALT = 0.0_rp ! ocean salinity [psu]
987  real(RP) :: OCN_UVEL = 0.0_rp ! ocean u-velocity [m/s]
988  real(RP) :: OCN_VVEL = 0.0_rp ! ocean v-velocity [m/s]
989  real(RP) :: ICE_TEMP ! ocean temperature [K]
990  real(RP) :: ICE_MASS = 0.0_rp ! ocean temperature [K]
991  real(RP) :: SFC_TEMP ! ocean skin temperature [K]
992  real(RP) :: SFC_albedo_LW = 0.04_rp ! ocean surface albedo for LW (0-1)
993  real(RP) :: SFC_albedo_SW = 0.05_rp ! ocean surface albedo for SW (0-1)
994  real(RP) :: SFC_Z0M = 1.e-4_rp ! ocean surface roughness length (momentum) [m]
995  real(RP) :: SFC_Z0H = 1.e-4_rp ! ocean surface roughness length (heat) [m]
996  real(RP) :: SFC_Z0E = 1.e-4_rp ! ocean surface roughness length (vapor) [m]
997 
998  namelist / param_mkinit_ocean / &
999  ocn_temp, &
1000  ocn_salt, &
1001  ocn_uvel, &
1002  ocn_vvel, &
1003  ice_temp, &
1004  ice_mass, &
1005  sfc_temp, &
1006  sfc_albedo_lw, &
1007  sfc_albedo_sw, &
1008  sfc_z0m, &
1009  sfc_z0h, &
1010  sfc_z0e
1011 
1012  integer :: ierr
1013  !---------------------------------------------------------------------------
1014 
1015  ocn_temp = thetastd
1016  ice_temp = 271.35_rp ! freezing point of the ocean
1017  sfc_temp = thetastd
1018 
1019  !--- read namelist
1020  rewind(io_fid_conf)
1021  read(io_fid_conf,nml=param_mkinit_ocean,iostat=ierr)
1022  if( ierr < 0 ) then !--- missing
1023  log_info("ocean_setup",*) 'Not found namelist. Default used.'
1024  elseif( ierr > 0 ) then !--- fatal error
1025  log_error("ocean_setup",*) 'Not appropriate names in namelist PARAM_MKINIT_OCEAN. Check!'
1026  call prc_abort
1027  endif
1028  log_nml(param_mkinit_ocean)
1029 
1030  ocean_temp(:,:,:) = ocn_temp
1031  ocean_salt(:,:,:) = ocn_salt
1032  ocean_uvel(:,:,:) = ocn_uvel
1033  ocean_vvel(:,:,:) = ocn_vvel
1034  ocean_ocn_z0m(:,:) = sfc_z0m
1035  ocean_ice_temp(:,:) = ice_temp
1036  ocean_ice_mass(:,:) = ice_mass
1037 
1038  ocean_sfc_temp(:,:) = sfc_temp
1039  ocean_sfc_albedo(:,:,:,i_r_ir) = sfc_albedo_lw
1040  ocean_sfc_albedo(:,:,:,i_r_nir) = sfc_albedo_sw
1041  ocean_sfc_albedo(:,:,:,i_r_vis) = sfc_albedo_sw
1042  ocean_sfc_z0m(:,:) = sfc_z0m
1043  ocean_sfc_z0h(:,:) = sfc_z0h
1044  ocean_sfc_z0e(:,:) = sfc_z0e
1045 
1046  return
1047  end subroutine ocean_setup
1048 
1049  !-----------------------------------------------------------------------------
1051  subroutine urban_setup
1052  use mod_urban_vars, only: &
1053  urban_tr, &
1054  urban_tb, &
1055  urban_tg, &
1056  urban_tc, &
1057  urban_qc, &
1058  urban_uc, &
1059  urban_trl, &
1060  urban_tbl, &
1061  urban_tgl, &
1062  urban_rainr, &
1063  urban_rainb, &
1064  urban_raing, &
1065  urban_roff, &
1066  urban_sfc_temp, &
1068  implicit none
1069 
1070  real(RP) :: URB_ROOF_TEMP ! Surface temperature of roof [K]
1071  real(RP) :: URB_BLDG_TEMP ! Surface temperature of building [K]
1072  real(RP) :: URB_GRND_TEMP ! Surface temperature of ground [K]
1073  real(RP) :: URB_CNPY_TEMP ! Diagnostic canopy air temperature [K]
1074  real(RP) :: URB_CNPY_HMDT = 0.0_rp ! Diagnostic canopy humidity [kg/kg]
1075  real(RP) :: URB_CNPY_WIND = 0.0_rp ! Diagnostic canopy wind [m/s]
1076  real(RP) :: URB_ROOF_LAYER_TEMP ! temperature in layer of roof [K]
1077  real(RP) :: URB_BLDG_LAYER_TEMP ! temperature in layer of building [K]
1078  real(RP) :: URB_GRND_LAYER_TEMP ! temperature in layer of ground [K]
1079  real(RP) :: URB_ROOF_RAIN = 0.0_rp ! temperature in layer of roof [K]
1080  real(RP) :: URB_BLDG_RAIN = 0.0_rp ! temperature in layer of building [K]
1081  real(RP) :: URB_GRND_RAIN = 0.0_rp ! temperature in layer of ground [K]
1082  real(RP) :: URB_RUNOFF = 0.0_rp ! temperature in layer of ground [K]
1083  real(RP) :: URB_SFC_TEMP ! Grid average of surface temperature [K]
1084  real(RP) :: URB_ALB_LW = 0.10_rp ! Grid average of surface albedo for LW (0-1)
1085  real(RP) :: URB_ALB_SW = 0.20_rp ! Grid average of surface albedo for SW (0-1)
1086 
1087  namelist / param_mkinit_urban / &
1088  urb_roof_temp, &
1089  urb_bldg_temp, &
1090  urb_grnd_temp, &
1091  urb_cnpy_temp, &
1092  urb_cnpy_hmdt, &
1093  urb_cnpy_wind, &
1094  urb_roof_layer_temp, &
1095  urb_bldg_layer_temp, &
1096  urb_grnd_layer_temp, &
1097  urb_roof_rain, &
1098  urb_bldg_rain, &
1099  urb_grnd_rain, &
1100  urb_runoff, &
1101  urb_sfc_temp, &
1102  urb_alb_lw, &
1103  urb_alb_sw
1104 
1105  integer :: ierr
1106  !---------------------------------------------------------------------------
1107 
1108  urb_roof_temp = thetastd
1109  urb_bldg_temp = thetastd
1110  urb_grnd_temp = thetastd
1111  urb_cnpy_temp = thetastd
1112  urb_roof_layer_temp = thetastd
1113  urb_bldg_layer_temp = thetastd
1114  urb_grnd_layer_temp = thetastd
1115  urb_sfc_temp = thetastd
1116 
1117  !--- read namelist
1118  rewind(io_fid_conf)
1119  read(io_fid_conf,nml=param_mkinit_urban,iostat=ierr)
1120  if( ierr < 0 ) then !--- missing
1121  log_info("urban_setup",*) 'Not found namelist. Default used.'
1122  elseif( ierr > 0 ) then !--- fatal error
1123  log_error("urban_setup",*) 'Not appropriate names in namelist PARAM_MKINIT_URBAN. Check!'
1124  call prc_abort
1125  endif
1126  log_nml(param_mkinit_urban)
1127 
1128  urban_trl(:,:,:) = urb_roof_layer_temp
1129  urban_tbl(:,:,:) = urb_bldg_layer_temp
1130  urban_tgl(:,:,:) = urb_grnd_layer_temp
1131 
1132  urban_tr(:,:) = urb_roof_temp
1133  urban_tb(:,:) = urb_bldg_temp
1134  urban_tg(:,:) = urb_grnd_temp
1135  urban_tc(:,:) = urb_cnpy_temp
1136  urban_qc(:,:) = urb_cnpy_hmdt
1137  urban_uc(:,:) = urb_cnpy_wind
1138  urban_rainr(:,:) = urb_roof_rain
1139  urban_rainb(:,:) = urb_bldg_rain
1140  urban_raing(:,:) = urb_grnd_rain
1141  urban_roff(:,:) = urb_runoff
1142  urban_sfc_temp(:,:) = urb_sfc_temp
1143  urban_sfc_albedo(:,:,:,i_r_ir) = urb_alb_lw
1144  urban_sfc_albedo(:,:,:,i_r_nir) = urb_alb_sw
1145  urban_sfc_albedo(:,:,:,i_r_vis) = urb_alb_sw
1146 
1147  return
1148  end subroutine urban_setup
1149 
1150  !-----------------------------------------------------------------------------
1152  subroutine tke_setup
1153  use scale_const, only: &
1154  eps => const_eps
1155  use mod_atmos_phy_tb_vars, only: &
1156  i_tke
1157  use mod_atmos_phy_bl_vars, only: &
1158  qs_bl => qs, &
1159  qe_bl => qe
1160  implicit none
1161 
1162  real(RP) :: TKE_CONST
1163 
1164  namelist / param_mkinit_tke / &
1165  tke_const
1166 
1167  integer :: k, i, j
1168  integer :: ierr
1169  !---------------------------------------------------------------------------
1170 
1171  tke_const = eps
1172 
1173  !--- read namelist
1174  rewind(io_fid_conf)
1175  read(io_fid_conf,nml=param_mkinit_tke,iostat=ierr)
1176  if( ierr < 0 ) then !--- missing
1177  log_info("tke_setup",*) 'Not found namelist. Default used.'
1178  elseif( ierr > 0 ) then !--- fatal error
1179  log_error("tke_setup",*) 'Not appropriate names in namelist PARAM_MKINIT_TKE. Check!'
1180  call prc_abort
1181  endif
1182  log_nml(param_mkinit_tke)
1183 
1184  if ( i_tke > 0 ) then
1185  do j = 1, ja
1186  do i = 1, ia
1187  do k = 1, ka
1188  qtrc(k,i,j,i_tke) = tke_const
1189  enddo
1190  enddo
1191  enddo
1192  end if
1193  if ( qs_bl > 0 ) then
1194  do j = 1, ja
1195  do i = 1, ia
1196  do k = 1, ka
1197  qtrc(k,i,j,qs_bl) = tke_const
1198  qtrc(k,i,j,qs_bl+1:qe_bl) = 0.0_rp
1199  enddo
1200  enddo
1201  enddo
1202  end if
1203 
1204  return
1205  end subroutine tke_setup
1206 
1207  !-----------------------------------------------------------------------------
1209  subroutine read_sounding( &
1210  DENS, VELX, VELY, POTT, QV )
1213  implicit none
1214 
1215  real(RP), intent(out) :: DENS(ka)
1216  real(RP), intent(out) :: VELX(ka)
1217  real(RP), intent(out) :: VELY(ka)
1218  real(RP), intent(out) :: POTT(ka)
1219  real(RP), intent(out) :: QV (ka)
1220 
1221  real(RP) :: TEMP(ka)
1222  real(RP) :: PRES(ka)
1223  real(RP) :: QC (ka)
1224 
1225  character(len=H_LONG) :: ENV_IN_SOUNDING_file = ''
1226 
1227  integer, parameter :: EXP_klim = 100
1228  integer :: EXP_kmax
1229 
1230  real(RP) :: SFC_THETA ! surface potential temperature [K]
1231  real(RP) :: SFC_PRES ! surface pressure [hPa]
1232  real(RP) :: SFC_QV ! surface watervapor [g/kg]
1233 
1234  real(RP) :: EXP_z (exp_klim+1) ! height [m]
1235  real(RP) :: EXP_pott(exp_klim+1) ! potential temperature [K]
1236  real(RP) :: EXP_qv (exp_klim+1) ! water vapor [g/kg]
1237  real(RP) :: EXP_u (exp_klim+1) ! velocity u [m/s]
1238  real(RP) :: EXP_v (exp_klim+1) ! velocity v [m/s]
1239 
1240  real(RP) :: fact1, fact2
1241  integer :: k, kref
1242  integer :: fid
1243  integer :: ierr
1244 
1245  namelist / param_mkinit_sounding / &
1246  env_in_sounding_file
1247 
1248  !--- read namelist
1249  rewind(io_fid_conf)
1250  read(io_fid_conf,nml=param_mkinit_sounding,iostat=ierr)
1251 
1252  if( ierr < 0 ) then !--- missing
1253  log_info("read_sounding",*) 'Not found namelist. Default used.'
1254  elseif( ierr > 0 ) then !--- fatal error
1255  log_error("read_sounding",*) 'Not appropriate names in namelist PARAM_MKINIT_SOUNDING. Check!'
1256  call prc_abort
1257  endif
1258  log_nml(param_mkinit_sounding)
1259 
1260  !--- prepare sounding profile
1261  log_info("read_sounding",*) 'Input sounding file:', trim(env_in_sounding_file)
1262  fid = io_get_available_fid()
1263  open( fid, &
1264  file = trim(env_in_sounding_file), &
1265  form = 'formatted', &
1266  status = 'old', &
1267  iostat = ierr )
1268 
1269  if ( ierr /= 0 ) then
1270  log_error("read_sounding",*) '[mod_mkinit/read_sounding] Input file not found!'
1271  endif
1272 
1273  !--- read sounding file till end
1274  read(fid,*) sfc_pres, sfc_theta, sfc_qv
1275 
1276  log_info("read_sounding",*) '+ Surface pressure [hPa]', sfc_pres
1277  log_info("read_sounding",*) '+ Surface pot. temp [K]', sfc_theta
1278  log_info("read_sounding",*) '+ Surface water vapor [g/kg]', sfc_qv
1279 
1280  do k = 2, exp_klim
1281  read(fid,*,iostat=ierr) exp_z(k), exp_pott(k), exp_qv(k), exp_u(k), exp_v(k)
1282  if ( ierr /= 0 ) exit
1283  enddo
1284 
1285  exp_kmax = k - 1
1286  close(fid)
1287 
1288  ! Boundary
1289  exp_z(1) = 0.0_rp
1290  exp_pott(1) = sfc_theta
1291  exp_qv(1) = sfc_qv
1292  exp_u(1) = exp_u(2)
1293  exp_v(1) = exp_v(2)
1294  exp_z(exp_kmax+1) = 100.e3_rp
1295  exp_pott(exp_kmax+1) = exp_pott(exp_kmax)
1296  exp_qv(exp_kmax+1) = exp_qv(exp_kmax)
1297  exp_u(exp_kmax+1) = exp_u(exp_kmax)
1298  exp_v(exp_kmax+1) = exp_v(exp_kmax)
1299 
1300  do k = 1, exp_kmax+1
1301  exp_qv(k) = exp_qv(k) * 1.e-3_rp ! [g/kg]->[kg/kg]
1302  enddo
1303 
1304  ! calc in dry condition
1305  pres_sfc(:,:) = sfc_pres * 1.e2_rp ! [hPa]->[Pa]
1306  pott_sfc(:,:) = sfc_theta
1307  if ( .not. atmos_hydrometeor_dry ) then
1308  qv_sfc(:,:) = sfc_qv * 1.e-3_rp ! [g/kg]->[kg/kg]
1309  end if
1310 
1311  !--- linear interpolate to model grid
1312  do k = ks, ke
1313  do kref = 2, exp_kmax+1
1314  if ( cz(k) > exp_z(kref-1) &
1315  .AND. cz(k) <= exp_z(kref ) ) then
1316 
1317  fact1 = ( exp_z(kref) - cz(k) ) / ( exp_z(kref)-exp_z(kref-1) )
1318  fact2 = ( cz(k) - exp_z(kref-1) ) / ( exp_z(kref)-exp_z(kref-1) )
1319 
1320  pott(k) = exp_pott(kref-1) * fact1 &
1321  + exp_pott(kref ) * fact2
1322  qv(k) = exp_qv(kref-1) * fact1 &
1323  + exp_qv(kref ) * fact2
1324  velx(k) = exp_u(kref-1) * fact1 &
1325  + exp_u(kref ) * fact2
1326  vely(k) = exp_v(kref-1) * fact1 &
1327  + exp_v(kref ) * fact2
1328  endif
1329  enddo
1330  enddo
1331  if ( atmos_hydrometeor_dry ) qv(:) = 0.0_rp
1332 
1333  qc(:) = 0.0_rp
1334 
1335  ! make density & pressure profile in moist condition
1336  call hydrostatic_buildrho( ka, ks, ke, &
1337  pott(:), qv(:), qc(:), & ! [IN]
1338  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
1339  cz(:), fz(:), & ! [IN]
1340  dens(:), temp(:), pres(:), temp_sfc(1,1) ) ! [OUT]
1341 
1342  return
1343  end subroutine read_sounding
1344 
1345  !-----------------------------------------------------------------------------
1347  subroutine mkinit_planestate
1348  use scale_atmos_hydrometeor, only: &
1350  implicit none
1351 
1352  ! Surface state
1353  real(RP) :: SFC_THETA ! surface potential temperature [K]
1354  real(RP) :: SFC_PRES ! surface pressure [Pa]
1355  real(RP) :: SFC_RH = 0.0_rp ! surface relative humidity [%]
1356  ! Environment state
1357  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1358  real(RP) :: ENV_TLAPS = 0.0_rp ! Lapse rate of THETA [K/m]
1359  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
1360  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1361  real(RP) :: ENV_RH = 0.0_rp ! relative humidity of environment [%]
1362  ! Disturbance
1363  real(RP) :: RANDOM_THETA = 0.0_rp ! amplitude of random disturbance theta
1364  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
1365  real(RP) :: RANDOM_V = 0.0_rp ! amplitude of random disturbance v
1366  real(RP) :: RANDOM_RH = 0.0_rp ! amplitude of random disturbance RH
1367 
1368  namelist / param_mkinit_planestate / &
1369  sfc_theta, &
1370  sfc_pres, &
1371  sfc_rh, &
1372  env_theta, &
1373  env_tlaps, &
1374  env_u, &
1375  env_v, &
1376  env_rh, &
1377  random_theta, &
1378  random_u, &
1379  random_v, &
1380  random_rh
1381 
1382  integer :: ierr
1383  integer :: k, i, j
1384  !---------------------------------------------------------------------------
1385 
1386  log_newline
1387  log_info("MKINIT_planestate",*) 'Setup initial state'
1388 
1389  sfc_theta = thetastd
1390  sfc_pres = pstd
1391  env_theta = thetastd
1392 
1393  !--- read namelist
1394  rewind(io_fid_conf)
1395  read(io_fid_conf,nml=param_mkinit_planestate,iostat=ierr)
1396 
1397  if( ierr < 0 ) then !--- missing
1398  log_info("MKINIT_planestate",*) 'Not found namelist. Default used.'
1399  elseif( ierr > 0 ) then !--- fatal error
1400  log_error_cont(*) 'Not appropriate names in namelist PARAM_MKINIT_PLANESTATE. Check!'
1401  call prc_abort
1402  endif
1403  log_nml(param_mkinit_planestate)
1404 
1405  ! calc in dry condition
1406  do j = jsb, jeb
1407  do i = isb, ieb
1408  pott_sfc(i,j) = sfc_theta
1409  pres_sfc(i,j) = sfc_pres
1410  enddo
1411  enddo
1412 
1413  if ( env_theta < 0.0_rp ) then ! use isa profile
1414 
1415  call profile_isa( ka, ks, ke, & ! [IN]
1416  ia, isb, ieb, & ! [IN]
1417  ja, jsb, jeb, & ! [IN]
1418  pott_sfc(:,:), & ! [IN]
1419  pres_sfc(:,:), & ! [IN]
1420  real_cz(:,:,:), & ! [IN]
1421  pott(:,:,:) ) ! [OUT]
1422 
1423  else
1424 
1425  do j = jsb, jeb
1426  do i = isb, ieb
1427  do k = ks, ke
1428  pott(k,i,j) = env_theta + env_tlaps * real_cz(k,i,j)
1429  enddo
1430  enddo
1431  enddo
1432 
1433  endif
1434 
1435  ! make density & pressure profile in moist condition
1436  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
1437  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
1438  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
1439  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
1440  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
1441 
1442  if ( .not. atmos_hydrometeor_dry ) then
1443  ! calc QV from RH
1444  call saturation_psat_all( ia, isb, ieb, ja, jsb, jeb, &
1445  temp_sfc(:,:), & ! [IN]
1446  psat_sfc(:,:) ) ! [OUT]
1447  qdry(:,:,:) = 1.0_rp - qv(:,:,:) - qc(:,:,:)
1448  call saturation_pres2qsat_all( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
1449  temp(:,:,:), pres(:,:,:), qdry(:,:,:), & ! [IN]
1450  qsat(:,:,:) ) ! [OUT]
1451 
1452  call random_uniform(rndm) ! make random
1453  do j = jsb, jeb
1454  do i = isb, ieb
1455  qsat_sfc(i,j) = epsvap * psat_sfc(i,j) / ( pres_sfc(i,j) - ( 1.0_rp-epsvap ) * psat_sfc(i,j) )
1456  qv_sfc(i,j) = ( sfc_rh + rndm(ks-1,i,j) * random_rh ) * 1.e-2_rp * qsat_sfc(i,j)
1457 
1458  do k = ks, ke
1459  qv(k,i,j) = ( env_rh + rndm(k,i,j) * random_rh ) * 1.e-2_rp * qsat(k,i,j)
1460  enddo
1461  enddo
1462  enddo
1463  end if
1464 
1465  call random_uniform(rndm) ! make random
1466  do j = jsb, jeb
1467  do i = isb, ieb
1468  pott_sfc(i,j) = pott_sfc(i,j) + rndm(ks-1,i,j) * random_theta
1469 
1470  do k = ks, ke
1471  pott(k,i,j) = pott(k,i,j) + rndm(k,i,j) * random_theta
1472  enddo
1473  enddo
1474  enddo
1475 
1476  ! make density & pressure profile in moist condition
1477  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
1478  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
1479  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
1480  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
1481  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
1482 
1483  call comm_vars8( dens(:,:,:), 1 )
1484  call comm_wait ( dens(:,:,:), 1 )
1485 
1486  call random_uniform(rndm) ! make random
1487  do j = jsb, jeb
1488  do i = isb, min(ieb,ia-1)
1489  do k = ks, ke
1490  momx(k,i,j) = ( env_u + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u ) &
1491  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
1492  enddo
1493  enddo
1494  enddo
1495 
1496  call random_uniform(rndm) ! make random
1497  do j = jsb, min(jeb,ja-1)
1498  do i = isb, ieb
1499  do k = ks, ke
1500  momy(k,i,j) = ( env_v + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_v ) &
1501  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
1502  enddo
1503  enddo
1504  enddo
1505 
1506  do j = jsb, jeb
1507  do i = isb, ieb
1508  do k = ks, ke
1509  momz(k,i,j) = 0.0_rp
1510  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
1511  enddo
1512  enddo
1513  enddo
1514 
1515  call flux_setup
1516 
1517  return
1518  end subroutine mkinit_planestate
1519 
1520  !-----------------------------------------------------------------------------
1522  subroutine mkinit_tracerbubble
1523  implicit none
1524 
1525 #ifndef DRY
1526  ! Surface state
1527  real(RP) :: SFC_THETA ! surface potential temperature [K]
1528  real(RP) :: SFC_PRES ! surface pressure [Pa]
1529  ! Environment state
1530  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1531  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
1532  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1533  ! Bubble
1534  character(len=H_SHORT) :: SHAPE_NC = 'BUBBLE' ! BUBBLE or RECT
1535  real(RP) :: BBL_NC = 1.0_rp ! extremum of NC in bubble [kg/kg]
1536 
1537  namelist / param_mkinit_tracerbubble / &
1538  sfc_theta, &
1539  sfc_pres, &
1540  env_theta, &
1541  env_u, &
1542  env_v, &
1543  shape_nc, &
1544  bbl_nc
1545 
1546  real(RP), pointer :: shapeFac(:,:,:) => null()
1547 
1548  integer :: k, i, j
1549  integer :: ierr
1550  !---------------------------------------------------------------------------
1551 
1552  log_newline
1553  log_info("MKINIT_tracerbubble",*) 'Setup initial state'
1554 
1555  sfc_theta = thetastd
1556  sfc_pres = pstd
1557  env_theta = thetastd
1558 
1559  !--- read namelist
1560  rewind(io_fid_conf)
1561  read(io_fid_conf,nml=param_mkinit_tracerbubble,iostat=ierr)
1562 
1563  if( ierr < 0 ) then !--- missing
1564  log_info("MKINIT_tracerbubble",*) 'Not found namelist. Default used.'
1565  elseif( ierr > 0 ) then !--- fatal error
1566  log_error("MKINIT_tracerbubble",*) 'Not appropriate names in namelist PARAM_MKINIT_TRACERBUBBLE. Check!'
1567  call prc_abort
1568  endif
1569  log_nml(param_mkinit_tracerbubble)
1570 
1571  ! calc in dry condition
1572  pres_sfc(1,1) = sfc_pres
1573  pott_sfc(1,1) = sfc_theta
1574 
1575  do k = ks, ke
1576  pott(k,1,1) = env_theta
1577  enddo
1578 
1579  ! make density & pressure profile in dry condition
1580  call hydrostatic_buildrho( ka, ks, ke, &
1581  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
1582  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
1583  cz(:), fz(:), & ! [IN]
1584  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
1585 
1586  do j = jsb, jeb
1587  do i = isb, ieb
1588  do k = ks, ke
1589  dens(k,i,j) = dens(k,1,1)
1590  momz(k,i,j) = 0.0_rp
1591  momx(k,i,j) = env_u * dens(k,1,1)
1592  momy(k,i,j) = env_v * dens(k,1,1)
1593  rhot(k,i,j) = pott(k,1,1) * dens(k,1,1)
1594  enddo
1595  enddo
1596  enddo
1597 
1598  ! make tracer bubble
1599  select case(shape_nc)
1600  case('BUBBLE')
1601  call bubble_setup
1602  shapefac => bubble
1603  case('RECT')
1604  call rect_setup
1605  shapefac => rect
1606  case default
1607  log_error("MKINIT_tracerbubble",*) 'SHAPE_NC=', trim(shape_nc), ' cannot be used on advect. Check!'
1608  call prc_abort
1609  end select
1610 
1611  do j = jsb, jeb
1612  do i = isb, ieb
1613  do k = ks, ke
1614  nc(k,i,j) = bbl_nc * shapefac(k,i,j)
1615  enddo
1616  enddo
1617  enddo
1618 
1619 #endif
1620 
1621  return
1622  end subroutine mkinit_tracerbubble
1623 
1624  !-----------------------------------------------------------------------------
1634  subroutine mkinit_coldbubble
1635  implicit none
1636 
1637  ! Surface state
1638  real(RP) :: SFC_THETA ! surface potential temperature [K]
1639  real(RP) :: SFC_PRES ! surface pressure [Pa]
1640  ! Environment state
1641  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1642  ! Bubble
1643  real(RP) :: BBL_TEMP = -15.0_rp ! extremum of temperature in bubble [K]
1644 
1645  namelist / param_mkinit_coldbubble / &
1646  sfc_theta, &
1647  sfc_pres, &
1648  env_theta, &
1649  bbl_temp
1650 
1651  real(RP) :: RovCP
1652 
1653  integer :: ierr
1654  integer :: k, i, j
1655  !---------------------------------------------------------------------------
1656 
1657  log_newline
1658  log_info("MKINIT_coldbubble",*) 'Setup initial state'
1659 
1660  sfc_theta = thetastd
1661  sfc_pres = pstd
1662  env_theta = thetastd
1663 
1664  !--- read namelist
1665  rewind(io_fid_conf)
1666  read(io_fid_conf,nml=param_mkinit_coldbubble,iostat=ierr)
1667 
1668  if( ierr < 0 ) then !--- missing
1669  log_info("MKINIT_coldbubble",*) 'Not found namelist. Default used.'
1670  elseif( ierr > 0 ) then !--- fatal error
1671  log_error("MKINIT_coldbubble",*) 'Not appropriate names in namelist PARAM_MKINIT_COLDBUBBLE. Check!'
1672  call prc_abort
1673  endif
1674  log_nml(param_mkinit_coldbubble)
1675 
1676  rovcp = rdry / cpdry
1677 
1678  ! calc in dry condition
1679  pres_sfc(1,1) = sfc_pres
1680  pott_sfc(1,1) = sfc_theta
1681 
1682  do k = ks, ke
1683  pott(k,1,1) = env_theta
1684  enddo
1685 
1686  ! make density & pressure profile in dry condition
1687  call hydrostatic_buildrho( ka, ks, ke, &
1688  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
1689  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
1690  cz(:), fz(:), & ! [IN]
1691  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
1692 
1693  do j = 1, ja
1694  do i = 1, ia
1695  do k = ks, ke
1696  dens(k,i,j) = dens(k,1,1)
1697  momz(k,i,j) = 0.0_rp
1698  momx(k,i,j) = 0.0_rp
1699  momy(k,i,j) = 0.0_rp
1700 
1701  ! make cold bubble
1702  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) &
1703  + bbl_temp * ( p00/pres(k,1,1) )**rovcp * bubble(k,i,j) )
1704  enddo
1705  enddo
1706  enddo
1707 
1708  return
1709  end subroutine mkinit_coldbubble
1710 
1711  !-----------------------------------------------------------------------------
1713  subroutine mkinit_lambwave
1714  implicit none
1715 
1716  ! Surface state
1717  real(RP) :: SFC_PRES ! surface pressure [Pa]
1718  ! Environment state
1719  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
1720  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1721  real(RP) :: ENV_TEMP = 300.0_rp ! temperature of environment [K]
1722  ! Bubble
1723  real(RP) :: BBL_PRES = 100._rp ! extremum of pressure in bubble [Pa]
1724 
1725  namelist / param_mkinit_lambwave / &
1726  sfc_pres, &
1727  env_u, &
1728  env_v, &
1729  env_temp, &
1730  bbl_pres
1731 
1732  real(RP) :: RovCP
1733 
1734  integer :: ierr
1735  integer :: k, i, j
1736  !---------------------------------------------------------------------------
1737 
1738  log_newline
1739  log_info("MKINIT_lambwave",*) 'Setup initial state'
1740 
1741  sfc_pres = pstd
1742 
1743  !--- read namelist
1744  rewind(io_fid_conf)
1745  read(io_fid_conf,nml=param_mkinit_lambwave,iostat=ierr)
1746 
1747  if( ierr < 0 ) then !--- missing
1748  log_info("MKINIT_lambwave",*) 'Not found namelist. Default used.'
1749  elseif( ierr > 0 ) then !--- fatal error
1750  log_error("MKINIT_lambwave",*) 'Not appropriate names in namelist PARAM_MKINIT_LAMBWAVE. Check!'
1751  call prc_abort
1752  endif
1753  log_nml(param_mkinit_lambwave)
1754 
1755  rovcp = rdry / cpdry
1756 
1757  do j = jsb, jeb
1758  do i = isb, ieb
1759  do k = ks, ke
1760  dens(k,i,j) = sfc_pres/(rdry*env_temp) * exp( - grav/(rdry*env_temp) * cz(k) )
1761  momz(k,i,j) = 0.0_rp
1762  momx(k,i,j) = env_u * dens(k,i,j)
1763  momy(k,i,j) = env_v * dens(k,i,j)
1764 
1765  ! make pressure bubble
1766  pres(k,i,j) = dens(k,i,j) * env_temp * rdry + bbl_pres * bubble(k,i,j)
1767 
1768  rhot(k,i,j) = dens(k,i,j) * env_temp * ( p00/pres(k,i,j) )**rovcp
1769  enddo
1770  enddo
1771  enddo
1772 
1773  return
1774  end subroutine mkinit_lambwave
1775 
1776  !-----------------------------------------------------------------------------
1779  subroutine mkinit_gravitywave
1780  implicit none
1781 
1782  ! Surface state
1783  real(RP) :: SFC_THETA ! surface potential temperature [K]
1784  real(RP) :: SFC_PRES ! surface pressure [Pa]
1785  ! Environment state
1786  real(RP) :: ENV_U = 20.0_rp ! velocity u of environment [m/s]
1787  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1788  real(RP) :: ENV_BVF = 0.01_rp ! Brunt Vaisala frequencies of environment [1/s]
1789  ! Bubble
1790  real(RP) :: BBL_THETA = 0.01_rp ! extremum of potential temperature in bubble [K]
1791 
1792  namelist / param_mkinit_gravitywave / &
1793  sfc_theta, &
1794  sfc_pres, &
1795  env_u, &
1796  env_v, &
1797  env_bvf, &
1798  bbl_theta
1799 
1800  integer :: ierr
1801  integer :: k, i, j
1802  !---------------------------------------------------------------------------
1803 
1804  log_newline
1805  log_info("MKINIT_gravitywave",*) 'Setup initial state'
1806 
1807  sfc_theta = thetastd
1808  sfc_pres = pstd
1809 
1810  !--- read namelist
1811  rewind(io_fid_conf)
1812  read(io_fid_conf,nml=param_mkinit_gravitywave,iostat=ierr)
1813 
1814  if( ierr < 0 ) then !--- missing
1815  log_info("MKINIT_gravitywave",*) 'Not found namelist. Default used.'
1816  elseif( ierr > 0 ) then !--- fatal error
1817  log_error("MKINIT_gravitywave",*) 'Not appropriate names in namelist PARAM_MKINIT_GRAVITYWAVE. Check!'
1818  call prc_abort
1819  endif
1820  log_nml(param_mkinit_gravitywave)
1821 
1822  ! calc in dry condition
1823  pres_sfc(1,1) = sfc_pres
1824  pott_sfc(1,1) = sfc_theta
1825 
1826  do k = ks, ke
1827  pott(k,1,1) = sfc_theta * exp( env_bvf*env_bvf / grav * cz(k) )
1828  enddo
1829 
1830  ! make density & pressure profile in dry condition
1831  call hydrostatic_buildrho( ka, ks, ke, &
1832  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
1833  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
1834  cz(:), fz(:), & ! [IN]
1835  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
1836 
1837  do j = jsb, jeb
1838  do i = isb, ieb
1839  do k = ks, ke
1840  dens(k,i,j) = dens(k,1,1)
1841  momz(k,i,j) = 0.0_rp
1842  momx(k,i,j) = env_u * dens(k,1,1)
1843  momy(k,i,j) = env_v * dens(k,1,1)
1844 
1845  ! make warm bubble
1846  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
1847 
1848  enddo
1849  enddo
1850  enddo
1851 
1852  return
1853  end subroutine mkinit_gravitywave
1854 
1855  !-----------------------------------------------------------------------------
1857  subroutine mkinit_khwave
1858  implicit none
1859 
1860  ! Surface state
1861  real(RP) :: SFC_THETA ! surface potential temperature [K]
1862  real(RP) :: SFC_PRES ! surface pressure [Pa]
1863  ! Environment state
1864  real(RP) :: ENV_L1_ZTOP = 1900.0_rp ! top height of the layer1 (low THETA) [m]
1865  real(RP) :: ENV_L3_ZBOTTOM = 2100.0_rp ! bottom height of the layer3 (high THETA) [m]
1866  real(RP) :: ENV_L1_THETA = 300.0_rp ! THETA in the layer1 (low THETA) [K]
1867  real(RP) :: ENV_L3_THETA = 301.0_rp ! THETA in the layer3 (high THETA) [K]
1868  real(RP) :: ENV_L1_U = 0.0_rp ! velocity u in the layer1 (low THETA) [K]
1869  real(RP) :: ENV_L3_U = 20.0_rp ! velocity u in the layer3 (high THETA) [K]
1870  ! Disturbance
1871  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
1872 
1873  namelist / param_mkinit_khwave / &
1874  sfc_theta, &
1875  sfc_pres, &
1876  env_l1_ztop, &
1877  env_l3_zbottom, &
1878  env_l1_theta, &
1879  env_l3_theta, &
1880  env_l1_u, &
1881  env_l3_u, &
1882  random_u
1883 
1884  real(RP) :: fact
1885 
1886  integer :: ierr
1887  integer :: k, i, j
1888  !---------------------------------------------------------------------------
1889 
1890  log_newline
1891  log_info("MKINIT_khwave",*) 'Setup initial state'
1892 
1893  sfc_theta = thetastd
1894  sfc_pres = pstd
1895 
1896  !--- read namelist
1897  rewind(io_fid_conf)
1898  read(io_fid_conf,nml=param_mkinit_khwave,iostat=ierr)
1899 
1900  if( ierr < 0 ) then !--- missing
1901  log_info("MKINIT_khwave",*) 'Not found namelist. Default used.'
1902  elseif( ierr > 0 ) then !--- fatal error
1903  log_error("MKINIT_khwave",*) 'Not appropriate names in namelist PARAM_MKINIT_KHWAVE. Check!'
1904  call prc_abort
1905  endif
1906  log_nml(param_mkinit_khwave)
1907 
1908  ! calc in dry condition
1909  pres_sfc(1,1) = sfc_pres
1910  pott_sfc(1,1) = sfc_theta
1911 
1912  do k = ks, ke
1913  fact = ( cz(k)-env_l1_ztop ) / ( env_l3_zbottom-env_l1_ztop )
1914  fact = max( min( fact, 1.0_rp ), 0.0_rp )
1915 
1916  pott(k,1,1) = env_l1_theta * ( 1.0_rp - fact ) &
1917  + env_l3_theta * ( fact )
1918  enddo
1919 
1920  ! make density & pressure profile in dry condition
1921  call hydrostatic_buildrho( ka, ks, ke, &
1922  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
1923  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
1924  cz(:), fz(:), & ! [IN]
1925  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
1926 
1927  do j = jsb, jeb
1928  do i = isb, ieb
1929  do k = ks, ke
1930  dens(k,i,j) = dens(k,1,1)
1931  momz(k,i,j) = 0.0_rp
1932  momy(k,i,j) = 0.0_rp
1933  rhot(k,i,j) = dens(k,1,1) * pott(k,1,1)
1934  enddo
1935  enddo
1936  enddo
1937 
1938  call random_uniform(rndm) ! make random
1939  do j = jsb, jeb
1940  do i = isb, ieb
1941  do k = ks, ke
1942  fact = ( cz(k)-env_l1_ztop ) / ( env_l3_zbottom-env_l1_ztop )
1943  fact = max( min( fact, 1.0_rp ), 0.0_rp )
1944 
1945  momx(k,i,j) = ( env_l1_u * ( 1.0_rp - fact ) &
1946  + env_l3_u * ( fact ) &
1947  + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u &
1948  ) * dens(k,i,j)
1949  enddo
1950  enddo
1951  enddo
1952 
1953  return
1954  end subroutine mkinit_khwave
1955 
1956  !-----------------------------------------------------------------------------
1958  subroutine mkinit_turbulence
1959  use scale_atmos_hydrometeor, only: &
1961  implicit none
1962 
1963  ! Surface state
1964  real(RP) :: SFC_THETA ! surface potential temperature [K]
1965  real(RP) :: SFC_PRES ! surface pressure [Pa]
1966  real(RP) :: SFC_RH = 0.0_rp ! surface relative humidity [%]
1967  ! Environment state
1968  real(RP) :: ENV_THETA ! potential temperature of environment [K]
1969  real(RP) :: ENV_TLAPS = 4.e-3_rp ! Lapse rate of THETA [K/m]
1970  real(RP) :: ENV_U = 5.0_rp ! velocity u of environment [m/s]
1971  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
1972  real(RP) :: ENV_RH = 0.0_rp ! relative humidity of environment [%]
1973  ! Disturbance
1974  real(RP) :: RANDOM_THETA = 1.0_rp ! amplitude of random disturbance theta
1975  real(RP) :: RANDOM_U = 0.0_rp ! amplitude of random disturbance u
1976  real(RP) :: RANDOM_V = 0.0_rp ! amplitude of random disturbance v
1977  real(RP) :: RANDOM_RH = 0.0_rp ! amplitude of random disturbance RH
1978 
1979  namelist / param_mkinit_turbulence / &
1980  sfc_theta, &
1981  sfc_pres, &
1982  sfc_rh, &
1983  env_theta, &
1984  env_tlaps, &
1985  env_u, &
1986  env_v, &
1987  env_rh, &
1988  random_theta, &
1989  random_u, &
1990  random_v, &
1991  random_rh
1992 
1993  integer :: ierr
1994  integer :: k, i, j
1995  !---------------------------------------------------------------------------
1996 
1997  log_newline
1998  log_info("MKINIT_turbulence",*) 'Setup initial state'
1999 
2000  sfc_theta = thetastd
2001  sfc_pres = pstd
2002  env_theta = thetastd
2003 
2004  !--- read namelist
2005  rewind(io_fid_conf)
2006  read(io_fid_conf,nml=param_mkinit_turbulence,iostat=ierr)
2007 
2008  if( ierr < 0 ) then !--- missing
2009  log_info("MKINIT_turbulence",*) 'Not found namelist. Default used.'
2010  elseif( ierr > 0 ) then !--- fatal error
2011  log_error("MKINIT_turbulence",*) 'Not appropriate names in namelist PARAM_MKINIT_TURBULENCE. Check!'
2012  call prc_abort
2013  endif
2014  log_nml(param_mkinit_turbulence)
2015 
2016  ! calc in dry condition
2017  pres_sfc(1,1) = sfc_pres
2018  pott_sfc(1,1) = sfc_theta
2019 
2020  do k = ks, ke
2021  pott(k,1,1) = env_theta + env_tlaps * cz(k)
2022  enddo
2023 
2024  ! make density & pressure profile in dry condition
2025  call hydrostatic_buildrho( ka, ks, ke, &
2026  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
2027  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
2028  cz(:), fz(:), & ! [IN]
2029  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
2030 
2031  if ( .not. atmos_hydrometeor_dry ) then
2032  ! calc QV from RH
2033  call saturation_psat_all( temp_sfc(1,1), & ! [IN]
2034  psat_sfc(1,1) ) ! [OUT]
2035  qdry(:,1,1) = 1.0_rp - qv(:,1,1) - qc(:,1,1)
2036  call saturation_pres2qsat_all( ka, ks, ke, &
2037  temp(:,1,1), pres(:,1,1), qdry(:,1,1), & ! [IN]
2038  qsat(:,1,1) ) ! [OUT]
2039 
2040  call random_uniform(rndm) ! make random
2041  do j = jsb, jeb
2042  do i = isb, ieb
2043  qsat_sfc(1,1) = epsvap * psat_sfc(i,j) / ( pres_sfc(i,j) - ( 1.0_rp-epsvap ) * psat_sfc(i,j) )
2044  qv_sfc(i,j) = ( sfc_rh + rndm(ks-1,i,j) * random_rh ) * 1.e-2_rp * qsat_sfc(1,1)
2045 
2046  do k = ks, ke
2047  qv(k,i,j) = ( env_rh + rndm(k,i,j) * random_rh ) * 1.e-2_rp * qsat(k,1,1)
2048  enddo
2049  enddo
2050  enddo
2051  end if
2052 
2053  call random_uniform(rndm) ! make random
2054  do j = jsb, jeb
2055  do i = isb, ieb
2056  pres_sfc(i,j) = sfc_pres
2057  pott_sfc(i,j) = sfc_theta + rndm(ks-1,i,j) * random_theta
2058 
2059  do k = ks, ke
2060  pott(k,i,j) = env_theta + env_tlaps * cz(k) + rndm(k,i,j) * random_theta
2061  enddo
2062  enddo
2063  enddo
2064 
2065  ! make density & pressure profile in moist condition
2066  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2067  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
2068  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
2069  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
2070  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
2071 
2072  call comm_vars8( dens(:,:,:), 1 )
2073  call comm_wait ( dens(:,:,:), 1 )
2074 
2075  call random_uniform(rndm) ! make random
2076  do j = jsb, jeb
2077  do i = isb, ieb
2078  do k = ks, ke
2079  momx(k,i,j) = ( env_u + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_u ) &
2080  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
2081  enddo
2082  enddo
2083  enddo
2084 
2085  call random_uniform(rndm) ! make random
2086  do j = jsb, jeb
2087  do i = isb, ieb
2088  do k = ks, ke
2089  momy(k,i,j) = ( env_v + ( rndm(k,i,j) - 0.5_rp ) * 2.0_rp * random_v ) &
2090  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
2091  enddo
2092  enddo
2093  enddo
2094 
2095  do j = jsb, jeb
2096  do i = isb, ieb
2097  do k = ks, ke
2098  momz(k,i,j) = 0.0_rp
2099  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
2100  enddo
2101  enddo
2102  enddo
2103 
2104  return
2105  end subroutine mkinit_turbulence
2106 
2107  !-----------------------------------------------------------------------------
2109  subroutine mkinit_cavityflow
2110  implicit none
2111 
2112  ! Nondimenstional numbers for a cavity flow problem
2113  real(RP) :: REYNOLDS_NUM = 1.d03
2114  real(RP) :: MACH_NUM = 3.d-2
2115  real(RP) :: Ulid = 1.d01
2116  real(RP) :: PRES0 = 1.d05
2117 
2118  namelist / param_mkinit_cavityflow / &
2119  ulid , &
2120  pres0 , &
2121  reynolds_num, &
2122  mach_num
2123 
2124  real(RP) :: DENS0
2125  real(RP) :: TEMP
2126  real(RP) :: Gam
2127  real(RP) :: Cs2
2128 
2129  integer :: k, i, j
2130  integer :: ierr
2131  !---------------------------------------------------------------------------
2132 
2133  log_newline
2134  log_info("MKINIT_cavityflow",*) 'Setup initial state'
2135 
2136  !--- read namelist
2137  rewind(io_fid_conf)
2138  read(io_fid_conf,nml=param_mkinit_cavityflow,iostat=ierr)
2139 
2140  if( ierr < 0 ) then !--- missing
2141  log_info("MKINIT_cavityflow",*) 'Not found namelist. Default used.'
2142  elseif( ierr > 0 ) then !--- fatal error
2143  log_error("MKINIT_cavityflow",*) 'Not appropriate names in namelist PARAM_MKINIT_CAVITYFLOW. Check!'
2144  call prc_abort
2145  endif
2146  log_nml(param_mkinit_cavityflow)
2147 
2148  gam = cpdry / ( cpdry - rdry )
2149  cs2 = ( ulid / mach_num )**2
2150  temp = cs2 / ( gam * rdry )
2151  dens0 = pres0 / ( rdry * temp )
2152 
2153  log_info("MKINIT_cavityflow",*) "DENS = ", dens0
2154  log_info("MKINIT_cavityflow",*) "PRES = ", pres0
2155  log_info("MKINIT_cavityflow",*) "TEMP = ", rhot(10,10,4)/dens0, temp
2156  log_info("MKINIT_cavityflow",*) "Ulid = ", ulid
2157  log_info("MKINIT_cavityflow",*) "Cs = ", sqrt(cs2)
2158 
2159  do j = 1, ja
2160  do i = 1, ia
2161  do k = ks, ke
2162  dens(k,i,j) = dens0
2163  momz(k,i,j) = 0.0_rp
2164  momx(k,i,j) = 0.0_rp
2165  momy(k,i,j) = 0.0_rp
2166  pres(k,i,j) = pres0
2167  rhot(k,i,j) = p00/rdry * (p00/pres0)**((rdry - cpdry)/cpdry)
2168  enddo
2169  enddo
2170  enddo
2171 
2172  momx(ke+1:ka,:,:) = dens0 * ulid
2173 
2174  return
2175  end subroutine mkinit_cavityflow
2176 
2177  !-----------------------------------------------------------------------------
2179  subroutine mkinit_mountainwave
2180  implicit none
2181 
2182  ! Surface state
2183  real(RP) :: SFC_THETA ! surface potential temperature [K]
2184  real(RP) :: SFC_PRES ! surface pressure [Pa]
2185  ! Environment state
2186  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
2187  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2188 
2189  real(RP) :: SCORER = 2.e-3_rp ! Scorer parameter (~=N/U) [1/m]
2190  real(RP) :: BBL_NC = 0.0_rp ! extremum of NC in bubble [kg/kg]
2191 
2192  namelist / param_mkinit_mountainwave / &
2193  sfc_theta, &
2194  sfc_pres, &
2195  env_u, &
2196  env_v, &
2197  scorer, &
2198  bbl_nc
2199 
2200  real(RP) :: Ustar2, N2
2201 
2202  integer :: ierr
2203  integer :: k, i, j
2204  !---------------------------------------------------------------------------
2205 
2206  log_newline
2207  log_info("MKINIT_mountainwave",*) 'Setup initial state'
2208 
2209  sfc_theta = thetastd
2210  sfc_pres = pstd
2211 
2212  !--- read namelist
2213  rewind(io_fid_conf)
2214  read(io_fid_conf,nml=param_mkinit_mountainwave,iostat=ierr)
2215 
2216  if( ierr < 0 ) then !--- missing
2217  log_info("MKINIT_mountainwave",*) 'Not found namelist. Default used.'
2218  elseif( ierr > 0 ) then !--- fatal error
2219  log_error("MKINIT_mountainwave",*) 'Not appropriate names in namelist PARAM_MKINIT_MOUNTAINWAVE. Check!'
2220  call prc_abort
2221  endif
2222  log_nml(param_mkinit_mountainwave)
2223 
2224  ! calc in dry condition
2225  do j = jsb, jeb
2226  do i = isb, ieb
2227  pres_sfc(i,j) = sfc_pres
2228  pott_sfc(i,j) = sfc_theta
2229  enddo
2230  enddo
2231 
2232  do j = jsb, jeb
2233  do i = isb, ieb
2234  do k = ks, ke
2235  ustar2 = env_u * env_u + env_v * env_v
2236  n2 = ustar2 * (scorer*scorer)
2237 
2238  pott(k,i,j) = sfc_theta * exp( n2 / grav * real_cz(k,i,j) )
2239  enddo
2240  enddo
2241  enddo
2242 
2243  ! make density & pressure profile in dry condition
2244  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2245  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
2246  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
2247  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
2248  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
2249 
2250  do j = jsb, jeb
2251  do i = isb, ieb
2252  do k = ks, ke
2253  dens(k,i,j) = dens(k,i,j)
2254  momz(k,i,j) = 0.0_rp
2255  momx(k,i,j) = env_u * dens(k,i,j)
2256  momy(k,i,j) = env_v * dens(k,i,j)
2257  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
2258  enddo
2259  enddo
2260  enddo
2261 
2262  ! optional : add tracer bubble
2263  if ( bbl_nc > 0.0_rp ) then
2264  do j = jsb, jeb
2265  do i = isb, ieb
2266  do k = ks, ke
2267  nc(k,i,j) = bbl_nc * bubble(k,i,j)
2268  enddo
2269  enddo
2270  enddo
2271  endif
2272 
2273  return
2274  end subroutine mkinit_mountainwave
2275 
2276  !-----------------------------------------------------------------------------
2281  subroutine mkinit_barocwave
2282  use scale_const, only: &
2283  ohm => const_ohm, &
2284  rplanet => const_radius, &
2285  grav => const_grav
2286  use scale_prc
2287  use scale_atmos_grid_cartesc, only: &
2289  fyg => atmos_grid_cartesc_fyg
2290  use scale_atmos_hydrometeor, only: &
2292 
2293  implicit none
2294 
2295  ! Parameters for global domain size
2296  real(RP) :: Ly ! The domain size in y-direction [m]
2297 
2298  ! Parameters for inital stratification
2299  real(RP) :: REF_TEMP = 288.e0_rp ! The reference temperature [K]
2300  real(RP) :: REF_PRES = 1.e5_rp ! The reference pressure [Pa]
2301  real(RP) :: LAPSE_RATE = 5.e-3_rp ! The lapse rate [K/m]
2302 
2303  ! Parameters associated with coriolis parameter on a beta-plane
2304  real(RP) :: Phi0Deg = 45.e0_rp ! The central latitude [degree_north]
2305 
2306  ! Parameters for background zonal jet
2307  real(RP) :: U0 = 35.e0_rp ! The parameter associated with zonal jet maximum amplitude [m/s]
2308  real(RP) :: b = 2.e0_rp ! The vertical half-width [1]
2309 
2310  ! Parameters for inital perturbation of zonal wind with a Gaussian profile
2311  !
2312  real(RP) :: Up = 1.e0_rp ! The maximum amplitude of zonal wind perturbation [m/s]
2313  real(RP) :: Lp = 600.e3_rp ! The width of Gaussian profile
2314  real(RP) :: Xc = 2000.e3_rp ! The center point (x) of inital perturbation
2315  real(RP) :: Yc = 2500.e3_rp ! The center point (y) of inital perturbation
2316 
2317  namelist / param_mkinit_barocwave / &
2318  ref_temp, ref_pres, lapse_rate, &
2319  phi0deg, &
2320  u0, b, &
2321  up, lp, xc, yc
2322 
2323  real(RP) :: f0, beta0
2324 
2325  real(RP) :: geopot(ka,ia,ja)
2326  real(RP) :: eta(ka,ia,ja)
2327  real(RP) :: temp(ka,ia,ja)
2328 
2329  real(RP) :: y
2330  real(RP) :: ln_eta
2331  real(RP) :: del_eta
2332  real(RP) :: yphase
2333  real(RP) :: temp_vfunc
2334  real(RP) :: geopot_hvari
2335 
2336  integer :: ierr
2337  integer :: k, i, j
2338 
2339  integer :: itr
2340 
2341  integer, parameter :: ITRMAX = 1000
2342  real(RP), parameter :: CONV_EPS = 1e-15_rp
2343  !---------------------------------------------------------------------------
2344 
2345  log_newline
2346  log_info("MKINIT_barocwave",*) 'Setup initial state'
2347 
2348  !--- read namelist
2349  rewind(io_fid_conf)
2350  read(io_fid_conf,nml=param_mkinit_barocwave,iostat=ierr)
2351 
2352  if( ierr < 0 ) then !--- missing
2353  log_info("MKINIT_barocwave",*) 'Not found namelist. Default used.'
2354  elseif( ierr > 0 ) then !--- fatal error
2355  log_error("MKINIT_barocwave",*) 'Not appropriate names in namelist PARAM_MKINIT_BAROCWAVE. Check!'
2356  call prc_abort
2357  endif
2358  log_nml(param_mkinit_barocwave)
2359 
2360  ly = fyg(jag-jhalo) - fyg(jhalo)
2361 
2362  ! Set coriolis parameters
2363  f0 = 2.0_rp*ohm*sin(phi0deg*pi/180.0_rp)
2364  beta0 = (2.0_rp*ohm/rplanet)*cos(phi0deg*pi/180.0_rp)
2365 
2366  ! Calculate eta(=p/p_s) level corresponding to z level of each (y,z) grid point
2367  ! using Newton's iteration method
2368 
2369  eta(:,:,:) = 1.0e-8_rp ! Set first guess of eta
2370 
2371  do j = jsb, jeb
2372  do i = isb, ieb ! Note that initial fields are zonaly symmetric
2373 
2374  y = cy(j)
2375  yphase = 2.0_rp*pi*y/ly
2376 
2377  ! Calc horizontal variation of geopotential height
2378  geopot_hvari = 0.5_rp*u0*( &
2379  (f0 - beta0*y0)*(y - 0.5_rp*ly*(1.0_rp + sin(yphase)/pi)) &
2380  + 0.5_rp*beta0*( y**2 - ly*y/pi*sin(yphase) - 0.5_rp*(ly/pi)**2*(cos(yphase) + 1.0_rp) &
2381  - ly**2/3.0_rp ) &
2382  )
2383 
2384  ! Set surface pressure and temperature
2385  pres_sfc(i,j) = ref_pres
2386  pott_sfc(i,j) = ref_temp - geopot_hvari/rdry
2387 
2388  do k = ks, ke
2389  del_eta = 1.0_rp
2390 
2391  !-- The loop for iteration
2392  itr = 0
2393  do while( abs(del_eta) > conv_eps )
2394  ln_eta = log(eta(k,i,j))
2395 
2396  temp_vfunc = eta(k,i,j)**(rdry*lapse_rate/grav)
2397  temp(k,i,j) = &
2398  ref_temp*temp_vfunc &
2399  + geopot_hvari/rdry*(2.0_rp*(ln_eta/b)**2 - 1.0_rp)*exp(-(ln_eta/b)**2)
2400  geopot(k,i,j) = &
2401  ref_temp*grav/lapse_rate*(1.0_rp - temp_vfunc) &
2402  + geopot_hvari*ln_eta*exp(-(ln_eta/b)**2)
2403 
2404  del_eta = - ( - grav*cz(k) + geopot(k,i,j) ) & ! <- F
2405  & *( - eta(k,i,j)/(rdry*temp(k,i,j)) ) ! <- (dF/deta)^-1
2406 
2407  eta(k,i,j) = eta(k,i,j) + del_eta
2408  itr = itr + 1
2409 
2410  if ( itr > itrmax ) then
2411  log_error("MKINIT_barocwave",*) "Fail the convergence of iteration. Check!"
2412  log_error_cont(*) "* (X,Y,Z)=", cx(i), cy(j), cz(k)
2413  log_error_cont(*) "itr=", itr, "del_eta=", del_eta, "eta=", eta(k,i,j), "temp=", temp(k,i,j)
2414  call prc_abort
2415  end if
2416  enddo !- End of loop for iteration ----------------------------
2417 
2418  pres(k,i,j) = eta(k,i,j)*ref_pres
2419  dens(k,i,j) = pres(k,i,j)/(rdry*temp(k,i,j))
2420  pott(k,i,j) = temp(k,i,j)*eta(k,i,j)**(-rdry/cpdry)
2421 
2422  enddo
2423 
2424  ! Make density & pressure profile in dry condition using the profile of
2425  ! potential temperature calculated above.
2426  call hydrostatic_buildrho( ka, ks, ke, &
2427  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
2428  pres_sfc(i,j), pott_sfc(i,j), qv_sfc(i,j), qc_sfc(i,j), & ! [IN]
2429  real_cz(:,i,j), real_fz(:,i,j), & ! [IN]
2430  dens(:,i,j), temp(:,i,j), pres(:,i,j), temp_sfc(i,j) ) ! [OUT]
2431  enddo
2432  enddo
2433 
2434  !-----------------------------------------------------------------------------------
2435 
2436  do j = jsb, jeb
2437  do k = ks, ke
2438 
2439  eta(k,is,j) = pres(k,is,j)/ref_pres
2440  ln_eta = log(eta(k,is,j))
2441  yphase = 2.0_rp*pi*cy(j)/ly
2442 !!$ PRES(k,IS:IE,j) = eta(k,IS,j)*REF_PRES
2443 !!$ DENS(k,IS:IE,j) = PRES(k,IS,j)/(Rdry*temp(k,IS,j))
2444  dens(k,is:ie,j) = dens(k,is,j)
2445  pres(k,is:ie,j) = pres(k,is,j)
2446  momx(k,is-1:ie,j) = dens(k,is,j)*(-u0*sin(0.5_rp*yphase)**2*ln_eta*exp(-(ln_eta/b)**2))
2447  rhot(k,is:ie,j) = dens(k,is,j)*pott(k,is,j) !temp(k,IS,j)*eta(k,IS,j)**(-Rdry/CPdry)
2448  enddo
2449  enddo
2450  momy(:,:,:) = 0.0_rp
2451  momz(:,:,:) = 0.0_rp
2452 
2453  !---------------------------------------------------------------------------------------
2454 
2455  ! Add the inital perturbation for zonal velocity
2456  do j = jsb, jeb
2457  do i = max(isb-1,1), ieb
2458  momx(ks:ke,i,j) = momx(ks:ke,i,j) &
2459  + dens(ks:ke,i,j)* up*exp( - ((fx(i) - xc)**2 + (cy(j) - yc)**2)/lp**2 )
2460  enddo
2461  enddo
2462 
2463  return
2464  end subroutine mkinit_barocwave
2465 
2466  !-----------------------------------------------------------------------------
2468  subroutine mkinit_warmbubble
2469  use scale_atmos_hydrometeor, only: &
2471  implicit none
2472 
2473  ! Surface state
2474  real(RP) :: SFC_THETA ! surface potential temperature [K]
2475  real(RP) :: SFC_PRES ! surface pressure [Pa]
2476  real(RP) :: SFC_RH = 80.0_rp ! surface relative humidity [%]
2477  ! Environment state
2478  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
2479  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
2480  real(RP) :: ENV_RH = 80.0_rp ! Relative Humidity of environment [%]
2481  real(RP) :: ENV_L1_ZTOP = 1.e3_rp ! top height of the layer1 (constant THETA) [m]
2482  real(RP) :: ENV_L2_ZTOP = 14.e3_rp ! top height of the layer2 (small THETA gradient) [m]
2483  real(RP) :: ENV_L2_TLAPS = 4.e-3_rp ! Lapse rate of THETA in the layer2 (small THETA gradient) [K/m]
2484  real(RP) :: ENV_L3_TLAPS = 3.e-2_rp ! Lapse rate of THETA in the layer3 (large THETA gradient) [K/m]
2485  ! Bubble
2486  real(RP) :: BBL_THETA = 1.0_rp ! extremum of temperature in bubble [K]
2487 
2488  namelist / param_mkinit_warmbubble / &
2489  sfc_theta, &
2490  sfc_pres, &
2491  env_u, &
2492  env_v, &
2493  env_rh, &
2494  env_l1_ztop, &
2495  env_l2_ztop, &
2496  env_l2_tlaps, &
2497  env_l3_tlaps, &
2498  bbl_theta
2499 
2500  integer :: ierr
2501  integer :: k, i, j
2502  !---------------------------------------------------------------------------
2503 
2504  log_newline
2505  log_info("MKINIT_warmbubble",*) 'Setup initial state'
2506 
2507  if ( atmos_hydrometeor_dry ) then
2508  log_error("MKINIT_warmbubble",*) 'QV is not registered'
2509  call prc_abort
2510  end if
2511 
2512  sfc_theta = thetastd
2513  sfc_pres = pstd
2514 
2515  !--- read namelist
2516  rewind(io_fid_conf)
2517  read(io_fid_conf,nml=param_mkinit_warmbubble,iostat=ierr)
2518 
2519  if( ierr < 0 ) then !--- missing
2520  log_info("MKINIT_warmbubble",*) 'Not found namelist. Default used.'
2521  elseif( ierr > 0 ) then !--- fatal error
2522  log_error("MKINIT_warmbubble",*) 'Not appropriate names in namelist PARAM_MKINIT_WARMBUBBLE. Check!'
2523  call prc_abort
2524  endif
2525  log_nml(param_mkinit_warmbubble)
2526 
2527  ! calc in dry condition
2528  pres_sfc(1,1) = sfc_pres
2529  pott_sfc(1,1) = sfc_theta
2530 
2531  do k = ks, ke
2532  if ( cz(k) <= env_l1_ztop ) then ! Layer 1
2533  pott(k,1,1) = sfc_theta
2534  elseif( cz(k) < env_l2_ztop ) then ! Layer 2
2535  pott(k,1,1) = pott(k-1,1,1) + env_l2_tlaps * ( cz(k)-cz(k-1) )
2536  else ! Layer 3
2537  pott(k,1,1) = pott(k-1,1,1) + env_l3_tlaps * ( cz(k)-cz(k-1) )
2538  endif
2539  enddo
2540 
2541  ! make density & pressure profile in dry condition
2542  call hydrostatic_buildrho( ka, ks, ke, &
2543  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
2544  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
2545  cz(:), fz(:), & ! [IN]
2546  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
2547 
2548  ! calc QV from RH
2549  call saturation_psat_all( temp_sfc(1,1), & ! [IN]
2550  psat_sfc(1,1) ) ! [OUT]
2551  qsat_sfc(1,1) = epsvap * psat_sfc(1,1) / ( pres_sfc(1,1) - ( 1.0_rp-epsvap ) * psat_sfc(1,1) )
2552  qv_sfc(1,1) = sfc_rh * 1.e-2_rp * qsat_sfc(1,1)
2553  qdry(:,1,1) = 1.0_rp - qv(:,1,1) - qc(:,1,1)
2554  call saturation_pres2qsat_all( ka, ks, ke, &
2555  temp(:,1,1), pres(:,1,1), qdry(:,1,1), & ! [IN]
2556  qsat(:,1,1) ) ! [OUT]
2557  do k = ks, ke
2558  if ( cz(k) <= env_l1_ztop ) then ! Layer 1
2559  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
2560  elseif( cz(k) <= env_l2_ztop ) then ! Layer 2
2561  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
2562  else ! Layer 3
2563  qv(k,1,1) = 0.0_rp
2564  endif
2565  enddo
2566 
2567  ! make density & pressure profile in moist condition
2568  call hydrostatic_buildrho( ka, ks, ke, &
2569  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
2570  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
2571  cz(:), fz(:), & ! [IN]
2572  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
2573 
2574  do j = jsb, jeb
2575  do i = isb, ieb
2576  do k = ks, ke
2577  dens(k,i,j) = dens(k,1,1)
2578  momz(k,i,j) = 0.0_rp
2579  momx(k,i,j) = env_u * dens(k,i,j)
2580  momy(k,i,j) = env_v * dens(k,i,j)
2581 
2582  ! make warm bubble
2583  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
2584 
2585  qv(k,i,j) = qv(k,1,1)
2586  enddo
2587  enddo
2588  enddo
2589 
2590  call flux_setup
2591 
2592  return
2593  end subroutine mkinit_warmbubble
2594 
2595  !-----------------------------------------------------------------------------
2597  subroutine mkinit_supercell
2598  use scale_atmos_hydrometeor, only: &
2600  implicit none
2601 
2602  real(RP) :: RHO(ka)
2603  real(RP) :: VELX(ka)
2604  real(RP) :: VELY(ka)
2605  real(RP) :: POTT(ka)
2606  real(RP) :: QV1D(ka)
2607 
2608  ! Bubble
2609  real(RP) :: BBL_THETA = 3.d0 ! extremum of temperature in bubble [K]
2610 
2611  namelist / param_mkinit_supercell / &
2612  bbl_theta
2613 
2614  integer :: ierr
2615  integer :: k, i, j
2616  !---------------------------------------------------------------------------
2617 
2618  log_newline
2619  log_info("MKINIT_supercell",*) 'Setup initial state'
2620 
2621  if ( atmos_hydrometeor_dry ) then
2622  log_error("MKINIT_supercell",*) 'QV is not registered'
2623  call prc_abort
2624  end if
2625 
2626  !--- read namelist
2627  rewind(io_fid_conf)
2628  read(io_fid_conf,nml=param_mkinit_supercell,iostat=ierr)
2629 
2630  if( ierr < 0 ) then !--- missing
2631  log_info("MKINIT_supercell",*) 'Not found namelist. Default used.'
2632  elseif( ierr > 0 ) then !--- fatal error
2633  log_error("MKINIT_supercell",*) 'Not appropriate names in namelist PARAM_MKINIT_SUPERCELL. Check!'
2634  call prc_abort
2635  endif
2636  log_nml(param_mkinit_supercell)
2637 
2638  call read_sounding( rho, velx, vely, pott, qv1d ) ! (out)
2639 
2640  do j = jsb, jeb
2641  do i = isb, ieb
2642  do k = ks, ke
2643  dens(k,i,j) = rho(k)
2644  momz(k,i,j) = 0.0_rp
2645  momx(k,i,j) = rho(k) * velx(k)
2646  momy(k,i,j) = rho(k) * vely(k)
2647 
2648  ! make warm bubble
2649  rhot(k,i,j) = rho(k) * ( pott(k) + bbl_theta * bubble(k,i,j) )
2650 
2651  qv(k,i,j) = qv1d(k)
2652  enddo
2653  enddo
2654  enddo
2655 
2656  call flux_setup
2657 
2658  return
2659  end subroutine mkinit_supercell
2660 
2661  !-----------------------------------------------------------------------------
2663  subroutine mkinit_squallline
2664  use scale_atmos_hydrometeor, only: &
2666  implicit none
2667 
2668  real(RP) :: RHO(ka)
2669  real(RP) :: VELX(ka)
2670  real(RP) :: VELY(ka)
2671  real(RP) :: POTT(ka)
2672  real(RP) :: QV1D(ka)
2673 
2674  real(RP) :: RANDOM_THETA = 0.01_rp
2675  real(RP) :: OFFSET_velx = 12.0_rp
2676  real(RP) :: OFFSET_vely = -2.0_rp
2677 
2678  namelist / param_mkinit_squallline / &
2679  random_theta, &
2680  offset_velx, &
2681  offset_vely
2682 
2683  integer :: ierr
2684  integer :: k, i, j
2685  !---------------------------------------------------------------------------
2686 
2687  log_newline
2688  log_info("MKINIT_squallline",*) 'Setup initial state'
2689 
2690  if ( atmos_hydrometeor_dry ) then
2691  log_error("MKINIT_squallline",*) 'QV is not registered'
2692  call prc_abort
2693  end if
2694 
2695  !--- read namelist
2696  rewind(io_fid_conf)
2697  read(io_fid_conf,nml=param_mkinit_squallline,iostat=ierr)
2698 
2699  if( ierr < 0 ) then !--- missing
2700  log_info("MKINIT_squallline",*) 'Not found namelist. Default used.'
2701  elseif( ierr > 0 ) then !--- fatal error
2702  log_error("MKINIT_squallline",*) 'Not appropriate names in namelist PARAM_MKINIT_SQUALLLINE. Check!'
2703  call prc_abort
2704  endif
2705  log_nml(param_mkinit_squallline)
2706 
2707  call read_sounding( rho, velx, vely, pott, qv1d ) ! (out)
2708 
2709  call random_uniform(rndm) ! make random
2710  do j = jsb, jeb
2711  do i = isb, ieb
2712  do k = ks, ke
2713  dens(k,i,j) = rho(k)
2714  momz(k,i,j) = 0.0_rp
2715  momx(k,i,j) = ( velx(k) - offset_velx ) * rho(k)
2716  momy(k,i,j) = ( vely(k) - offset_vely ) * rho(k)
2717  rhot(k,i,j) = rho(k) * ( pott(k) + rndm(k,i,j) * random_theta )
2718  qv(k,i,j) = qv1d(k)
2719  enddo
2720  enddo
2721  enddo
2722 
2723  call flux_setup
2724 
2725  return
2726  end subroutine mkinit_squallline
2727 
2728  !-----------------------------------------------------------------------------
2730  subroutine mkinit_wk1982
2731  use scale_atmos_hydrometeor, only: &
2733  implicit none
2734 
2735  ! Surface state
2736  real(RP) :: SFC_THETA = 300.0_rp ! surface pot. temperature [K]
2737  real(RP) :: SFC_PRES ! surface pressure [Pa]
2738  ! Parameter in Weisman and Klemp (1982)
2739  real(RP) :: TR_Z = 12000.0_rp ! height of tropopause [m]
2740  real(RP) :: TR_THETA = 343.0_rp ! pot. temperature at tropopause [K]
2741  real(RP) :: TR_TEMP = 213.0_rp ! temperature at tropopause [K]
2742  real(RP) :: SHEAR_Z = 3000.0_rp ! center height of shear layer [m]
2743  real(RP) :: SHEAR_U = 15.0_rp ! velocity u over the shear layer [m/s]
2744  real(RP) :: QV0 = 14.0_rp ! maximum vapor mixing ration [g/kg]
2745  ! Bubble
2746  real(RP) :: BBL_THETA = 3.d0 ! extremum of temperature in bubble [K]
2747 
2748  namelist / param_mkinit_wk1982 / &
2749  sfc_theta, &
2750  sfc_pres, &
2751  tr_z, &
2752  tr_theta, &
2753  tr_temp, &
2754  shear_z, &
2755  shear_u, &
2756  qv0, &
2757  bbl_theta
2758 
2759  real(RP) :: rh (ka,ia,ja)
2760  real(RP) :: rh_sfc( ia,ja)
2761 
2762  integer :: ierr
2763  integer :: k, i, j
2764  !---------------------------------------------------------------------------
2765 
2766  log_newline
2767  log_info("MKINIT_wk1982",*) 'Setup initial state'
2768 
2769  if ( atmos_hydrometeor_dry ) then
2770  log_error("MKINIT_wk1982",*) 'QV is not registered'
2771  call prc_abort
2772  end if
2773 
2774  sfc_pres = pstd
2775 
2776  rewind(io_fid_conf)
2777  read(io_fid_conf,nml=param_mkinit_wk1982,iostat=ierr)
2778  if( ierr < 0 ) then !--- missing
2779  log_info("MKINIT_wk1982",*) 'Not found namelist. Default used.'
2780  elseif( ierr > 0 ) then !--- fatal error
2781  log_error("MKINIT_wk1982",*) 'Not appropriate names in namelist PARAM_MKINIT_WK1982. Check!'
2782  call prc_abort
2783  endif
2784  log_nml(param_mkinit_wk1982)
2785 
2786  ! calc in dry condition
2787  do j = jsb, jeb
2788  do i = isb, ieb
2789  pres_sfc(i,j) = sfc_pres
2790  pott_sfc(i,j) = sfc_theta
2791 
2792  do k = ks, ke
2793  if ( real_cz(k,i,j) <= tr_z ) then ! below initial cloud top
2794  pott(k,i,j) = pott_sfc(i,j) &
2795  + ( tr_theta - pott_sfc(i,j) ) * ( real_cz(k,i,j) / tr_z )**1.25_rp
2796  else
2797  pott(k,i,j) = tr_theta * exp( grav * ( real_cz(k,i,j) - tr_z ) / cpdry / tr_temp )
2798  endif
2799  enddo
2800  enddo
2801  enddo
2802 
2803  ! make density & pressure profile in dry condition
2804  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2805  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
2806  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
2807  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
2808  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
2809 
2810  ! calc QV from RH
2811  do j = jsb, jeb
2812  do i = isb, ieb
2813  rh_sfc(i,j) = 1.0_rp - 0.75_rp * ( real_fz(ks-1,i,j) / tr_z )**1.25_rp
2814 
2815  do k = ks, ke
2816  if ( real_cz(k,i,j) <= tr_z ) then ! below initial cloud top
2817  rh(k,i,j) = 1.0_rp - 0.75_rp * ( real_cz(k,i,j) / tr_z )**1.25_rp
2818  else
2819  rh(k,i,j) = 0.25_rp
2820  endif
2821  enddo
2822  enddo
2823  enddo
2824 
2825  call saturation_psat_all( ia, isb, ieb, ja, jsb, jeb, &
2826  temp_sfc(:,:), & ! [IN]
2827  psat_sfc(:,:) ) ! [OUT]
2828  qdry(:,:,:) = 1.0_rp - qv(:,:,:) - qc(:,:,:)
2829  call saturation_pres2qsat_all( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2830  temp(:,:,:), pres(:,:,:), qdry(:,:,:), & ! [IN]
2831  qsat(:,:,:) ) ! [OUT]
2832 
2833  qv0 = qv0 * 1e-3_rp ! g/kg to kg/kg
2834  qv0 = qv0 / ( 1.0_rp + qv0 ) ! mixing ratio to specicic humidity
2835 
2836  do j = jsb, jeb
2837  do i = isb, ieb
2838  qsat_sfc(i,j) = epsvap * psat_sfc(i,j) / ( pres_sfc(i,j) - ( 1.0_rp-epsvap ) * psat_sfc(i,j) )
2839  qv_sfc(i,j) = min( rh_sfc(i,j) * qsat_sfc(i,j), qv0 )
2840  do k = ks, ke
2841  qv(k,i,j) = min( rh(k,i,j) * qsat(k,i,j), qv0 )
2842  enddo
2843  enddo
2844  enddo
2845 
2846  ! make density & pressure profile in moist condition
2847  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2848  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
2849  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
2850  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
2851  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
2852 
2853  do k = ks, ke
2854  log_info("MKINIT_wk1982",*) k, real_cz(k,is,js), pres(k,is,js), pott(k,is,js), rh(k,is,js), qv(k,is,js)*1000
2855  enddo
2856 
2857  call comm_vars8( dens(:,:,:), 1 )
2858  call comm_wait ( dens(:,:,:), 1 )
2859 
2860  do j = jsb, jeb
2861  do i = isb, ieb
2862  do k = ks, ke
2863  momx(k,i,j) = shear_u * tanh( real_cz(k,i,j) / shear_z ) &
2864  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
2865  enddo
2866  enddo
2867  enddo
2868 
2869  do j = jsb, jeb
2870  do i = isb, ieb
2871  do k = ks, ke
2872  momy(k,i,j) = 0.0_rp
2873  momz(k,i,j) = 0.0_rp
2874  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
2875 
2876  ! make warm bubble
2877  rhot(k,i,j) = dens(k,i,j) * ( pott(k,i,j) + bbl_theta * bubble(k,i,j) )
2878  enddo
2879  enddo
2880  enddo
2881 
2882  call flux_setup
2883 
2884  return
2885  end subroutine mkinit_wk1982
2886 
2887  !-----------------------------------------------------------------------------
2889  subroutine mkinit_dycoms2_rf01
2890  use scale_atmos_hydrometeor, only: &
2892  implicit none
2893 
2894 #ifndef DRY
2895  real(RP) :: PERTURB_AMP = 0.0_rp
2896  integer :: RANDOM_LIMIT = 5
2897  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
2898  ! 1 -> petrurbation for pt
2899  ! 2 -> perturbation for u, v, w
2900  logical :: USE_LWSET = .false. ! use liq. water. static energy temp.?
2901 
2902  namelist / param_mkinit_rf01 / &
2903  perturb_amp, &
2904  random_limit, &
2905  random_flag, &
2906  use_lwset
2907 
2908  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
2909  real(RP) :: LHV (ka,ia,ja) ! latent heat of vaporization [J/kg]
2910 
2911  real(RP) :: qall ! QV+QC
2912  real(RP) :: fact
2913  real(RP) :: pi2
2914  real(RP) :: sint
2915  real(RP) :: GEOP_sw ! switch for geopotential energy correction
2916 
2917  integer :: ierr
2918  integer :: k, i, j
2919  !---------------------------------------------------------------------------
2920 
2921  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
2922 
2923  log_newline
2924  log_info("MKINIT_DYCOMS2_RF01",*) 'Setup initial state'
2925 
2926  rewind(io_fid_conf)
2927 
2928  if ( atmos_hydrometeor_dry ) then
2929  log_error("MKINIT_DYCOMS2_RF01",*) 'QV is not registered'
2930  call prc_abort
2931  end if
2932 
2933  read(io_fid_conf,nml=param_mkinit_rf01,iostat=ierr)
2934  if( ierr < 0 ) then !--- missing
2935  log_info("MKINIT_DYCOMS2_RF01",*) 'Not found namelist. Default used.'
2936  elseif( ierr > 0 ) then !--- fatal error
2937  log_error("MKINIT_DYCOMS2_RF01",*) 'Not appropriate names in namelist PARAM_MKINIT_RF01. Check!'
2938  call prc_abort
2939  endif
2940  log_nml(param_mkinit_rf01)
2941 
2942  if ( use_lwset ) then
2943  geop_sw = 1.0_rp
2944  else
2945  geop_sw = 0.0_rp
2946  endif
2947 
2948  ! calc in dry condition
2949  do j = jsb, jeb
2950  do i = isb, ieb
2951 
2952  pres_sfc(i,j) = 1017.8e2_rp ! [Pa]
2953  pott_sfc(i,j) = 289.0_rp ! [K]
2954 
2955  do k = ks, ke
2956  velx(k,i,j) = 7.0_rp
2957  vely(k,i,j) = -5.5_rp
2958  if ( cz(k) < 820.0_rp ) then ! below initial cloud top
2959  potl(k,i,j) = 289.0_rp - grav / cpdry * cz(k) * geop_sw
2960  elseif( cz(k) <= 860.0_rp ) then
2961  sint = sin( pi2 * ( cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
2962  potl(k,i,j) = ( 289.0_rp - grav / cpdry * cz(k) * geop_sw ) * (0.5_rp-sint) &
2963  + ( 297.5_rp+sign(abs(cz(k)-840.0_rp)**(1.0_rp/3.0_rp),cz(k)-840.0_rp) &
2964  - grav / cpdry * cz(k) * geop_sw ) * (0.5_rp+sint)
2965  else
2966  potl(k,i,j) = 297.5_rp + ( cz(k)-840.0_rp )**(1.0_rp/3.0_rp) &
2967  - grav / cpdry * cz(k) * geop_sw
2968  endif
2969  enddo
2970 
2971  enddo
2972  enddo
2973 
2974  ! make density & pressure profile in dry condition
2975  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
2976  potl(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
2977  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
2978  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
2979  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
2980 
2981  ! calc in moist condition
2982  do j = jsb, jeb
2983  do i = isb, ieb
2984  qv_sfc(i,j) = 9.0e-3_rp ! [kg/kg]
2985 
2986  do k = ks, ke
2987  if ( cz(k) < 820.0_rp ) then ! below initial cloud top
2988  qall = 9.0e-3_rp
2989  elseif( cz(k) <= 860.0_rp ) then ! boundary
2990  sint = sin( pi2 * ( cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
2991  qall = 9.0e-3_rp * (0.5_rp-sint) &
2992  + 1.5e-3_rp * (0.5_rp+sint)
2993  elseif( cz(k) <= 5000.0_rp ) then
2994  qall = 1.5e-3_rp
2995  else
2996  qall = 0.0_rp
2997  endif
2998 
2999  if ( cz(k) <= 600.0_rp ) then
3000  qc(k,i,j) = 0.0_rp
3001  elseif( cz(k) < 820.0_rp ) then ! in the cloud
3002  fact = ( cz(k)-600.0_rp ) / ( 840.0_rp-600.0_rp )
3003  qc(k,i,j) = 0.45e-3_rp * fact
3004  elseif( cz(k) <= 860.0_rp ) then ! boundary
3005  sint = sin( pi2 * ( cz(k)-840.0_rp ) / 20.0_rp ) * 0.5_rp
3006  fact = ( cz(k)-600.0_rp ) / ( 840.0_rp-600.0_rp )
3007  qc(k,i,j) = 0.45e-3_rp * fact * (0.5_rp-sint)
3008  else
3009  qc(k,i,j) = 0.0_rp
3010  endif
3011 
3012  qv(k,i,j) = qall - qc(k,i,j)
3013  enddo
3014 
3015  enddo
3016  enddo
3017 
3018  call hydrometeor_lhv( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3019  temp(:,:,:), lhv(:,:,:) )
3020 
3021  do j = jsb, jeb
3022  do i = isb, ieb
3023  do k = ks, ke
3024  temp(k,i,j) = temp(k,i,j) + lhv(k,i,j) / cpdry * qc(k,i,j)
3025  enddo
3026  enddo
3027  enddo
3028 
3029  ! make density & pressure profile in moist condition
3030  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3031  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3032  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3033  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3034  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3035 
3036  do j = jsb, jeb
3037  do i = isb, ieb
3038  dens( 1:ks-1,i,j) = dens(ks,i,j)
3039  dens(ke+1:ka, i,j) = dens(ke,i,j)
3040  enddo
3041  enddo
3042 
3043  call comm_vars8( dens(:,:,:), 1 )
3044  call comm_wait ( dens(:,:,:), 1 )
3045 
3046  call random_uniform(rndm) ! make random
3047  do j = jsb, jeb
3048  do i = isb, ieb
3049  do k = ks, ke
3050  if ( random_flag == 2 .and. k <= random_limit ) then ! below initial cloud top
3051  momz(k,i,j) = ( 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3052  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3053  else
3054  momz(k,i,j) = 0.0_rp
3055  endif
3056  enddo
3057  enddo
3058  enddo
3059 
3060  call random_uniform(rndm) ! make random
3061  do j = jsb, jeb
3062  do i = isb, ieb
3063  do k = ks, ke
3064  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
3065  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3066  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3067  else
3068  momx(k,i,j) = velx(k,i,j) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3069  endif
3070  enddo
3071  enddo
3072  enddo
3073 
3074  call random_uniform(rndm) ! make random
3075  do j = jsb, jeb
3076  do i = isb, ieb
3077  do k = ks, ke
3078  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
3079  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3080  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3081  else
3082  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3083  endif
3084  enddo
3085  enddo
3086  enddo
3087 
3088  call random_uniform(rndm) ! make random
3089  do j = jsb, jeb
3090  do i = isb, ieb
3091  do k = ks, ke
3092  if ( random_flag == 1 .and. k <= random_limit ) then ! below initial cloud top
3093  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
3094  * dens(k,i,j)
3095  else
3096  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3097  endif
3098  enddo
3099  enddo
3100  enddo
3101 
3102  do j = jsb, jeb
3103  do i = isb, ieb
3104  do k = ks, ke
3105  if ( qc(k,i,j) > 0.0_rp ) then
3106  nc(k,i,j) = 120.e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3107  end if
3108  enddo
3109  enddo
3110  enddo
3111 
3112 #endif
3113  return
3114  end subroutine mkinit_dycoms2_rf01
3115 
3116  !-----------------------------------------------------------------------------
3118  subroutine mkinit_dycoms2_rf02
3119  use scale_atmos_hydrometeor, only: &
3121  implicit none
3122 
3123 #ifndef DRY
3124  real(RP) :: PERTURB_AMP = 0.0_rp
3125  integer :: RANDOM_LIMIT = 5
3126  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
3127  ! 1 -> perturbation for PT
3128  ! 2 -> perturbation for u,v,w
3129 
3130  namelist / param_mkinit_rf02 / &
3131  perturb_amp, &
3132  random_limit, &
3133  random_flag
3134 
3135  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3136  real(RP) :: LHV (ka,ia,ja) ! latent heat of vaporization [J/kg]
3137 
3138  real(RP) :: qall ! QV+QC
3139  real(RP) :: fact
3140  real(RP) :: pi2
3141  real(RP) :: sint
3142 
3143  integer :: ierr
3144  integer :: k, i, j
3145  !---------------------------------------------------------------------------
3146 
3147  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
3148  log_newline
3149  log_info("MKINIT_DYCOMS2_RF02",*) 'Setup initial state'
3150 
3151  if ( atmos_hydrometeor_dry ) then
3152  log_error("MKINIT_DYCOMS2_RF02",*) 'QV is not registered'
3153  call prc_abort
3154  end if
3155 
3156  rewind(io_fid_conf)
3157  read(io_fid_conf,nml=param_mkinit_rf02,iostat=ierr)
3158  if( ierr < 0 ) then !--- missing
3159  log_info("MKINIT_DYCOMS2_RF02",*) 'Not found namelist. Default used.'
3160  elseif( ierr > 0 ) then !--- fatal error
3161  log_error("MKINIT_DYCOMS2_RF02",*) 'Not appropriate names in namelist PARAM_MKINIT_RF02. Check!'
3162  call prc_abort
3163  endif
3164  log_nml(param_mkinit_rf02)
3165 
3166  ! calc in dry condition
3167  call random_uniform(rndm) ! make random
3168  do j = jsb, jeb
3169  do i = isb, ieb
3170 
3171  pres_sfc(i,j) = 1017.8e2_rp ! [Pa]
3172  pott_sfc(i,j) = 288.3_rp ! [K]
3173 
3174  do k = ks, ke
3175  velx(k,i,j) = 3.0_rp + 4.3 * cz(k)*1.e-3_rp
3176  vely(k,i,j) = -9.0_rp + 5.6 * cz(k)*1.e-3_rp
3177 
3178  if ( cz(k) < 775.0_rp ) then ! below initial cloud top
3179  potl(k,i,j) = 288.3_rp ! [K]
3180  else if ( cz(k) <= 815.0_rp ) then
3181  sint = sin( pi2 * (cz(k) - 795.0_rp)/20.0_rp )
3182  potl(k,i,j) = 288.3_rp * (1.0_rp-sint)*0.5_rp &
3183  + ( 295.0_rp+sign(abs(cz(k)-795.0_rp)**(1.0_rp/3.0_rp),cz(k)-795.0_rp) ) &
3184  * (1.0_rp+sint)*0.5_rp
3185  else
3186  potl(k,i,j) = 295.0_rp + ( cz(k)-795.0_rp )**(1.0_rp/3.0_rp)
3187  endif
3188  enddo
3189  enddo
3190  enddo
3191 
3192  ! make density & pressure profile in dry condition
3193  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3194  potl(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3195  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3196  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3197  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3198 
3199  ! calc in moist condition
3200  do j = jsb, jeb
3201  do i = isb, ieb
3202  qv_sfc(i,j) = 9.45e-3_rp
3203 
3204  do k = ks, ke
3205  if ( cz(k) < 775.0_rp ) then ! below initial cloud top
3206  qall = 9.45e-3_rp ! [kg/kg]
3207  else if ( cz(k) <= 815.0_rp ) then
3208  sint = sin( pi2 * (cz(k) - 795.0_rp)/20.0_rp )
3209  qall = 9.45e-3_rp * (1.0_rp-sint)*0.5_rp + &
3210  ( 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-cz(k))/500.0_rp ) ) ) * (1.0_rp+sint)*0.5_rp
3211  else
3212  qall = 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-cz(k))/500.0_rp ) ) ! [kg/kg]
3213  endif
3214 
3215  if( cz(k) < 400.0_rp ) then
3216  qc(k,i,j) = 0.0_rp
3217  elseif( cz(k) < 775.0_rp ) then
3218  fact = ( cz(k)-400.0_rp ) / ( 795.0_rp-400.0_rp )
3219  qc(k,i,j) = 0.65e-3_rp * fact
3220  elseif( cz(k) <= 815.0_rp ) then
3221  sint = sin( pi2 * ( cz(k)-795.0_rp )/20.0_rp )
3222  fact = ( cz(k)-400.0_rp ) / ( 795.0_rp-400.0_rp )
3223  qc(k,i,j) = 0.65e-3_rp * fact * (1.0_rp-sint) * 0.5_rp
3224  else
3225  qc(k,i,j) = 0.0_rp
3226  endif
3227  qv(k,i,j) = qall - qc(k,i,j)
3228  enddo
3229 
3230  enddo
3231  enddo
3232 
3233  call hydrometeor_lhv( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3234  temp(:,:,:), lhv(:,:,:) )
3235 
3236  do j = jsb, jeb
3237  do i = isb, ieb
3238  do k = ks, ke
3239  temp(k,i,j) = temp(k,i,j) + lhv(k,i,j) / cpdry * qc(k,i,j)
3240  enddo
3241  enddo
3242  enddo
3243 
3244  ! make density & pressure profile in moist condition
3245  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3246  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3247  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3248  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3249  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3250 
3251  do j = jsb, jeb
3252  do i = isb, ieb
3253  dens( 1:ks-1,i,j) = dens(ks,i,j)
3254  dens(ke+1:ka, i,j) = dens(ke,i,j)
3255  enddo
3256  enddo
3257 
3258  call comm_vars8( dens(:,:,:), 1 )
3259  call comm_wait ( dens(:,:,:), 1 )
3260 
3261  call random_uniform(rndm) ! make random
3262  do j = jsb, jeb
3263  do i = isb, ieb
3264  do k = ks, ke
3265  if( random_flag == 2 .and. k <= random_limit ) then
3266  momz(k,i,j) = ( 0.0_rp + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3267  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3268  else
3269  momz(k,i,j) = 0.0_rp
3270  endif
3271  enddo
3272  enddo
3273  enddo
3274 
3275  call random_uniform(rndm) ! make random
3276  do j = jsb, jeb
3277  do i = isb, ieb
3278  do k = ks, ke
3279  if( random_flag == 2 .and. k <= random_limit ) then
3280  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3281  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3282  else
3283  momx(k,i,j) = ( velx(k,i,j) ) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3284  endif
3285  enddo
3286  enddo
3287  enddo
3288 
3289  call random_uniform(rndm) ! make random
3290  do j = jsb, jeb
3291  do i = isb, ieb
3292  do k = ks, ke
3293  if( random_flag == 2 .and. k <= random_limit ) then
3294  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3295  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3296  else
3297  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3298  endif
3299  enddo
3300  enddo
3301  enddo
3302 
3303  call random_uniform(rndm) ! make random
3304  do j = jsb, jeb
3305  do i = isb, ieb
3306  do k = ks, ke
3307  if( random_flag == 1 .and. k <= random_limit ) then
3308  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3309  * dens(k,i,j)
3310  else
3311  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3312  endif
3313  enddo
3314  enddo
3315  enddo
3316 
3317  do j = jsb, jeb
3318  do i = isb, ieb
3319  do k = ks, ke
3320  if ( qc(k,i,j) > 0.0_rp ) then
3321  nc(k,i,j) = 55.0e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3322  endif
3323  enddo
3324  enddo
3325  enddo
3326 
3327 #endif
3328  return
3329  end subroutine mkinit_dycoms2_rf02
3330 
3331  !-----------------------------------------------------------------------------
3333  subroutine mkinit_dycoms2_rf02_dns
3334  use scale_atmos_hydrometeor, only: &
3336  implicit none
3337 
3338 #ifndef DRY
3339  real(RP) :: ZB = 750.0_rp ! domain bottom
3340 ! real(RP) :: ZT = 900.0_RP ! domain top
3341  real(RP) :: CONST_U = 0.0_rp
3342  real(RP) :: CONST_V = 0.0_rp
3343  real(RP) :: PRES_ZB = 93060.0_rp
3344  real(RP) :: PERTURB_AMP = 0.0_rp
3345  integer :: RANDOM_LIMIT = 5
3346  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
3347  ! 1 -> perturbation for PT
3348  ! 2 -> perturbation for u,v,w
3349 
3350  namelist / param_mkinit_rf02_dns / &
3351  zb, const_u, const_v,pres_zb,&
3352  perturb_amp, &
3353  random_limit, &
3354  random_flag
3355 
3356  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3357  real(RP) :: LHV (ka,ia,ja) ! latent heat of vaporization [J/kg]
3358 
3359  real(RP) :: qall ! QV+QC
3360  real(RP) :: fact
3361  real(RP) :: pi2
3362  real(RP) :: RovCP
3363 
3364  integer :: ierr
3365  integer :: k, i, j
3366  !---------------------------------------------------------------------------
3367 
3368  pi2 = atan(1.0_rp) * 2.0_rp ! pi/2
3369 
3370  log_newline
3371  log_info("MKINIT_DYCOMS2_RF02_DNS",*) 'Setup initial state'
3372 
3373  if ( atmos_hydrometeor_dry ) then
3374  log_error("MKINIT_DYCOMS2_RF02_DNS",*) 'QV is not registered'
3375  call prc_abort
3376  end if
3377 
3378  rewind(io_fid_conf)
3379  read(io_fid_conf,nml=param_mkinit_rf02_dns,iostat=ierr)
3380  if( ierr < 0 ) then !--- missing
3381  log_info("MKINIT_DYCOMS2_RF02_DNS",*) 'Not found namelist. Default used.'
3382  elseif( ierr > 0 ) then !--- fatal error
3383  log_error("MKINIT_DYCOMS2_RF02_DNS",*) 'Not appropriate names in namelist PARAM_MKINIT_RF02_DNS. Check!'
3384  call prc_abort
3385  endif
3386  log_nml(param_mkinit_rf02_dns)
3387 
3388  ! calc in dry condition
3389  call random_uniform(rndm) ! make random
3390  do j = jsb, jeb
3391  do i = isb, ieb
3392 
3393  pres_sfc(i,j) = pres_zb
3394 ! pott_sfc(i,j) = 288.3_RP ! [K]
3395 ! qv_sfc (i,j) = 9.45E-3_RP
3396 
3397  do k = ks, ke
3398 
3399  velx(k,i,j) = const_u
3400  vely(k,i,j) = const_v
3401 
3402 ! if ( ZB+CZ(k) < 775.0_RP ) then ! below initial cloud top
3403  if ( zb+cz(k) <= 795.0_rp ) then ! below initial cloud top
3404  potl(k,i,j) = 288.3_rp ! [K]
3405  qall = 9.45e-3_rp ! [kg/kg]
3406 ! necessary?
3407 ! else if ( CZ(k) <= 815.0_RP ) then
3408 ! sint = sin( pi2 * (CZ(k) - 795.0_RP)/20.0_RP )
3409 ! potl(k,i,j) = 288.3_RP * (1.0_RP-sint)*0.5_RP + &
3410 ! ( 295.0_RP+sign(abs(CZ(k)-795.0_RP)**(1.0_RP/3.0_RP),CZ(k)-795.0_RP) ) * (1.0_RP+sint)*0.5_RP
3411 ! qall = 9.45E-3_RP * (1.0_RP-sint)*0.5_RP + &
3412 ! ( 5.E-3_RP - 3.E-3_RP * ( 1.0_RP - exp( (795.0_RP-CZ(k))/500.0_RP ) ) ) * (1.0_RP+sint)*0.5_RP
3413  else
3414  potl(k,i,j) = 295.0_rp + ( zb+cz(k)-795.0_rp )**(1.0_rp/3.0_rp)
3415  qall = 5.e-3_rp - 3.e-3_rp * ( 1.0_rp - exp( (795.0_rp-(zb+cz(k)))/500.0_rp ) ) ! [kg/kg]
3416  endif
3417 
3418  if( zb+cz(k) < 400.0_rp ) then
3419  qc(k,i,j) = 0.0_rp
3420  elseif( zb+cz(k) <= 795.0_rp ) then
3421  fact = ( (zb+cz(k))-400.0_rp ) / ( 795.0_rp-400.0_rp )
3422  qc(k,i,j) = 0.8e-3_rp * fact
3423  else
3424  qc(k,i,j) = 0.0_rp
3425  endif
3426  qv(k,i,j) = qall - qc(k,i,j)
3427 
3428  !if(i==is.and.j==js)LOG_INFO("MKINIT_DYCOMS2_RF02_DNS",*)'chkk',k,cz(k)+zb,qc(k,i,j),qv(k,i,j)
3429  enddo
3430  enddo
3431  enddo
3432 
3433  !LOG_INFO("MKINIT_DYCOMS2_RF02_DNS",*)'chk3',ks,ke
3434  ! extrapolation (temtative)
3435  pott_sfc(:,:) = potl(ks,:,:)-0.5*(potl(ks+1,:,:)-potl(ks,:,:))
3436  qv_sfc(:,:) = qv(ks,:,:)-0.5*(qv(ks+1,:,:)-qv(ks,:,:))
3437  qc_sfc(:,:) = qc(ks,:,:)-0.5*(qc(ks+1,:,:)-qc(ks,:,:))
3438 
3439  ! make density & pressure profile in moist condition
3440  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3441  potl(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3442  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3443  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3444  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3445 
3446  call hydrometeor_lhv( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3447  temp(:,:,:), lhv(:,:,:) )
3448 
3449  rovcp = rdry / cpdry
3450  do j = jsb, jeb
3451  do i = isb, ieb
3452  do k = ks, ke
3453  pott(k,i,j) = potl(k,i,j) + lhv(k,i,j) / cpdry * qc(k,i,j) * ( p00/pres(k,i,j) )**rovcp
3454  enddo
3455  enddo
3456  enddo
3457 
3458  ! make density & pressure profile in moist condition
3459  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3460  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3461  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3462  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3463  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3464 
3465  do j = jsb, jeb
3466  do i = isb, ieb
3467  dens( 1:ks-1,i,j) = dens(ks,i,j)
3468  dens(ke+1:ka, i,j) = dens(ke,i,j)
3469  enddo
3470  enddo
3471 
3472  call comm_vars8( dens(:,:,:), 1 )
3473  call comm_wait ( dens(:,:,:), 1 )
3474 
3475  call random_uniform(rndm) ! make random
3476  do j = jsb, jeb
3477  do i = isb, ieb
3478  do k = ks, ke
3479  if( random_flag == 2 .and. k <= random_limit ) then
3480  momz(k,i,j) = ( 0.0_rp + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3481  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
3482  else
3483  momz(k,i,j) = 0.0_rp
3484  endif
3485  enddo
3486  enddo
3487  enddo
3488 
3489  !LOG_INFO("MKINIT_DYCOMS2_RF02_DNS",*)'chk8'
3490  call random_uniform(rndm) ! make random
3491  do j = jsb, jeb
3492  do i = isb, ieb
3493  do k = ks, ke
3494  if( random_flag == 2 .and. k <= random_limit ) then
3495  momx(k,i,j) = ( velx(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3496  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3497  else
3498  momx(k,i,j) = ( velx(k,i,j) ) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3499  endif
3500  enddo
3501  enddo
3502  enddo
3503  !LOG_INFO("MKINIT_DYCOMS2_RF02_DNS",*)'chk9'
3504 
3505  call random_uniform(rndm) ! make random
3506  do j = jsb, jeb
3507  do i = isb, ieb
3508  do k = ks, ke
3509  if( random_flag == 2 .and. k <= random_limit ) then
3510  momy(k,i,j) = ( vely(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3511  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3512  else
3513  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3514  endif
3515  enddo
3516  enddo
3517  enddo
3518 
3519  call random_uniform(rndm) ! make random
3520  do j = jsb, jeb
3521  do i = isb, ieb
3522  do k = ks, ke
3523  if( random_flag == 1 .and. k <= random_limit ) then
3524  rhot(k,i,j) = ( pott(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp ) &
3525  * dens(k,i,j)
3526  else
3527  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3528  endif
3529  enddo
3530  enddo
3531  enddo
3532 
3533  do j = jsb, jeb
3534  do i = isb, ieb
3535  do k = ks, ke
3536  if ( qc(k,i,j) > 0.0_rp ) then
3537  nc(k,i,j) = 55.0e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3538  endif
3539  enddo
3540  enddo
3541  enddo
3542 
3543 #endif
3544  return
3545  end subroutine mkinit_dycoms2_rf02_dns
3546 
3547  !-----------------------------------------------------------------------------
3549  subroutine mkinit_rico
3550  use scale_atmos_hydrometeor, only: &
3552  implicit none
3553 
3554 #ifndef DRY
3555  real(RP):: PERTURB_AMP_PT = 0.1_rp
3556  real(RP):: PERTURB_AMP_QV = 2.5e-5_rp
3557 
3558  namelist / param_mkinit_rico / &
3559  perturb_amp_pt, &
3560  perturb_amp_qv
3561 
3562  real(RP) :: LHV (ka,ia,ja) ! latent heat of vaporization [J/kg]
3563  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3564  real(RP) :: qall ! QV+QC
3565  real(RP) :: fact
3566 
3567  integer :: ierr
3568  integer :: k, i, j
3569  !---------------------------------------------------------------------------
3570 
3571  log_newline
3572  log_info("MKINIT_RICO",*) 'Setup initial state'
3573 
3574  if ( atmos_hydrometeor_dry ) then
3575  log_error("MKINIT_RICO",*) 'QV is not registered'
3576  call prc_abort
3577  end if
3578 
3579  rewind(io_fid_conf)
3580  read(io_fid_conf,nml=param_mkinit_rico,iostat=ierr)
3581  if( ierr < 0 ) then !--- missing
3582  log_info("MKINIT_RICO",*) 'Not found namelist. Default used.'
3583  elseif( ierr > 0 ) then !--- fatal error
3584  log_error("MKINIT_RICO",*) 'Not appropriate names in namelist PARAM_MKINIT_RICO. Check!'
3585  call prc_abort
3586  endif
3587  log_nml(param_mkinit_rico)
3588 
3589  ! calc in moist condition
3590  do j = jsb, jeb
3591  do i = isb, ieb
3592 
3593  pres_sfc(i,j) = 1015.4e2_rp ! [Pa]
3594  pott_sfc(i,j) = 297.9_rp
3595 
3596  do k = ks, ke
3597  !--- potential temperature
3598  if ( cz(k) < 740.0_rp ) then ! below initial cloud top
3599  potl(k,i,j) = 297.9_rp
3600  else
3601  fact = ( cz(k)-740.0_rp ) * ( 317.0_rp-297.9_rp ) / ( 4000.0_rp-740.0_rp )
3602  potl(k,i,j) = 297.9_rp + fact
3603  endif
3604 
3605  !--- horizontal wind velocity
3606  if ( cz(k) <= 4000.0_rp ) then ! below initial cloud top
3607  fact = ( cz(k)-0.0_rp ) * ( -1.9_rp+9.9_rp ) / ( 4000.0_rp-0.0_rp )
3608  velx(k,i,j) = -9.9_rp + fact
3609  vely(k,i,j) = -3.8_rp
3610  else
3611  velx(k,i,j) = -1.9_rp
3612  vely(k,i,j) = -3.8_rp
3613  endif
3614  enddo
3615 
3616  enddo
3617  enddo
3618 
3619  ! make density & pressure profile in moist condition
3620  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3621  potl(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3622  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3623  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3624  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3625 
3626 
3627  do j = jsb, jeb
3628  do i = isb, ieb
3629  qv_sfc(i,j) = 16.0e-3_rp ! [kg/kg]
3630 
3631  do k = ks, ke
3632  !--- mixing ratio of vapor
3633  if ( cz(k) <= 740.0_rp ) then ! below initial cloud top
3634  fact = ( cz(k)-0.0_rp ) * ( 13.8e-3_rp-16.0e-3_rp ) / ( 740.0_rp-0.0_rp )
3635  qall = 16.0e-3_rp + fact
3636  elseif ( cz(k) <= 3260.0_rp ) then ! boundary
3637  fact = ( cz(k)-740.0_rp ) * ( 2.4e-3_rp-13.8e-3_rp ) / ( 3260.0_rp-740.0_rp )
3638  qall = 13.8e-3_rp + fact
3639  elseif( cz(k) <= 4000.0_rp ) then
3640  fact = ( cz(k)-3260.0_rp ) * ( 1.8e-3_rp-2.4e-3_rp ) / ( 4000.0_rp-3260.0_rp )
3641  qall = 2.4e-3_rp + fact
3642  else
3643  qall = 0.0_rp
3644  endif
3645 
3646  qv(k,i,j) = qall - qc(k,i,j)
3647  enddo
3648 
3649  enddo
3650  enddo
3651 
3652  call hydrometeor_lhv( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3653  temp(:,:,:), lhv(:,:,:) )
3654 
3655  do j = jsb, jeb
3656  do i = isb, ieb
3657  do k = ks, ke
3658  temp(k,i,j) = temp(k,i,j) + lhv(k,i,j) / cpdry * qc(k,i,j)
3659  enddo
3660  enddo
3661  enddo
3662 
3663  ! make density & pressure profile in moist condition
3664  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3665  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3666  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3667  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3668  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3669 
3670 
3671  do j = jsb, jeb
3672  do i = isb, ieb
3673  dens( 1:ks-1,i,j) = dens(ks,i,j)
3674  dens(ke+1:ka ,i,j) = dens(ke,i,j)
3675  enddo
3676  enddo
3677 
3678  call comm_vars8( dens(:,:,:), 1 )
3679  call comm_wait ( dens(:,:,:), 1 )
3680 
3681  do j = jsb, jeb
3682  do i = isb, ieb
3683  do k = ks, ke
3684  momz(k,i,j) = 0.0_rp
3685  enddo
3686  enddo
3687  enddo
3688 
3689  do j = jsb, jeb
3690  do i = isb, ieb
3691  do k = ks, ke
3692  momx(k,i,j) = velx(k,i,j) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3693  enddo
3694  enddo
3695  enddo
3696 
3697  do j = jsb, jeb
3698  do i = isb, ieb
3699  do k = ks, ke
3700  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3701  enddo
3702  enddo
3703  enddo
3704 
3705  call random_uniform(rndm) ! make random
3706  do j = jsb, jeb
3707  do i = isb, ieb
3708  do k = ks, ke
3709  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)
3710  enddo
3711  enddo
3712  enddo
3713 
3714  call random_uniform(rndm) ! make random
3715  do j = jsb, jeb
3716  do i = isb, ieb
3717  do k = ks, ke
3718  qv(k,i,j) = qv(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp_qv
3719  enddo
3720  enddo
3721  enddo
3722 
3723  do j = jsb, jeb
3724  do i = isb, ieb
3725  do k = ks, ke
3726  if ( qc(k,i,j) > 0.0_rp ) then
3727  nc(k,i,j) = 70.e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3728  endif
3729  enddo
3730  enddo
3731  enddo
3732 
3733 #endif
3734  return
3735  end subroutine mkinit_rico
3736 
3737  !-----------------------------------------------------------------------------
3739  subroutine mkinit_bomex
3740  use scale_const, only: &
3741  rdry => const_rdry, &
3742  rvap => const_rvap, &
3743  cpdry => const_cpdry, &
3744  cpvap => const_cpvap, &
3745  cl => const_cl
3746  use scale_atmos_hydrometeor, only: &
3748  implicit none
3749 
3750 #ifndef DRY
3751  real(RP):: PERTURB_AMP_PT = 0.1_rp
3752  real(RP):: PERTURB_AMP_QV = 2.5e-5_rp
3753 
3754  namelist / param_mkinit_bomex / &
3755  perturb_amp_pt, &
3756  perturb_amp_qv
3757 
3758  real(RP) :: LHV (ka,ia,ja) ! latent heat of vaporization [J/kg]
3759  real(RP) :: potl(ka,ia,ja) ! liquid potential temperature
3760  real(RP) :: qall ! QV+QC
3761  real(RP) :: fact
3762 
3763  real(RP) :: qdry, Rtot, CPtot
3764 
3765  integer :: ierr
3766  integer :: k, i, j
3767  !---------------------------------------------------------------------------
3768 
3769  log_newline
3770  log_info("MKINIT_BOMEX",*) 'Setup initial state'
3771 
3772  if ( atmos_hydrometeor_dry ) then
3773  log_error("MKINIT_BOMEX",*) 'QV is not registered'
3774  call prc_abort
3775  end if
3776 
3777  rewind(io_fid_conf)
3778  read(io_fid_conf,nml=param_mkinit_bomex,iostat=ierr)
3779  if( ierr < 0 ) then !--- missing
3780  log_info("MKINIT_BOMEX",*) 'Not found namelist. Default used.'
3781  elseif( ierr > 0 ) then !--- fatal error
3782  log_error("MKINIT_BOMEX",*) 'Not appropriate names in namelist PARAM_MKINIT_BOMEX. Check!'
3783  call prc_abort
3784  endif
3785  log_nml(param_mkinit_bomex)
3786 
3787  ! calc in moist condition
3788  do j = jsb, jeb
3789  do i = isb, ieb
3790 
3791  pres_sfc(i,j) = 1015.e2_rp ! [Pa]
3792  pott_sfc(i,j) = 299.1_rp
3793 
3794  do k = ks, ke
3795  !--- potential temperature
3796  if ( cz(k) < 520.0_rp ) then ! below initial cloud top
3797  potl(k,i,j) = 298.7_rp
3798  elseif( cz(k) < 1480.0_rp ) then
3799  fact = ( cz(k)-520.0_rp ) * ( 302.4_rp-298.7_rp ) / ( 1480.0_rp-520.0_rp )
3800  potl(k,i,j) = 298.7_rp + fact
3801  elseif( cz(k) < 2000.0_rp ) then
3802  fact = ( cz(k)-1480.0_rp ) * ( 308.2_rp-302.4_rp ) / ( 2000.0_rp-1480.0_rp )
3803  potl(k,i,j) = 302.4_rp + fact
3804  else
3805  fact = ( cz(k)-2000.0_rp ) * 3.65e-3_rp
3806  potl(k,i,j) = 308.2_rp + fact
3807  endif
3808 
3809  !--- horizontal wind velocity
3810  if ( cz(k) <= 700.0_rp ) then ! below initial cloud top
3811  velx(k,i,j) = -8.75_rp
3812  vely(k,i,j) = 0.0_rp
3813  else
3814  fact = 1.8e-3_rp * ( cz(k)-700.0_rp )
3815  velx(k,i,j) = -8.75_rp + fact
3816  vely(k,i,j) = 0.0_rp
3817  endif
3818  enddo
3819 
3820  enddo
3821  enddo
3822 
3823  ! make density & pressure profile in moist condition
3824  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3825  potl(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3826  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3827  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3828  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3829 
3830 
3831  do j = jsb, jeb
3832  do i = isb, ieb
3833  qv_sfc(i,j) = 22.45e-3_rp ! [kg/kg]
3834 
3835  do k = ks, ke
3836  !--- mixing ratio of vapor
3837  if ( cz(k) <= 520.0_rp ) then ! below initial cloud top
3838  fact = ( cz(k)-0.0_rp ) * ( 16.3e-3_rp-17.0e-3_rp ) / ( 520.0_rp-0.0_rp )
3839  qall = 17.0e-3_rp + fact
3840  elseif ( cz(k) <= 1480.0_rp ) then ! boundary
3841  fact = ( cz(k)-520.0_rp ) * ( 10.7e-3_rp-16.3e-3_rp ) / ( 1480.0_rp-520.0_rp )
3842  qall = 16.3e-3_rp + fact
3843  elseif( cz(k) <= 2000.0_rp ) then
3844  fact = ( cz(k)-1480.0_rp ) * ( 4.2e-3_rp-10.7e-3_rp ) / ( 2000.0_rp-1480.0_rp )
3845  qall = 10.7e-3_rp + fact
3846  else
3847  fact = ( cz(k)-2000.0_rp ) * ( -1.2e-6_rp )
3848  qall = 4.2e-3_rp + fact
3849  endif
3850 
3851  qv(k,i,j) = qall - qc(k,i,j)
3852  enddo
3853 
3854  enddo
3855  enddo
3856 
3857  call hydrometeor_lhv( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3858  temp(:,:,:), lhv(:,:,:) )
3859 
3860  do j = jsb, jeb
3861  do i = isb, ieb
3862  do k = ks, ke
3863  qdry = 1.0_rp - qv(k,i,j) - qc(k,i,j)
3864  rtot = rdry * qdry + rvap * qv(k,i,j)
3865  cptot = cpdry * qdry + cpvap * qv(k,i,j) + cl * qc(k,i,j)
3866  pott(k,i,j) = ( temp(k,i,j) + lhv(k,i,j) / cpdry * qc(k,i,j) ) * ( p00 / pres(k,i,j) )**(rtot/cptot)
3867  enddo
3868  enddo
3869  enddo
3870 
3871  ! make density & pressure profile in moist condition
3872  call hydrostatic_buildrho( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
3873  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
3874  pres_sfc(:,:), pott_sfc(:,:), qv_sfc(:,:), qc_sfc(:,:), & ! [IN]
3875  real_cz(:,:,:), real_fz(:,:,:), area(:,:), & ! [IN]
3876  dens(:,:,:), temp(:,:,:), pres(:,:,:), temp_sfc(:,:) ) ! [OUT]
3877 
3878 
3879  do j = jsb, jeb
3880  do i = isb, ieb
3881  dens( 1:ks-1,i,j) = dens(ks,i,j)
3882  dens(ke+1:ka ,i,j) = dens(ke,i,j)
3883  enddo
3884  enddo
3885 
3886  call comm_vars8( dens(:,:,:), 1 )
3887  call comm_wait ( dens(:,:,:), 1 )
3888 
3889  do j = jsb, jeb
3890  do i = isb, ieb
3891  do k = ks, ke
3892  momz(k,i,j) = 0.0_rp
3893  enddo
3894  enddo
3895  enddo
3896 
3897  do j = jsb, jeb
3898  do i = isb, ieb
3899  do k = ks, ke
3900  momx(k,i,j) = velx(k,i,j) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
3901  enddo
3902  enddo
3903  enddo
3904 
3905  do j = jsb, jeb
3906  do i = isb, ieb
3907  do k = ks, ke
3908  momy(k,i,j) = vely(k,i,j) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
3909  enddo
3910  enddo
3911  enddo
3912 
3913  call random_uniform(rndm) ! make random
3914  do j = jsb, jeb
3915  do i = isb, ieb
3916  do k = ks, ke
3917  if( cz(k) <= 1600.0_rp ) then !--- lowest 40 model layer when dz=40m
3918  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)
3919  else
3920  rhot(k,i,j) = pott(k,i,j) * dens(k,i,j)
3921  endif
3922  enddo
3923  enddo
3924  enddo
3925 
3926  call random_uniform(rndm) ! make random
3927  do j = jsb, jeb
3928  do i = isb, ieb
3929  do k = ks, ke
3930  if( cz(k) <= 1600.0_rp ) then !--- lowest 40 model layer when dz=40m
3931  qv(k,i,j) = qv(k,i,j) + 2.0_rp * ( rndm(k,i,j)-0.50_rp ) * perturb_amp_qv
3932  endif
3933  enddo
3934  enddo
3935  enddo
3936 
3937  do j = jsb, jeb
3938  do i = isb, ieb
3939  do k = ks, ke
3940  if ( qc(k,i,j) > 0.0_rp ) then
3941  nc(k,i,j) = 70.e6_rp / dens(k,i,j) ! [number/m3] / [kg/m3]
3942  endif
3943  enddo
3944  enddo
3945  enddo
3946 
3947 #endif
3948  return
3949  end subroutine mkinit_bomex
3950 
3951  !-----------------------------------------------------------------------------
3953  subroutine mkinit_oceancouple
3954  implicit none
3955 
3956  log_newline
3957  log_info("MKINIT_oceancouple",*) 'Setup initial state'
3958 
3959  call flux_setup
3960 
3961  call ocean_setup
3962 
3963  return
3964  end subroutine mkinit_oceancouple
3965 
3966  !-----------------------------------------------------------------------------
3968  subroutine mkinit_landcouple
3969  implicit none
3970 
3971  log_newline
3972  log_info("MKINIT_landcouple",*) 'Setup initial state'
3973 
3974  call flux_setup
3975 
3976  call land_setup
3977 
3978  return
3979  end subroutine mkinit_landcouple
3980 
3981  !-----------------------------------------------------------------------------
3983  subroutine mkinit_urbancouple
3984  implicit none
3985 
3986  log_newline
3987  log_info("MKINIT_urbancouple",*) 'Setup initial state'
3988 
3989  call flux_setup
3990 
3991  call urban_setup
3992 
3993  return
3994  end subroutine mkinit_urbancouple
3995 
3996  !-----------------------------------------------------------------------------
3998  subroutine mkinit_seabreeze
3999  use scale_prc_cartesc, only: &
4000  prc_num_x
4001  use scale_landuse, only: &
4004  landuse_fillhalo, &
4006  use scale_atmos_grid_cartesc, only: &
4007  domain_center_x => atmos_grid_cartesc_domain_center_x
4008  implicit none
4009 
4010  real(RP) :: LAND_SIZE
4011 
4012  namelist / param_mkinit_seabreeze / &
4013  land_size
4014 
4015  integer :: ierr
4016  integer :: i, j
4017  !---------------------------------------------------------------------------
4018 
4019  log_newline
4020  log_info("MKINIT_seabreeze",*) 'Setup initial state'
4021 
4022  land_size = 0.0_rp
4023 
4024  !--- read namelist
4025  rewind(io_fid_conf)
4026  read(io_fid_conf,nml=param_mkinit_seabreeze,iostat=ierr)
4027 
4028  if( ierr < 0 ) then !--- missing
4029  log_info("MKINIT_seabreeze",*) 'Not found namelist. Default used.'
4030  elseif( ierr > 0 ) then !--- fatal error
4031  log_error("MKINIT_seabreeze",*) 'Not appropriate names in namelist PARAM_MKINIT_SEABREEZE. Check!'
4032  call prc_abort
4033  endif
4034  log_nml(param_mkinit_seabreeze)
4035 
4036  call flux_setup
4037 
4038  call land_setup
4039 
4040  call ocean_setup
4041 
4042  ! make landuse conditions
4043  do j = jsb, jeb
4044  do i = isb, ieb
4045  if ( abs( cx(i) - domain_center_x ) < land_size ) then
4046  landuse_frac_land(i,j) = 1.0_rp
4047  else
4048  landuse_frac_land(i,j) = 0.0_rp
4049  endif
4050  enddo
4051  enddo
4052 
4053  ! calculate landuse factors
4054  call landuse_fillhalo( fill_bnd=.true. )
4055  call landuse_calc_fact
4056 
4057  ! output landuse file
4058  call landuse_write
4059 
4060  return
4061  end subroutine mkinit_seabreeze
4062 
4063  !-----------------------------------------------------------------------------
4065  subroutine mkinit_heatisland
4066  use scale_prc_cartesc, only: &
4067  prc_num_x
4068  use scale_landuse, only: &
4072  landuse_fillhalo, &
4074  implicit none
4075 
4076  real(RP) :: dist
4077 
4078  integer :: i, j
4079  !---------------------------------------------------------------------------
4080 
4081  log_newline
4082  log_info("MKINIT_heatisland",*) 'Setup initial state'
4083 
4084  call flux_setup
4085 
4086  call land_setup
4087 
4088  call urban_setup
4089 
4090  ! 1/9 size of domain
4091  dist = ( cxg(imax*prc_num_x) - cxg(1) ) / 9.0_rp
4092 
4093  ! make landuse conditions
4094  do j = jsb, jeb
4095  do i = isb, ieb
4096  if ( cx(i) >= dist * 4.0_rp &
4097  .AND. cx(i) < dist * 5.0_rp ) then
4098  landuse_frac_land(i,j) = 1.0_rp
4099  landuse_frac_urban(i,j) = 1.0_rp
4100  else
4101  landuse_frac_land(i,j) = 1.0_rp
4102  landuse_frac_urban(i,j) = 0.0_rp
4103  endif
4104  enddo
4105  enddo
4106 
4107  ! calculate landuse factors
4108  call landuse_fillhalo( fill_bnd=.true. )
4109  call landuse_calc_fact
4110 
4111  ! output landuse file
4112  call landuse_write
4113 
4114  return
4115  end subroutine mkinit_heatisland
4116 
4117  !-----------------------------------------------------------------------------
4119  subroutine mkinit_grayzone
4120  use scale_atmos_hydrometeor, only: &
4122  implicit none
4123 
4124  real(RP) :: RHO(ka)
4125  real(RP) :: VELX(ka)
4126  real(RP) :: VELY(ka)
4127  real(RP) :: POTT(ka)
4128  real(RP) :: QV1D(ka)
4129 
4130  real(RP) :: PERTURB_AMP = 0.0_rp
4131  integer :: RANDOM_LIMIT = 0
4132  integer :: RANDOM_FLAG = 0 ! 0 -> no perturbation
4133  ! 1 -> petrurbation for pt
4134  ! 2 -> perturbation for u, v, w
4135 
4136  namelist / param_mkinit_grayzone / &
4137  perturb_amp, &
4138  random_limit, &
4139  random_flag
4140 
4141  integer :: ierr
4142  integer :: k, i, j
4143  !---------------------------------------------------------------------------
4144 
4145  log_newline
4146  log_info("MKINIT_grayzone",*) 'Setup initial state'
4147 
4148  if ( atmos_hydrometeor_dry ) then
4149  log_error("MKINIT_grayzone",*) 'QV is not registered'
4150  call prc_abort
4151  end if
4152 
4153  !--- read namelist
4154  rewind(io_fid_conf)
4155  read(io_fid_conf,nml=param_mkinit_grayzone,iostat=ierr)
4156 
4157  if( ierr < 0 ) then !--- missing
4158  log_info("MKINIT_grayzone",*) 'Not found namelist. Default used.'
4159  elseif( ierr > 0 ) then !--- fatal error
4160  log_error("MKINIT_grayzone",*) 'Not appropriate names in namelist PARAM_MKINIT_GRAYZONE. Check!'
4161  call prc_abort
4162  endif
4163  log_nml(param_mkinit_grayzone)
4164 
4165  call read_sounding( rho, velx, vely, pott, qv1d ) ! (out)
4166 
4167 ! do j = JS, JE
4168 ! do i = IS, IE
4169  do j = 1, ja
4170  do i = 1, ia
4171  do k = ks, ke
4172  dens(k,i,j) = rho(k)
4173 ! MOMZ(k,i,j) = 0.0_RP
4174 ! MOMX(k,i,j) = RHO(k) * VELX(k)
4175 ! MOMY(k,i,j) = RHO(k) * VELY(k)
4176 
4177 ! RHOT(k,i,j) = RHO(k) * POTT(k)
4178  qv(k,i,j) = qv1d(k)
4179  enddo
4180  enddo
4181  enddo
4182 
4183  do j = jsb, jeb
4184  do i = isb, ieb
4185  dens( 1:ks-1,i,j) = dens(ks,i,j)
4186  dens(ke+1:ka, i,j) = dens(ke,i,j)
4187  enddo
4188  enddo
4189 
4190  call random_uniform(rndm) ! make random
4191  do j = jsb, jeb
4192  do i = isb, ieb
4193  do k = ks, ke
4194  if ( random_flag == 2 .and. k <= random_limit ) then ! below initial cloud top
4195  momz(k,i,j) = ( 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4196  * 0.5_rp * ( dens(k+1,i,j) + dens(k,i,j) )
4197  else
4198  momz(k,i,j) = 0.0_rp
4199  endif
4200  enddo
4201  enddo
4202  enddo
4203 
4204  call random_uniform(rndm) ! make random
4205  do j = jsb, jeb
4206  do i = isb, ieb
4207  do k = ks, ke
4208  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
4209  momx(k,i,j) = ( velx(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4210  * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
4211  else
4212  momx(k,i,j) = velx(k) * 0.5_rp * ( dens(k,i+1,j) + dens(k,i,j) )
4213  endif
4214  enddo
4215  enddo
4216  enddo
4217 
4218  call random_uniform(rndm) ! make random
4219  do j = jsb, jeb
4220  do i = isb, ieb
4221  do k = ks, ke
4222  if ( random_flag == 2 .AND. k <= random_limit ) then ! below initial cloud top
4223  momy(k,i,j) = ( vely(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4224  * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
4225  else
4226  momy(k,i,j) = vely(k) * 0.5_rp * ( dens(k,i,j+1) + dens(k,i,j) )
4227  endif
4228  enddo
4229  enddo
4230  enddo
4231 
4232  call random_uniform(rndm) ! make random
4233  do j = jsb, jeb
4234  do i = isb, ieb
4235  do k = ks, ke
4236  if ( random_flag == 1 .and. k <= random_limit ) then ! below initial cloud top
4237  rhot(k,i,j) = ( pott(k) + 2.0_rp * ( rndm(k,i,j)-0.5_rp ) * perturb_amp ) &
4238  * dens(k,i,j)
4239  else
4240  rhot(k,i,j) = pott(k) * dens(k,i,j)
4241  endif
4242  enddo
4243  enddo
4244  enddo
4245 
4246  return
4247  end subroutine mkinit_grayzone
4248 
4249  !-----------------------------------------------------------------------------
4251  subroutine mkinit_boxaero
4252  use scale_const, only: &
4253  rdry => const_rdry, &
4254  rvap => const_rvap, &
4255  cvdry => const_cvdry, &
4256  cvvap => const_cvvap, &
4257  cpdry => const_cpdry, &
4258  cpvap => const_cpvap
4259  use scale_atmos_hydrometeor, only: &
4261  use scale_atmos_thermodyn, only: &
4262  atmos_thermodyn_rhot2temp_pres
4263  use mod_atmos_admin, only: &
4265  implicit none
4266 
4267  real(RP) :: init_dens = 1.12_rp ![kg/m3]
4268  real(RP) :: init_temp = 298.18_rp ![K]
4269  real(RP) :: init_pres = 1.e+5_rp ![Pa]
4270  real(RP) :: init_ssliq = 0.01_rp ![%]
4271 
4272  namelist / param_mkinit_boxaero / &
4273  init_dens, &
4274  init_temp, &
4275  init_pres, &
4276  init_ssliq
4277 
4278  real(RP) :: rtot (ka,ia,ja)
4279  real(RP) :: cvtot(ka,ia,ja)
4280  real(RP) :: cptot(ka,ia,ja)
4281  real(RP) :: qdry
4282  real(RP) :: psat, qsat
4283  integer :: i, j, k, ierr
4284  !---------------------------------------------------------------------------
4285 
4286  if ( atmos_phy_ae_type /= 'KAJINO13' ) then
4287  log_info("MKINIT_boxaero",*) 'For [Box model of aerosol],'
4288  log_info("MKINIT_boxaero",*) 'ATMOS_PHY_AE_TYPE should be KAJINO13. Stop! ', trim(atmos_phy_ae_type)
4289  call prc_abort
4290  endif
4291 
4292  if ( atmos_hydrometeor_dry ) then
4293  log_error("MKINIT_boxaero",*) 'QV is not registered'
4294  call prc_abort
4295  end if
4296 
4297  log_newline
4298  log_info("MKINIT_boxaero",*) 'Setup initial state'
4299 
4300  !--- read namelist
4301  rewind(io_fid_conf)
4302  read(io_fid_conf,nml=param_mkinit_boxaero,iostat=ierr)
4303  if( ierr < 0 ) then !--- missing
4304  log_info("MKINIT_boxaero",*) 'Not found namelist. Default used.'
4305  elseif( ierr > 0 ) then !--- fatal error
4306  log_error("MKINIT_boxaero",*) 'Not appropriate names in namelist PARAM_MKINIT_BOXAERO. Check!'
4307  call prc_abort
4308  endif
4309  log_nml(param_mkinit_boxaero)
4310 
4311  call saturation_psat_all( init_temp, psat )
4312  qsat = epsvap * psat / ( init_pres - ( 1.0_rp-epsvap ) * psat )
4313 
4314  do j = 1, ja
4315  do i = 1, ia
4316  do k = 1, ka
4317  dens(k,i,j) = init_dens
4318  momx(k,i,j) = 0.0_rp
4319  momy(k,i,j) = 0.0_rp
4320  momz(k,i,j) = 0.0_rp
4321  pott(k,i,j) = init_temp * ( p00/init_pres )**(rdry/cpdry)
4322  rhot(k,i,j) = init_dens * pott(k,i,j)
4323 
4324  qv(k,i,j) = ( init_ssliq + 1.0_rp ) * qsat
4325 
4326  qdry = 1.0 - qv(k,i,j)
4327  rtot(k,i,j) = rdry * qdry + rvap * qv(i,i,j)
4328  cvtot(k,i,j) = cvdry * qdry + cvvap * qv(i,i,j)
4329  cptot(k,i,j) = cpdry * qdry + cpvap * qv(i,i,j)
4330  enddo
4331  enddo
4332  enddo
4333 
4334  call atmos_thermodyn_rhot2temp_pres( ka, 1, ka, ia, 1, ia, ja, 1, ja, &
4335  dens(:,:,:), rhot(:,:,:), & ! (in)
4336  rtot(:,:,:), cvtot(:,:,:), cptot(:,:,:), & ! (in)
4337  temp(:,:,:), pres(:,:,:) ) ! (out)
4338 
4339  return
4340  end subroutine mkinit_boxaero
4341 
4342  !-----------------------------------------------------------------------------
4344  subroutine mkinit_warmbubbleaero
4345  use scale_atmos_hydrometeor, only: &
4347  implicit none
4348 
4349  ! Surface state
4350  real(RP) :: SFC_THETA ! surface potential temperature [K]
4351  real(RP) :: SFC_PRES ! surface pressure [Pa]
4352  real(RP) :: SFC_RH = 80.0_rp ! surface relative humidity [%]
4353  ! Environment state
4354  real(RP) :: ENV_U = 0.0_rp ! velocity u of environment [m/s]
4355  real(RP) :: ENV_V = 0.0_rp ! velocity v of environment [m/s]
4356  real(RP) :: ENV_RH = 80.0_rp ! Relative Humidity of environment [%]
4357  real(RP) :: ENV_L1_ZTOP = 1.e3_rp ! top height of the layer1 (constant THETA) [m]
4358  real(RP) :: ENV_L2_ZTOP = 14.e3_rp ! top height of the layer2 (small THETA gradient) [m]
4359  real(RP) :: ENV_L2_TLAPS = 4.e-3_rp ! Lapse rate of THETA in the layer2 (small THETA gradient) [K/m]
4360  real(RP) :: ENV_L3_TLAPS = 3.e-2_rp ! Lapse rate of THETA in the layer3 (large THETA gradient) [K/m]
4361  ! Bubble
4362  real(RP) :: BBL_THETA = 1.0_rp ! extremum of temperature in bubble [K]
4363 
4364  namelist / param_mkinit_warmbubble / &
4365  sfc_theta, &
4366  sfc_pres, &
4367  env_u, &
4368  env_v, &
4369  env_rh, &
4370  env_l1_ztop, &
4371  env_l2_ztop, &
4372  env_l2_tlaps, &
4373  env_l3_tlaps, &
4374  bbl_theta
4375 
4376  integer :: ierr
4377  integer :: k, i, j
4378  !---------------------------------------------------------------------------
4379 
4380  log_newline
4381  log_info("MKINIT_warmbubbleaero",*) 'Setup initial state'
4382 
4383  if ( atmos_hydrometeor_dry ) then
4384  log_error("MKINIT_warmbubbleaero",*) 'QV is not registerd'
4385  call prc_abort
4386  end if
4387 
4388 
4389  sfc_theta = thetastd
4390  sfc_pres = pstd
4391 
4392  !--- read namelist
4393  rewind(io_fid_conf)
4394  read(io_fid_conf,nml=param_mkinit_warmbubble,iostat=ierr)
4395 
4396  if( ierr < 0 ) then !--- missing
4397  log_info("MKINIT_warmbubbleaero",*) 'Not found namelist. Default used.'
4398  elseif( ierr > 0 ) then !--- fatal error
4399  log_error("MKINIT_warmbubbleaero",*) 'Not appropriate names in namelist PARAM_MKINIT_WARMBUBBLE. Check!'
4400  call prc_abort
4401  endif
4402  log_nml(param_mkinit_warmbubble)
4403 
4404  ! calc in dry condition
4405  pres_sfc(1,1) = sfc_pres
4406  pott_sfc(1,1) = sfc_theta
4407 
4408  do k = ks, ke
4409  if ( cz(k) <= env_l1_ztop ) then ! Layer 1
4410  pott(k,1,1) = sfc_theta
4411  elseif( cz(k) < env_l2_ztop ) then ! Layer 2
4412  pott(k,1,1) = pott(k-1,1,1) + env_l2_tlaps * ( cz(k)-cz(k-1) )
4413  else ! Layer 3
4414  pott(k,1,1) = pott(k-1,1,1) + env_l3_tlaps * ( cz(k)-cz(k-1) )
4415  endif
4416  enddo
4417 
4418  ! make density & pressure profile in dry condition
4419  call hydrostatic_buildrho( ka, ks, ke, &
4420  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
4421  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
4422  cz(:), fz(:), & ! [IN]
4423  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
4424 
4425  ! calc QV from RH
4426  call saturation_psat_all( temp_sfc(1,1), psat_sfc(1,1) ) ! [IN], [OUT]
4427  qsat_sfc(1,1) = epsvap * psat_sfc(1,1) / ( pres_sfc(1,1) - ( 1.0_rp-epsvap ) * psat_sfc(1,1) )
4428 
4429  qdry(:,1,1) = 1.0_rp - qv(:,1,1) - qc(:,1,1)
4430  call saturation_pres2qsat_all( ka, ks, ke, &
4431  temp(:,1,1), pres(:,1,1), qdry(:,1,1), & ! [IN]
4432  qsat(:,1,1) ) ! [OUT]
4433  qv_sfc(1,1) = sfc_rh * 1.e-2_rp * qsat_sfc(1,1)
4434  do k = ks, ke
4435  if ( cz(k) <= env_l1_ztop ) then ! Layer 1
4436  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
4437  elseif( cz(k) <= env_l2_ztop ) then ! Layer 2
4438  qv(k,1,1) = env_rh * 1.e-2_rp * qsat(k,1,1)
4439  else ! Layer 3
4440  qv(k,1,1) = 0.0_rp
4441  endif
4442  enddo
4443 
4444  ! make density & pressure profile in moist condition
4445  call hydrostatic_buildrho( ka, ks, ke, &
4446  pott(:,1,1), qv(:,1,1), qc(:,1,1), & ! [IN]
4447  pres_sfc(1,1), pott_sfc(1,1), qv_sfc(1,1), qc_sfc(1,1), & ! [IN]
4448  cz(:), fz(:), & ! [IN]
4449  dens(:,1,1), temp(:,1,1), pres(:,1,1), temp_sfc(1,1) ) ! [OUT]
4450 
4451  do j = jsb, jeb
4452  do i = isb, ieb
4453  do k = ks, ke
4454  dens(k,i,j) = dens(k,1,1)
4455  momz(k,i,j) = 0.0_rp
4456  momx(k,i,j) = env_u * dens(k,i,j)
4457  momy(k,i,j) = env_v * dens(k,i,j)
4458 
4459  ! make warm bubble
4460  rhot(k,i,j) = dens(k,1,1) * ( pott(k,1,1) + bbl_theta * bubble(k,i,j) )
4461 
4462  qv(k,i,j) = qv(k,1,1)
4463  enddo
4464  enddo
4465  enddo
4466 
4467  call flux_setup
4468 
4469  return
4470  end subroutine mkinit_warmbubbleaero
4471 
4472  !-----------------------------------------------------------------------------
4474  subroutine mkinit_real
4475  use mod_realinput, only: &
4476  realinput_atmos, &
4478  implicit none
4479 
4480  call prof_rapstart('__Real_Atmos',2)
4481 
4482  call realinput_atmos
4483 
4484  call prof_rapend ('__Real_Atmos',2)
4485  call prof_rapstart('__Real_Surface',2)
4486 
4487  call realinput_surface
4488 
4489  call prof_rapend ('__Real_Surface',2)
4490 
4491  call flux_setup
4492 
4493  return
4494  end subroutine mkinit_real
4495 
4496 end module mod_mkinit
module ATMOS admin
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cxg
center coordinate [m]: x, global
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:57
real(rp), dimension(:,:,:,:), allocatable, public atmos_phy_rd_sflx_down
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
module atmosphere / saturation
real(rp), dimension(:,:), allocatable, public urban_qc
module coupler / surface-atmospehre
integer, parameter, public i_oceancouple
Definition: mod_mkinit.F90:112
real(rp), dimension(:,:,:), allocatable, target, public rhot
integer, public imax
of computational cells: x, local
integer, public ia
of whole cells: x, local, with HALO
module Atmosphere / Physics Cloud Microphysics
subroutine, public landuse_fillhalo(FILL_BND)
HALO Communication.
integer, parameter, public i_r_vis
integer, parameter, public i_boxaero
Definition: mod_mkinit.F90:125
integer, public iag
of computational grids
subroutine land_setup
Land setup.
Definition: mod_mkinit.F90:919
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:66
module ATMOSPHERIC Variables
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:44
real(rp), dimension(:,:,:), allocatable, target, public momx
real(rp), dimension(:,:), allocatable, public ocean_sfc_temp
ocean surface skin temperature [K]
subroutine ocean_setup
Ocean setup.
Definition: mod_mkinit.F90:970
subroutine, public landuse_calc_fact
real(rp), dimension(:,:), allocatable, public ocean_ice_mass
sea ice mass [kg]
integer, parameter, public i_dycoms2_rf02_dns
Definition: mod_mkinit.F90:120
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fxg
face coordinate [m]: x, global
integer, parameter, public i_turbulence
Definition: mod_mkinit.F90:98
integer, parameter, public i_urbancouple
Definition: mod_mkinit.F90:113
integer, parameter, public i_gravitywave
Definition: mod_mkinit.F90:96
real(rp), dimension(:,:,:), allocatable, public atmos_grid_cartesc_real_fz
geopotential height [m] (wxy)
integer, public ja
of whole cells: y, local, with HALO
module process / cartesC
character(len=h_short), public atmos_phy_ae_type
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
real(rp) function faero(f0, r0, x, alpha, rhoa)
Definition: mod_mkinit.F90:826
real(rp), dimension(:,:), allocatable, public urban_tb
integer, parameter, public i_grayzone
Definition: mod_mkinit.F90:124
real(rp), dimension(:,:,:), allocatable, target, public dens
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
Definition: scale_const.F90:65
real(rp), dimension(:,:), allocatable, public landuse_frac_urban
urban fraction
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fz
face coordinate [m]: z, local
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:111
real(rp), dimension(:,:), allocatable, public atmos_phy_mp_sflx_rain
real(rp), dimension(:,:), allocatable, public urban_uc
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_sw_up
real(rp), dimension(:,:,:,:), allocatable, public urban_sfc_albedo
module COMMUNICATION
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_lw_dn
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cy
center coordinate [m]: y, local
module Atmosphere / Physics Radiation
subroutine, public atmos_phy_ae_kajino13_mkinit(KA, KS, KE, IA, IS, IE, JA, JS, JE, QA_AE, DENS, TEMP, PRES, QDRY, QV, m0_init, dg_init, sg_init, d_min_inp, d_max_inp, k_min_inp, k_max_inp, n_kap_inp, QTRC, CCN)
real(rp), dimension(:,:), allocatable, public land_sfc_temp
land surface skin temperature [K]
integer, public is
start point of inner domain: x, local
real(rp), public const_ohm
angular velocity of the planet [1/s]
Definition: scale_const.F90:45
integer, public ie
end point of inner domain: x, local
module TRACER
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_sw_up
subroutine tke_setup
TKE setup.
real(rp), dimension(:,:), allocatable, public urban_tr
module REAL input
module atmosphere / vertical profile
integer, parameter, public i_seabreeze
Definition: mod_mkinit.F90:117
real(rp), dimension(:,:,:), allocatable, public urban_tgl
module atmosphere / hydrometeor
real(rp), dimension(:,:,:,:), allocatable, public land_sfc_albedo
land surface albedo (direct/diffuse,IR/near-IR/VIS) (0-1)
subroutine rect_setup
Bubble.
Definition: mod_mkinit.F90:612
subroutine, public atmos_phy_mp_driver_qhyd2qtrc(KA, KS, KE, IA, IS, IE, JA, JS, JE, QV, QHYD, QTRC, QNUM)
real(rp), dimension(:,:,:), allocatable, public ocean_salt
ocean salinity [PSU]
module LANDUSE
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0h
ocean surface roughness length for heat [m]
module atmosphere / grid / cartesC index
integer, public ke
end point of inner domain: z, local
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_sflx_sw_dn
integer, parameter, public i_warmbubble
Definition: mod_mkinit.F90:101
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
real(rp), dimension(:,:,:), allocatable, public ocean_temp
ocean temperature [K]
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_io.F90:313
module PROCESS
Definition: scale_prc.F90:11
integer, public je
end point of inner domain: y, local
subroutine read_sounding(DENS, VELX, VELY, POTT, QV)
Read sounding data from file.
integer, parameter, public i_mountainwave
Definition: mod_mkinit.F90:99
subroutine, public landuse_write
Write landuse data.
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cx
center coordinate [m]: x, local
real(rp), dimension(:,:,:), allocatable, public land_temp
temperature of each soil layer [K]
module atmosphere / hydrostatic barance
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
module atmosphere / physics / PBL
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_lw_dn
module atmosphere / grid / cartesC
real(rp), dimension(:,:), allocatable, public urban_roff
integer, public ks
start point of inner domain: z, local
real(rp), dimension(:,:,:), allocatable, public atmos_phy_ae_ccn
integer, parameter, public i_real
Definition: mod_mkinit.F90:122
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:69
module Atmosphere / Physics Turbulence
integer, parameter, public i_heatisland
Definition: mod_mkinit.F90:118
integer, parameter, public i_triplecouple
Definition: mod_mkinit.F90:114
module LAND Variables
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
integer, public jag
of computational grids
module CONSTANT
Definition: scale_const.F90:11
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fy
face coordinate [m]: y, local
integer, parameter, public i_hc
liquid water cloud
integer, public js
start point of inner domain: y, local
integer, parameter, public i_planestate
Definition: mod_mkinit.F90:91
integer, parameter, public i_r_direct
character(len=h_short), public atmos_phy_mp_type
integer, parameter, public i_khwave
Definition: mod_mkinit.F90:97
integer, parameter, public i_squallline
Definition: mod_mkinit.F90:103
real(rp), dimension(:,:), allocatable, public urban_tc
real(rp), dimension(:,:,:), allocatable, target, public momy
real(rp), dimension(:,:), allocatable, public urban_rainr
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:157
real(dp), parameter, public const_undef8
undefined value (REAL8)
Definition: scale_const.F90:40
real(rp), dimension(:,:,:), allocatable, public land_water
moisture of each soil layer [m3/m3]
integer, public mkinit_type
Definition: mod_mkinit.F90:88
subroutine, public mkinit(output)
Driver.
Definition: mod_mkinit.F90:343
real(rp), dimension(:,:), allocatable, public atmos_phy_mp_sflx_snow
integer, parameter, public i_rico
Definition: mod_mkinit.F90:107
real(rp), dimension(:,:,:), allocatable, public atmos_grid_cartesc_real_cz
geopotential height [m] (zxy)
module profiler
Definition: scale_prof.F90:11
integer, parameter, public i_ignore
Definition: mod_mkinit.F90:89
real(rp), public const_eps
small number
Definition: scale_const.F90:33
integer, parameter, public i_r_nir
module atmosphere / thermodyn
integer, parameter, public i_lambwave
Definition: mod_mkinit.F90:95
integer, parameter, public i_supercell
Definition: mod_mkinit.F90:102
real(rp), dimension(:,:,:), allocatable, public ocean_vvel
ocean meridional velocity [m/s]
module Atmosphere GRID CartesC Real(real space)
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_sw_dn
real(rp), dimension(:,:), allocatable, public atmos_grid_cartesc_real_area
horizontal area ( xy, normal z) [m2]
integer, parameter, public i_bomex
Definition: mod_mkinit.F90:130
module PRECISION
module atmosphere / physics / aerosol / Kajino13
integer, parameter, public i_dycoms2_rf01
Definition: mod_mkinit.F90:105
integer, public ka
of whole cells: z, local, with HALO
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fx
face coordinate [m]: x, local
integer, parameter, public i_interporation
Definition: mod_mkinit.F90:109
subroutine, public realinput_atmos
real(rp), public const_pi
pi
Definition: scale_const.F90:31
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fyg
face coordinate [m]: y, global
real(rp), dimension(:,:), allocatable, public urban_tg
real(rp), public atmos_grid_cartesc_domain_center_x
center position of global domain [m]: x
real(rp), dimension(:,:), allocatable, public urban_sfc_temp
subroutine, public random_uniform(var)
Get uniform random number.
real(rp), dimension(:,:,:), allocatable, public urban_trl
integer, parameter, public i_r_ir
integer, parameter, public i_warmbubbleaero
Definition: mod_mkinit.F90:126
integer, parameter, public i_wk1982
Definition: mod_mkinit.F90:104
real(rp), dimension(:,:), allocatable, public ocean_ice_temp
sea ice temperature [K]
module STDIO
Definition: scale_io.F90:10
real(rp), dimension(:,:), allocatable, public ocean_ocn_z0m
surface roughness length for momentum, open ocean [m]
module INITIAL
Definition: mod_mkinit.F90:12
integer, parameter, public i_bubblecouple
Definition: mod_mkinit.F90:115
integer, parameter, public i_r_diffuse
module RANDOM
integer, parameter, public i_tracerbubble
Definition: mod_mkinit.F90:92
integer, parameter, public n_hyd
real(rp), dimension(:,:,:,:), allocatable, public ocean_sfc_albedo
ocean surface albedo (direct/diffuse,IR/near-IR/VIS) (0-1)
integer, parameter, public i_coldbubble
Definition: mod_mkinit.F90:93
subroutine, public realinput_surface
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:64
integer, parameter, public i_cavityflow
Definition: mod_mkinit.F90:128
module Spectran Bin Microphysics
module atmosphere / physics / cloud microphysics
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:210
real(rp), dimension(:,:,:), allocatable, public urban_tbl
module ATMOSPHERE / Physics Aerosol Microphysics
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cz
center coordinate [m]: z, local
real(rp), public atmos_grid_cartesc_domain_center_y
center position of global domain [m]: y
integer, parameter, public i_barocwave
Definition: mod_mkinit.F90:129
real(rp), dimension(:,:), allocatable, public landuse_frac_land
land fraction
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0m
ocean surface roughness length for momentum [m]
subroutine, public mkinit_setup
Setup.
Definition: mod_mkinit.F90:210
integer, parameter, public i_dycoms2_rf02
Definition: mod_mkinit.F90:106
real(rp), dimension(:,:,:), allocatable, public ocean_uvel
ocean zonal velocity [m/s]
module OCEAN Variables
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:87
real(rp), dimension(:,:), allocatable, public urban_rainb
real(rp), dimension(:,:), allocatable, public ocean_sfc_z0e
ocean surface roughness length for vapor [m]
real(rp), dimension(:,:), allocatable, public atmos_phy_rd_toaflx_lw_up
subroutine urban_setup
Urban setup.
subroutine flux_setup
flux setup
Definition: mod_mkinit.F90:845
real(rp), dimension(:,:,:,:), allocatable, target, public qtrc
integer, public prc_num_x
x length of 2D processor topology