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