SCALE-RM
scale_atmos_sub_refstate.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
15 !-------------------------------------------------------------------------------
17  !-----------------------------------------------------------------------------
18  !
19  !++ used modules
20  !
21  use scale_precision
22  use scale_stdio
23  use scale_prof
25  use scale_tracer
26  !-----------------------------------------------------------------------------
27  implicit none
28  private
29  !-----------------------------------------------------------------------------
30  !
31  !++ Public procedure
32  !
33  public :: atmos_refstate_setup
34  public :: atmos_refstate_resume
35  public :: atmos_refstate_read
36  public :: atmos_refstate_write
37  public :: atmos_refstate_update
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Public parameters & variables
42  !
43  logical, public :: atmos_refstate_update_flag = .false.
44 
45  real(RP), public, allocatable :: atmos_refstate_pres(:,:,:)
46  real(RP), public, allocatable :: atmos_refstate_temp(:,:,:)
47  real(RP), public, allocatable :: atmos_refstate_dens(:,:,:)
48  real(RP), public, allocatable :: atmos_refstate_pott(:,:,:)
49  real(RP), public, allocatable :: atmos_refstate_qv (:,:,:)
50 
51  !-----------------------------------------------------------------------------
52  !
53  !++ Private procedure
54  !
55  private :: atmos_refstate_generate_isa
56  private :: atmos_refstate_generate_uniform
57  private :: atmos_refstate_generate_zero
58  private :: atmos_refstate_generate_frominit
59 
60  !-----------------------------------------------------------------------------
61  !
62  !++ Private parameters & variables
63  !
64  character(len=H_LONG), private :: atmos_refstate_in_basename = ''
65  logical, private :: atmos_refstate_in_check_coordinates = .true.
66  character(len=H_LONG), private :: atmos_refstate_out_basename = ''
67  character(len=H_MID), private :: atmos_refstate_out_title = 'SCALE-RM RefState'
68  character(len=H_SHORT), private :: atmos_refstate_out_dtype = 'DEFAULT'
69 
70  character(len=H_SHORT), private :: atmos_refstate_type = 'UNIFORM'
71  real(RP), private :: atmos_refstate_temp_sfc = 300.0_rp
72  real(RP), private :: atmos_refstate_rh = 0.0_rp
73  real(RP), private :: atmos_refstate_pott_uniform = 300.0_rp
74  real(DP), private :: atmos_refstate_update_dt = 0.0_dp
75 
76  real(DP), private :: last_updated
77 
78  real(RP), private, allocatable :: atmos_refstate1d_pres(:)
79  real(RP), private, allocatable :: atmos_refstate1d_temp(:)
80  real(RP), private, allocatable :: atmos_refstate1d_dens(:)
81  real(RP), private, allocatable :: atmos_refstate1d_pott(:)
82  real(RP), private, allocatable :: atmos_refstate1d_qv (:)
83 
84  !-----------------------------------------------------------------------------
85 contains
86  !-----------------------------------------------------------------------------
88  subroutine atmos_refstate_setup
89  use scale_process, only: &
91  implicit none
92 
93  namelist / param_atmos_refstate / &
94  atmos_refstate_in_basename, &
95  atmos_refstate_out_basename, &
96  atmos_refstate_out_title, &
97  atmos_refstate_out_dtype, &
98  atmos_refstate_type, &
99  atmos_refstate_temp_sfc, &
100  atmos_refstate_rh, &
101  atmos_refstate_pott_uniform, &
103  atmos_refstate_update_dt
104 
105  integer :: ierr
106  !---------------------------------------------------------------------------
107 
108  if( io_l ) write(io_fid_log,*)
109  if( io_l ) write(io_fid_log,*) '++++++ Module[REFSTATE] / Categ[ATMOS SHARE] / Origin[SCALElib]'
110 
111  allocate( atmos_refstate_pres(ka,ia,ja) )
112  allocate( atmos_refstate_temp(ka,ia,ja) )
113  allocate( atmos_refstate_dens(ka,ia,ja) )
114  allocate( atmos_refstate_pott(ka,ia,ja) )
115  allocate( atmos_refstate_qv(ka,ia,ja) )
116 
117  allocate( atmos_refstate1d_pres(ka) )
118  allocate( atmos_refstate1d_temp(ka) )
119  allocate( atmos_refstate1d_dens(ka) )
120  allocate( atmos_refstate1d_pott(ka) )
121  allocate( atmos_refstate1d_qv(ka) )
122 
123  !--- read namelist
124  rewind(io_fid_conf)
125  read(io_fid_conf,nml=param_atmos_refstate,iostat=ierr)
126  if( ierr < 0 ) then !--- missing
127  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
128  elseif( ierr > 0 ) then !--- fatal error
129  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_REFSTATE. Check!'
130  call prc_mpistop
131  endif
132  if( io_nml ) write(io_fid_nml,nml=param_atmos_refstate)
133 
134  if( io_l ) write(io_fid_log,*)
135  if ( atmos_refstate_in_basename /= '' ) then
136  if( io_l ) write(io_fid_log,*) '*** Input file of reference state : ', trim(atmos_refstate_in_basename)
137  else
138  if( io_l ) write(io_fid_log,*) '*** Input file of reference state : Nothing, generate internally'
139  endif
140 
141  ! input or generate reference profile
142  if ( atmos_refstate_in_basename /= '' ) then
144  else
145  if ( atmos_refstate_type == 'ISA' ) then
146 
147  if( io_l ) write(io_fid_log,*) '*** Reference type : ISA'
148  if( io_l ) write(io_fid_log,*) '*** Surface temperature [K] : ', atmos_refstate_temp_sfc
149  if( io_l ) write(io_fid_log,*) '*** Surface & environment RH [%] : ', atmos_refstate_rh
150  call atmos_refstate_generate_isa
152 
153  elseif ( atmos_refstate_type == 'UNIFORM' ) then
154 
155  if( io_l ) write(io_fid_log,*) '*** Reference type : UNIFORM POTT'
156  if( io_l ) write(io_fid_log,*) '*** Potential temperature : ', atmos_refstate_pott_uniform
157  call atmos_refstate_generate_uniform
159 
160  elseif ( atmos_refstate_type == 'ZERO' ) then
161 
162  if( io_l ) write(io_fid_log,*) '*** Reference type : ZERO'
163  call atmos_refstate_generate_zero
165 
166  elseif ( atmos_refstate_type == 'INIT' ) then
167 
168  if( io_l ) write(io_fid_log,*) '*** Reference type : Generate from initial data'
169  if( io_l ) write(io_fid_log,*) '*** Update state? : ', atmos_refstate_update_flag
170  if( io_l ) write(io_fid_log,*) '*** Update interval [sec] : ', atmos_refstate_update_dt
171 
172  else
173  write(*,*) 'xxx ATMOS_REFSTATE_TYPE must be "ISA" or "UNIFORM". Check! : ', trim(atmos_refstate_type)
174  call prc_mpistop
175  endif
176 
177  endif
178 
179  return
180  end subroutine atmos_refstate_setup
181 
182  !-----------------------------------------------------------------------------
184  subroutine atmos_refstate_resume( &
185  DENS, RHOT, QTRC )
186  use scale_process, only: &
188  use scale_grid, only: &
189  cz => grid_cz
190  implicit none
191 
192  real(RP), intent(in) :: dens(ka,ia,ja)
193  real(RP), intent(in) :: rhot(ka,ia,ja)
194  real(RP), intent(in) :: qtrc(ka,ia,ja,qa)
195 
196  integer :: k
197 
198  ! input or generate reference profile
199  if ( atmos_refstate_in_basename == '' ) then
200 
201  if ( atmos_refstate_type == 'INIT' ) then
202 
203  if( io_l ) write(io_fid_log,*) '*** Reference type : make from initial data'
204  if( io_l ) write(io_fid_log,*) '*** Update state? : ', atmos_refstate_update_flag
205  if( io_l ) write(io_fid_log,*) '*** Update interval [sec] : ', atmos_refstate_update_dt
206  call atmos_refstate_generate_frominit( dens, rhot, qtrc ) ! (in)
207 
208  endif
209 
210  if( io_l ) write(io_fid_log,*)
211  if( io_l ) write(io_fid_log,*) '###### Generated Reference State of Atmosphere ######'
212  if( io_l ) write(io_fid_log,*) ' z*-coord.: pressure: temperature: density: pot.temp.: water vapor'
213  do k = ks, ke
214  if( io_l ) write(io_fid_log,'(6F13.5)') cz(k), &
215  atmos_refstate1d_pres(k), &
216  atmos_refstate1d_temp(k), &
217  atmos_refstate1d_dens(k), &
218  atmos_refstate1d_pott(k), &
219  atmos_refstate1d_qv(k)
220  enddo
221  if( io_l ) write(io_fid_log,*) '####################################################'
222  endif
223 
224  if ( atmos_refstate_out_basename /= '' ) then
225  if( io_l ) write(io_fid_log,*) '*** Reference state output? : ', trim(atmos_refstate_out_basename)
226  else
227  if( io_l ) write(io_fid_log,*) '*** Reference state output? : NO'
228  endif
229 
230  ! output reference profile
232 
233  return
234  end subroutine atmos_refstate_resume
235 
236  !-----------------------------------------------------------------------------
238  subroutine atmos_refstate_read
239  use scale_fileio, only: &
240  fileio_open, &
241  fileio_check_coordinates, &
242  fileio_read, &
244  use scale_process, only: &
246  implicit none
247 
248  integer :: fid
249  !---------------------------------------------------------------------------
250 
251  if( io_l ) write(io_fid_log,*)
252  if( io_l ) write(io_fid_log,*) '*** Input reference state profile ***'
253 
254  if ( atmos_refstate_in_basename /= '' ) then
255 
256  call fileio_open( fid, atmos_refstate_in_basename )
257 
258  if ( atmos_refstate_in_check_coordinates ) then
259  call fileio_check_coordinates( fid, atmos=.true. )
260  end if
261 
262  ! 1D
263  call fileio_read( atmos_refstate1d_pres(:), fid, 'PRES_ref', 'Z', step=1 )
264  call fileio_read( atmos_refstate1d_temp(:), fid, 'TEMP_ref', 'Z', step=1 )
265  call fileio_read( atmos_refstate1d_dens(:), fid, 'DENS_ref', 'Z', step=1 )
266  call fileio_read( atmos_refstate1d_pott(:), fid, 'POTT_ref', 'Z', step=1 )
267  call fileio_read( atmos_refstate1d_qv(:), fid, 'QV_ref', 'Z', step=1 )
268 
269  ! 3D
270  call fileio_read( atmos_refstate_pres(:,:,:), fid, 'PRES_ref3D', 'ZXY', step=1 )
271  call fileio_read( atmos_refstate_temp(:,:,:), fid, 'TEMP_ref3D', 'ZXY', step=1 )
272  call fileio_read( atmos_refstate_dens(:,:,:), fid, 'DENS_ref3D', 'ZXY', step=1 )
273  call fileio_read( atmos_refstate_pott(:,:,:), fid, 'POTT_ref3D', 'ZXY', step=1 )
274  call fileio_read( atmos_refstate_qv(:,:,:), fid, 'QV_ref3D', 'ZXY', step=1 )
275 
276  else
277  write(*,*) 'xxx [ATMOS_REFSTATE_read] refstate file is not specified.'
278  call prc_mpistop
279  endif
280 
282 
283 
284  return
285  end subroutine atmos_refstate_read
286 
287  !-----------------------------------------------------------------------------
289  subroutine atmos_refstate_write
290  use scale_fileio, only: &
291  fileio_write
292  implicit none
293  !---------------------------------------------------------------------------
294 
295  if ( atmos_refstate_out_basename /= '' ) then
296 
297  if( io_l ) write(io_fid_log,*)
298  if( io_l ) write(io_fid_log,*) '*** Output reference state profile ***'
299 
300  ! 1D
301  call fileio_write( atmos_refstate1d_pres(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
302  'PRES_ref', 'Reference profile of pres.', 'Pa', 'Z', atmos_refstate_out_dtype ) ! [IN]
303  call fileio_write( atmos_refstate1d_temp(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
304  'TEMP_ref', 'Reference profile of temp.', 'K', 'Z', atmos_refstate_out_dtype ) ! [IN]
305  call fileio_write( atmos_refstate1d_dens(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
306  'DENS_ref', 'Reference profile of rho', 'kg/m3', 'Z', atmos_refstate_out_dtype ) ! [IN]
307  call fileio_write( atmos_refstate1d_pott(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
308  'POTT_ref', 'Reference profile of theta', 'K', 'Z', atmos_refstate_out_dtype ) ! [IN]
309  call fileio_write( atmos_refstate1d_qv(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
310  'QV_ref', 'Reference profile of qv', 'kg/kg', 'Z', atmos_refstate_out_dtype ) ! [IN]
311 
312  ! 3D
313  call fileio_write( atmos_refstate_pres(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
314  'PRES_ref3D', 'Reference profile of pres.', 'Pa', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
315  call fileio_write( atmos_refstate_temp(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
316  'TEMP_ref3D', 'Reference profile of temp.', 'K', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
317  call fileio_write( atmos_refstate_dens(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
318  'DENS_ref3D', 'Reference profile of rho', 'kg/m3', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
319  call fileio_write( atmos_refstate_pott(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
320  'POTT_ref3D', 'Reference profile of theta', 'K', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
321  call fileio_write( atmos_refstate_qv(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
322  'QV_ref3D', 'Reference profile of qv', 'kg/kg', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
323 
324  endif
325 
326  return
327  end subroutine atmos_refstate_write
328 
329  !-----------------------------------------------------------------------------
331  subroutine atmos_refstate_generate_isa
332  use scale_const, only: &
333  epsvap => const_epsvap, &
334  pstd => const_pstd
335  use scale_grid_real, only: &
336  real_cz
337  use scale_comm, only: &
339  use scale_atmos_profile, only: &
340  profile_isa => atmos_profile_isa
341  use scale_atmos_hydrostatic, only: &
342  hydrostatic_buildrho => atmos_hydrostatic_buildrho
343  use scale_atmos_saturation, only: &
344  saturation_psat_all => atmos_saturation_psat_all, &
345  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
346  implicit none
347 
348  real(RP) :: z(ka)
349 
350  real(RP) :: temp(ka)
351  real(RP) :: pres(ka)
352  real(RP) :: dens(ka)
353  real(RP) :: pott(ka)
354  real(RP) :: qv (ka)
355  real(RP) :: qc (ka)
356 
357  real(RP) :: temp_sfc
358  real(RP) :: pres_sfc
359  real(RP) :: pott_sfc
360  real(RP) :: qv_sfc
361  real(RP) :: qc_sfc
362 
363  real(RP) :: qsat(ka)
364  real(RP) :: psat_sfc
365 
366  integer :: k
367  !---------------------------------------------------------------------------
368 
369  pott_sfc = atmos_refstate_temp_sfc
370  pres_sfc = pstd
371 
372  call comm_horizontal_mean( z(:), real_cz(:,:,:) )
373 
374  call profile_isa( ka, ks, ke, & ! [IN]
375  pott_sfc, & ! [IN]
376  pres_sfc, & ! [IN]
377  z(:), & ! [IN]
378  pott(:) ) ! [OUT]
379 
380  qv(:) = 0.0_rp
381  qc(:) = 0.0_rp
382  qv_sfc = 0.0_rp
383  qc_sfc = 0.0_rp
384 
385  ! make density & pressure profile in dry condition
386  call hydrostatic_buildrho( dens(:), & ! [OUT]
387  temp(:), & ! [OUT]
388  pres(:), & ! [OUT]
389  pott(:), & ! [IN]
390  qv(:), & ! [IN]
391  qc(:), & ! [IN]
392  temp_sfc, & ! [OUT]
393  pres_sfc, & ! [IN]
394  pott_sfc, & ! [IN]
395  qv_sfc, & ! [IN]
396  qc_sfc ) ! [IN]
397 
398  ! calc QV from RH
399  call saturation_psat_all( psat_sfc, temp_sfc )
400  call saturation_dens2qsat_all( qsat(:), temp(:), dens(:) )
401 
402  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
403  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp-epsvap) * psat_sfc )
404  do k = ks, ke
405  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
406  enddo
407 
408  ! make density & pressure profile in moist condition
409  call hydrostatic_buildrho( dens(:), & ! [OUT]
410  temp(:), & ! [OUT]
411  pres(:), & ! [OUT]
412  pott(:), & ! [IN]
413  qv(:), & ! [IN]
414  qc(:), & ! [IN]
415  temp_sfc, & ! [OUT]
416  pres_sfc, & ! [IN]
417  pott_sfc, & ! [IN]
418  qv_sfc, & ! [IN]
419  qc_sfc ) ! [IN]
420 
421  atmos_refstate1d_pres(:) = pres(:)
422  atmos_refstate1d_temp(:) = temp(:)
423  atmos_refstate1d_dens(:) = dens(:)
424  atmos_refstate1d_pott(:) = pott(:)
425  atmos_refstate1d_qv(:) = qv(:)
426 
428 
429  return
430  end subroutine atmos_refstate_generate_isa
431 
432  !-----------------------------------------------------------------------------
434  subroutine atmos_refstate_generate_uniform
435  use scale_const, only: &
436  epsvap => const_epsvap, &
437  pstd => const_pstd
438  use scale_atmos_hydrostatic, only: &
439  hydrostatic_buildrho => atmos_hydrostatic_buildrho
440  use scale_atmos_saturation, only: &
441  saturation_psat_all => atmos_saturation_psat_all, &
442  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
443  implicit none
444 
445  real(RP) :: temp(ka)
446  real(RP) :: pres(ka)
447  real(RP) :: dens(ka)
448  real(RP) :: pott(ka)
449  real(RP) :: qv (ka)
450  real(RP) :: qc (ka)
451 
452  real(RP) :: temp_sfc
453  real(RP) :: pres_sfc
454  real(RP) :: pott_sfc
455  real(RP) :: qv_sfc
456  real(RP) :: qc_sfc
457 
458  real(RP) :: qsat(ka)
459  real(RP) :: psat_sfc
460 
461  integer :: k
462  !---------------------------------------------------------------------------
463 
464  pres_sfc = pstd
465  pott_sfc = atmos_refstate_temp_sfc
466  qv_sfc = 0.0_rp
467  qc_sfc = 0.0_rp
468 
469  do k = 1, ka
470  pott(k) = atmos_refstate_pott_uniform
471  qv(k) = 0.0_rp
472  qc(k) = 0.0_rp
473  enddo
474 
475  ! make density & pressure profile in dry condition
476  call hydrostatic_buildrho( dens(:), & ! [OUT]
477  temp(:), & ! [OUT]
478  pres(:), & ! [OUT]
479  pott(:), & ! [IN]
480  qv(:), & ! [IN]
481  qc(:), & ! [IN]
482  temp_sfc, & ! [OUT]
483  pres_sfc, & ! [IN]
484  pott_sfc, & ! [IN]
485  qv_sfc, & ! [IN]
486  qc_sfc ) ! [IN]
487 
488  ! calc QV from RH
489  call saturation_psat_all( psat_sfc, temp_sfc )
490  call saturation_dens2qsat_all( qsat(:), temp(:), pres(:) )
491 
492  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
493  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp - epsvap) * psat_sfc )
494  do k = ks, ke
495  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
496  enddo
497 
498  ! make density & pressure profile in moist condition
499  call hydrostatic_buildrho( dens(:), & ! [OUT]
500  temp(:), & ! [OUT]
501  pres(:), & ! [OUT]
502  pott(:), & ! [IN]
503  qv(:), & ! [IN]
504  qc(:), & ! [IN]
505  temp_sfc, & ! [OUT]
506  pres_sfc, & ! [IN]
507  pott_sfc, & ! [IN]
508  qv_sfc, & ! [IN]
509  qc_sfc ) ! [IN]
510 
511  atmos_refstate1d_pres(:) = pres(:)
512  atmos_refstate1d_temp(:) = temp(:)
513  atmos_refstate1d_dens(:) = dens(:)
514  atmos_refstate1d_pott(:) = pott(:)
515  atmos_refstate1d_qv(:) = qv(:)
516 
518 
519  return
520  end subroutine atmos_refstate_generate_uniform
521 
522  !-----------------------------------------------------------------------------
524  subroutine atmos_refstate_generate_zero
525  implicit none
526 
527  integer :: k, i, j
528  !---------------------------------------------------------------------------
529 
530  do k = 1, ka
531  do i = 1, ia
532  do j = 1, ja
533  atmos_refstate_dens(k,i,j) = 0.0_rp
534  atmos_refstate_temp(k,i,j) = 0.0_rp
535  atmos_refstate_pres(k,i,j) = 0.0_rp
536  atmos_refstate_pott(k,i,j) = 0.0_rp
537  atmos_refstate_qv(k,i,j) = 0.0_rp
538  enddo
539  enddo
540  enddo
541 
542  return
543  end subroutine atmos_refstate_generate_zero
544 
545  !-----------------------------------------------------------------------------
547  subroutine atmos_refstate_generate_frominit( &
548  DENS, RHOT, QTRC )
549  use scale_time, only: &
551  implicit none
552 
553  real(RP), intent(in) :: dens(ka,ia,ja)
554  real(RP), intent(in) :: rhot(ka,ia,ja)
555  real(RP), intent(in) :: qtrc(ka,ia,ja,qa)
556  !---------------------------------------------------------------------------
557 
558  last_updated = time_nowsec - atmos_refstate_update_dt
559 
560  call atmos_refstate_update( dens, rhot, qtrc ) ! (in)
561 
562  return
563  end subroutine atmos_refstate_generate_frominit
564 
565  !-----------------------------------------------------------------------------
567  subroutine atmos_refstate_update( &
568  DENS, RHOT, QTRC )
569  use scale_time, only: &
571  use scale_comm, only: &
573  use scale_interpolation, only: &
575  use scale_atmos_thermodyn, only: &
576  thermodyn_temp_pres => atmos_thermodyn_temp_pres
577  use scale_atmos_hydrometeor, only: &
578  i_qv
579  implicit none
580  real(RP), intent(in) :: dens(ka,ia,ja)
581  real(RP), intent(in) :: rhot(ka,ia,ja)
582  real(RP), intent(in) :: qtrc(ka,ia,ja,qa)
583 
584  real(RP) :: temp(ka,ia,ja)
585  real(RP) :: pres(ka,ia,ja)
586  real(RP) :: pott(ka,ia,ja)
587  real(RP) :: work(ka,ia,ja)
588 
589  integer :: k, i, j
590  !---------------------------------------------------------------------------
591 
592  if ( time_nowsec - last_updated >= atmos_refstate_update_dt ) then
593 
594  if( io_l ) write(io_fid_log,*) '*** [REFSTATE] update reference state'
595 
596  call thermodyn_temp_pres( temp, & ! [OUT]
597  pres, & ! [OUT]
598  dens, & ! [IN]
599  rhot, & ! [IN]
600  qtrc, & ! [IN]
601  tracer_cv, & ! [IN]
602  tracer_cp, & ! [IN]
603  tracer_mass ) ! [IN]
604 
605  do j = 1, ja
606  do i = 1, ia
607  do k = ks, ke
608  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
609  enddo
610  enddo
611  enddo
612 
613  call interp_vertical_xi2z( temp(:,:,:), & ! [IN]
614  work(:,:,:) ) ! [OUT]
615 
616  call comm_horizontal_mean( atmos_refstate1d_temp(:), work(:,:,:) )
617 
618  call interp_vertical_xi2z( pres(:,:,:), & ! [IN]
619  work(:,:,:) ) ! [OUT]
620 
621  call comm_horizontal_mean( atmos_refstate1d_pres(:), work(:,:,:) )
622 
623  call interp_vertical_xi2z( dens(:,:,:), & ! [IN]
624  work(:,:,:) ) ! [OUT]
625 
626  call comm_horizontal_mean( atmos_refstate1d_dens(:), work(:,:,:) )
627 
628  call interp_vertical_xi2z( pott(:,:,:), & ! [IN]
629  work(:,:,:) ) ! [OUT]
630 
631  call comm_horizontal_mean( atmos_refstate1d_pott(:), work(:,:,:) )
632 
633  if ( i_qv > 0 ) then
634  call interp_vertical_xi2z( qtrc(:,:,:,i_qv), & ! [IN]
635  work(:,:,:) ) ! [OUT]
636 
637  call comm_horizontal_mean( atmos_refstate1d_qv(:), work(:,:,:) )
638  else
639  atmos_refstate1d_qv(:) = 0.0_rp
640  endif
641 
642  do k = ke-1, ks, -1 ! fill undefined value
643  if( atmos_refstate1d_dens(k) <= 0.0_rp ) atmos_refstate1d_dens(k) = atmos_refstate1d_dens(k+1)
644  if( atmos_refstate1d_temp(k) <= 0.0_rp ) atmos_refstate1d_temp(k) = atmos_refstate1d_temp(k+1)
645  if( atmos_refstate1d_pres(k) <= 0.0_rp ) atmos_refstate1d_pres(k) = atmos_refstate1d_pres(k+1)
646  if( atmos_refstate1d_pott(k) <= 0.0_rp ) atmos_refstate1d_pott(k) = atmos_refstate1d_pott(k+1)
647  if( atmos_refstate1d_qv(k) <= 0.0_rp ) atmos_refstate1d_qv(k) = atmos_refstate1d_qv(k+1)
648  enddo
649  call smoothing( atmos_refstate1d_pott(:) )
650  if ( i_qv > 0 ) call smoothing( atmos_refstate1d_qv(:) )
651 
653 
654  last_updated = time_nowsec
655 
656  endif
657 
658  return
659  end subroutine atmos_refstate_update
660 
661  !-----------------------------------------------------------------------------
663  subroutine atmos_refstate_calc3d
664  use scale_const, only: &
665  undef => const_undef, &
666  rdry => const_rdry, &
667  cpdry => const_cpdry, &
668  p00 => const_pre00
669  use scale_comm, only: &
670  comm_vars8, &
671  comm_wait
672  use scale_grid, only: &
673  grid_cz, &
674  grid_fz
675  use scale_grid_real, only: &
676  real_phi, &
677  real_cz, &
678  real_fz
679  use scale_interpolation, only: &
681  use scale_atmos_hydrostatic, only: &
682  hydrostatic_buildrho_atmos_0d => atmos_hydrostatic_buildrho_atmos_0d, &
683  hydrostatic_buildrho_atmos_rev_2d => atmos_hydrostatic_buildrho_atmos_rev_2d, &
684  hydrostatic_buildrho_atmos_rev_3d => atmos_hydrostatic_buildrho_atmos_rev_3d
685  implicit none
686 
687 
688  real(RP) :: dens(KA,IA,JA)
689  real(RP) :: temp(KA,IA,JA)
690  real(RP) :: pres(KA,IA,JA)
691  real(RP) :: pott(KA,IA,JA)
692  real(RP) :: qv (KA,IA,JA)
693  real(RP) :: qc (KA,IA,JA)
694  real(RP) :: dz (KA,IA,JA)
695 
696  real(RP) :: dens_toa_1D
697  real(RP) :: temp_toa_1D
698  real(RP) :: pres_toa_1D
699  real(RP) :: qc_1D
700  real(RP) :: dz_1D
701 
702  real(RP) :: work(KA,IA,JA)
703  real(RP) :: RovCP
704  integer :: k, i, j
705  !---------------------------------------------------------------------------
706 
707  rovcp = rdry / cpdry
708 
709  !--- potential temperature
710  do j = jsb, jeb
711  do i = isb, ieb
712  work(:,i,j) = atmos_refstate1d_pott(:)
713  enddo
714  enddo
715 
716  call interp_vertical_z2xi( work(:,:,:), & ! [IN]
717  pott(:,:,:) ) ! [OUT]
718 
719  !--- water vapor
720  do j = jsb, jeb
721  do i = isb, ieb
722  work(:,i,j) = atmos_refstate1d_qv(:)
723  enddo
724  enddo
725 
726  call interp_vertical_z2xi( work(:,:,:), & ! [IN]
727  qv(:,:,:) ) ! [OUT]
728 
729 
730 
731  !--- build up density to TOA (1D)
732  qc_1d = 0.0_rp
733  dz_1d = grid_fz(ke) - grid_cz(ke)
734 
735  call hydrostatic_buildrho_atmos_0d( dens_toa_1d, & ! [OUT]
736  temp_toa_1d, & ! [OUT]
737  pres_toa_1d, & ! [OUT]
738  atmos_refstate1d_pott(ke), & ! [IN]
739  atmos_refstate1d_qv(ke), & ! [IN]
740  qc_1d, & ! [IN]
741  atmos_refstate1d_dens(ke), & ! [IN]
742  atmos_refstate1d_pott(ke), & ! [IN]
743  atmos_refstate1d_qv(ke), & ! [IN]
744  qc_1d, & ! [IN]
745  dz_1d, & ! [IN]
746  ke+1 ) ! [IN]
747 
748  ! build down density from TOA (3D)
749  do j = jsb, jeb
750  do i = isb, ieb
751  dz(ks,i,j) = real_cz(ks,i,j) - real_fz(ks-1,i,j) ! distance from surface to cell center
752  do k = ks+1, ke
753  dz(k,i,j) = real_cz(k,i,j) - real_cz(k-1,i,j) ! distance from cell center to cell center
754  enddo
755  dz(ke+1,i,j) = real_fz(ke,i,j) - real_cz(ke,i,j) ! distance from cell center to TOA
756  enddo
757  enddo
758 
759  do j = jsb, jeb
760  do i = isb, ieb
761  dens(ke+1,i,j) = dens_toa_1d
762  temp(ke+1,i,j) = temp_toa_1d
763  pres(ke+1,i,j) = pres_toa_1d
764  pott(ke+1,i,j) = pott(ke,i,j)
765  qv(ke+1,i,j) = qv(ke,i,j)
766  enddo
767  enddo
768 
769  do j = jsb, jeb
770  do i = isb, ieb
771  pott(ks-1,i,j) = pott(ks,i,j)
772  qv(ks-1,i,j) = qv(ks,i,j)
773  enddo
774  enddo
775 
776  qc(:,:,:) = 0.0_rp
777 
778  call hydrostatic_buildrho_atmos_rev_2d( dens(ke ,:,:), & ! [OUT]
779  temp(ke ,:,:), & ! [OUT]
780  pres(ke ,:,:), & ! [OUT]
781  pott(ke ,:,:), & ! [IN]
782  qv(ke ,:,:), & ! [IN]
783  qc(ke ,:,:), & ! [IN]
784  dens(ke+1,:,:), & ! [IN]
785  pott(ke+1,:,:), & ! [IN]
786  qv(ke+1,:,:), & ! [IN]
787  qc(ke+1,:,:), & ! [IN]
788  dz(ke+1,:,:), & ! [IN]
789  ke+1 ) ! [IN]
790 
791  call hydrostatic_buildrho_atmos_rev_3d( dens(:,:,:), & ! [INOUT]
792  temp(:,:,:), & ! [OUT]
793  pres(:,:,:), & ! [OUT]
794  pott(:,:,:), & ! [IN]
795  qv(:,:,:), & ! [IN]
796  qc(:,:,:), & ! [IN]
797  dz(:,:,:) ) ! [IN]
798 
799 ! call HYDROSTATIC_buildrho_atmos_rev_2D( dens(KS-1,:,:), & ! [OUT]
800 ! temp(KS-1,:,:), & ! [OUT]
801 ! pres(KS-1,:,:), & ! [OUT]
802 ! pott(KS-1,:,:), & ! [IN]
803 ! qv (KS-1,:,:), & ! [IN]
804 ! qc (KS-1,:,:), & ! [IN]
805 ! dens(KS ,:,:), & ! [IN]
806 ! pott(KS ,:,:), & ! [IN]
807 ! qv (KS ,:,:), & ! [IN]
808 ! qc (KS ,:,:), & ! [IN]
809 ! dz (KS ,:,:), & ! [IN]
810 ! KS ) ! [IN]
811  do j = jsb, jeb
812  do i = isb, ieb
813  do k = ks, ke
814  atmos_refstate_dens(k,i,j) = dens(k,i,j)
815  atmos_refstate_temp(k,i,j) = temp(k,i,j)
816  atmos_refstate_pres(k,i,j) = pres(k,i,j)
817  atmos_refstate_pott(k,i,j) = pott(k,i,j)
818  atmos_refstate_qv(k,i,j) = qv(k,i,j)
819  enddo
820  enddo
821  enddo
822 
823  ! boundary condition
824  do j = jsb, jeb
825  do i = isb, ieb
826 
827  atmos_refstate_temp(1:ks-1,i,j) = temp(ks,i,j)
828  atmos_refstate_temp(ke+1:ka,i,j) = temp_toa_1d
829 
830  atmos_refstate_qv(1:ks-1,i,j) = qv(ks,i,j)
831  atmos_refstate_qv(ke+1:ka,i,j) = qv(ke,i,j)
832 
833  atmos_refstate_pres(1:ks-2,i,j) = undef
834  atmos_refstate_pres(ks-1,i,j) = atmos_refstate_pres(ks+1,i,j) &
835  - atmos_refstate_dens(ks ,i,j) * ( real_phi(ks-1,i,j) - real_phi(ks+1,i,j) )
836  atmos_refstate_pres(ke+1,i,j) = atmos_refstate_pres(ke-1,i,j) &
837  - atmos_refstate_dens(ke ,i,j) * ( real_phi(ke+1,i,j) - real_phi(ke-1,i,j) )
838  atmos_refstate_pres(ke+2:ka,i,j) = undef
839 
840  atmos_refstate_dens(1:ks-2,i,j) = undef
841  atmos_refstate_dens(ks-1,i,j) = atmos_refstate_pres(ks-1,i,j) / ( atmos_refstate_temp(ks-1,i,j) * rdry )
842  atmos_refstate_dens(ke+1,i,j) = atmos_refstate_pres(ke+1,i,j) / ( atmos_refstate_temp(ke+1,i,j) * rdry )
843  atmos_refstate_dens(ke+2:ka,i,j) = undef
844 
845  atmos_refstate_pott(1:ks-2,i,j) = undef
846  atmos_refstate_pott(ks-1,i,j) = atmos_refstate_temp(ks-1,i,j) * ( p00 / atmos_refstate_pres(ks-1,i,j) )**rovcp
847  atmos_refstate_pott(ke+1,i,j) = atmos_refstate_temp(ke+1,i,j) * ( p00 / atmos_refstate_pres(ke+1,i,j) )**rovcp
848  atmos_refstate_pott(ke+2:ka,i,j) = undef
849  enddo
850  enddo
851 
852  call comm_vars8( atmos_refstate_dens(:,:,:), 1 )
853  call comm_vars8( atmos_refstate_temp(:,:,:), 2 )
854  call comm_vars8( atmos_refstate_pres(:,:,:), 3 )
855  call comm_vars8( atmos_refstate_pott(:,:,:), 4 )
856  call comm_vars8( atmos_refstate_qv(:,:,:), 5 )
857  call comm_wait ( atmos_refstate_dens(:,:,:), 1, .false. )
858  call comm_wait ( atmos_refstate_temp(:,:,:), 2, .false. )
859  call comm_wait ( atmos_refstate_pres(:,:,:), 3, .false. )
860  call comm_wait ( atmos_refstate_pott(:,:,:), 4, .false. )
861  call comm_wait ( atmos_refstate_qv(:,:,:), 5, .false. )
862 
863  return
864  end subroutine atmos_refstate_calc3d
865 
866  !-----------------------------------------------------------------------------
867  subroutine smoothing( &
868  phi )
869  use scale_process, only: &
871  use scale_const, only: &
872  eps => const_eps
873  use scale_grid, only: &
874  fdz => grid_fdz, &
875  rcdz => grid_rcdz
876  implicit none
877 
878  real(RP), intent(inout) :: phi(KA)
879 
880  real(RP) :: dev (KA)
881  real(RP) :: flux(KA)
882  real(RP) :: fact(KA)
883 
884  integer, parameter :: iter_max = 100
885  real(RP) :: sig0, sig1, zerosw
886  logical :: updated
887  integer :: k, iter
888  !---------------------------------------------------------------------------
889 
890  dev(ks) = 0.0_rp
891  dev(ke) = 0.0_rp
892 
893  flux(ks-1:ks+1) = 0.0_rp
894  flux(ke-1:ke+1) = 0.0_rp
895 
896  fact(ks-1:ks+1) = 0.0_rp
897  fact(ke-1:ke+1) = 0.0_rp
898 
899  do iter = 1, iter_max
900  updated = .false.
901 
902  do k = ks+1, ke-1
903  dev(k) = phi(k) - ( fdz(k-1)*phi(k+1) + fdz(k)*phi(k-1) ) / ( fdz(k) + fdz(k-1) )
904  enddo
905 
906  do k = ks+2, ke-2
907  sig0 = dev(k) * dev(k-1)
908  sig1 = dev(k) * dev(k+1)
909  ! if (sig0>0 .OR. sig1>0) then flux(k) = 0.0
910  flux(k) = dev(k) &
911  / ( 2.0_rp*rcdz(k) + ( fdz(k-1)*rcdz(k+1) + fdz(k)*rcdz(k-1) ) / ( fdz(k) + fdz(k-1) ) ) &
912  * ( sign(0.5_rp ,sig0) + sign(0.5_rp ,sig1) ) &
913  * ( sign(0.25_rp,sig0) + sign(0.25_rp,sig1) - 0.5_rp )
914  updated = updated .OR. ( sig0 < -eps .AND. sig1 < -eps )
915  enddo
916 
917  sig1 = dev(ks+1) * dev(ks+2)
918  flux(ks+1) = dev(ks+1) &
919  / ( 2.0_rp*rcdz(ks+1) + (fdz(ks)*rcdz(ks+2)+fdz(ks+1)*rcdz(ks))/(fdz(ks+1)+fdz(ks)) ) &
920  * ( 0.5_rp - sign(0.5_rp ,sig1) )
921  updated = updated .OR. ( sig1 < -eps )
922 
923  sig0 = dev(ke-1) * dev(ke-2)
924  flux(ke-1) = dev(ke-1) &
925  / ( 2.0_rp*rcdz(ke-1) + (fdz(ke-2)*rcdz(ke)+fdz(ke-1)*rcdz(ke-2))/(fdz(ke-1)+fdz(ke-2)) ) &
926  * ( 0.5_rp - sign(0.5_rp ,sig0) )
927  updated = updated .OR. ( sig0 < -eps )
928 
929  if ( .NOT. updated ) exit
930 
931  do k = ks+1, ke-1
932  zerosw = 0.5_rp - sign( 0.5_rp, abs(flux(k))-eps ) ! if flux(k) == 0 then fact(k) = 0.0
933  fact(k) = flux(k) / ( flux(k) - flux(k+1) - flux(k-1) + zerosw )
934  enddo
935 
936  do k = ks, ke
937  phi(k) = phi(k) + ( flux(k+1) * fact(k+1) &
938  - flux(k ) * fact(k ) * 2.0_rp &
939  + flux(k-1) * fact(k-1) ) * rcdz(k)
940  enddo
941 
942  if ( iter == iter_max ) then
943  if( io_l ) write(io_fid_log,*) "*** [scale_atmos_refstate/smoothing] iteration not converged!", phi
944  endif
945  enddo
946 
947  return
948  end subroutine smoothing
949 
950 end module scale_atmos_refstate
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_temp
refernce temperature [K]
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
module ATMOSPHERE / Saturation adjustment
subroutine, public prc_mpistop
Abort MPI.
integer, public jeb
subroutine, public atmos_hydrostatic_buildrho_atmos_rev_3d(dens, temp, pres, pott, qv, qc, dz, kref_in)
Build up density from lowermost atmosphere (3D)
real(dp), public time_nowsec
subday part of current time [sec]
Definition: scale_time.F90:68
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
real(rp), dimension(:), allocatable, public grid_cz
center coordinate [m]: z, local=global
module ATMOSPHERE / Reference state
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pott
refernce potential temperature [K]
module STDIO
Definition: scale_stdio.F90:12
integer, public ke
end point of inner domain: z, local
integer, public qa
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
subroutine, public interp_vertical_xi2z(var, var_Z)
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
subroutine, public atmos_refstate_write
Write reference state profile.
real(rp), dimension(qa_max), public tracer_cv
module FILE I/O (netcdf)
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
real(rp), public const_undef
Definition: scale_const.F90:43
integer, public ieb
module grid index
real(rp), dimension(qa_max), public tracer_cp
logical, public io_nml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:62
module TRACER
module ATMOSPHERE / Typical vertical profile
integer, public ia
of whole cells: x, local, with HALO
logical, public atmos_refstate_update_flag
subroutine, public comm_horizontal_mean(varmean, var)
calculate horizontal mean (global total with communication)
Definition: scale_comm.F90:482
real(rp), dimension(:), allocatable, public grid_fdz
z-length of grid(k+1) to grid(k) [m]
module GRID (real space)
integer, public ka
of whole cells: z, local, with HALO
subroutine, public atmos_refstate_update(DENS, RHOT, QTRC)
Update reference state profile (Horizontal average)
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_dens
refernce density [kg/m3]
subroutine, public atmos_refstate_setup
Setup.
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:90
real(rp), dimension(:), allocatable, public grid_fz
face coordinate [m]: z, local=global
module COMMUNICATION
Definition: scale_comm.F90:23
module ATMOSPHERE / Hydrostatic barance
module TIME
Definition: scale_time.F90:15
module PROCESS
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:71
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pres
refernce pressure [Pa]
module CONSTANT
Definition: scale_const.F90:14
integer, public ks
start point of inner domain: z, local
module GRID (cartesian)
subroutine, public atmos_refstate_read
Read reference state profile.
module profiler
Definition: scale_prof.F90:10
real(rp), public const_eps
small number
Definition: scale_const.F90:36
real(rp), dimension(:,:,:), allocatable, public real_phi
geopotential [m2/s2] (cell center)
module ATMOSPHERE / Thermodynamics
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_qv
refernce vapor [kg/kg]
subroutine, public fileio_open(fid, basename)
open a netCDF file for read
module PRECISION
integer, public isb
subroutine, public fileio_close(fid)
Close a netCDF file.
subroutine, public interp_vertical_z2xi(var, var_Xi)
subroutine atmos_refstate_calc3d
apply 1D reference to 3D (terrain-following) with re-calc hydrostatic balance
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
subroutine, public atmos_refstate_resume(DENS, RHOT, QTRC)
Resume.
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
integer, public jsb
subroutine, public atmos_hydrostatic_buildrho_atmos_rev_2d(dens_L1, temp_L1, pres_L1, pott_L1, qv_L1, qc_L1, dens_L2, pott_L2, qv_L2, qc_L2, dz, k)
Build up density (2D)
integer, public io_fid_nml
Log file ID (only for output namelist)
Definition: scale_stdio.F90:57
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:89
subroutine, public atmos_hydrostatic_buildrho_atmos_0d(dens_L2, temp_L2, pres_L2, pott_L2, qv_L2, qc_L2, dens_L1, pott_L1, qv_L1, qc_L1, dz, k)
Build up density (0D)
real(rp), dimension(qa_max), public tracer_mass
module INTERPOLATION
integer, public ja
of whole cells: y, local, with HALO