SCALE-RM
scale_atmos_refstate.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: atmos_refstate_setup
28  public :: atmos_refstate_read
29  public :: atmos_refstate_write
30  public :: atmos_refstate_update
31 
32  !-----------------------------------------------------------------------------
33  !
34  !++ Public parameters & variables
35  !
36  logical, public :: atmos_refstate_update_flag = .false.
37 
38  real(RP), public, allocatable :: atmos_refstate_pres(:,:,:)
39  real(RP), public, allocatable :: atmos_refstate_temp(:,:,:)
40  real(RP), public, allocatable :: atmos_refstate_dens(:,:,:)
41  real(RP), public, allocatable :: atmos_refstate_pott(:,:,:)
42  real(RP), public, allocatable :: atmos_refstate_qv (:,:,:)
43 
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  private :: atmos_refstate_generate_isa
49  private :: atmos_refstate_generate_uniform
50  private :: atmos_refstate_generate_zero
51 
52  !-----------------------------------------------------------------------------
53  !
54  !++ Private parameters & variables
55  !
56  character(len=H_LONG), private :: atmos_refstate_in_basename = ''
57  logical, private :: atmos_refstate_in_check_coordinates = .true.
58  character(len=H_LONG), private :: atmos_refstate_out_basename = ''
59  character(len=H_MID), private :: atmos_refstate_out_title = 'SCALE-RM RefState'
60  character(len=H_SHORT), private :: atmos_refstate_out_dtype = 'DEFAULT'
61 
62  character(len=H_SHORT), private :: atmos_refstate_type = 'UNIFORM'
63  real(RP), private :: atmos_refstate_temp_sfc = 300.0_rp
64  real(RP), private :: atmos_refstate_rh = 0.0_rp
65  real(RP), private :: atmos_refstate_pott_uniform = 300.0_rp
66  real(DP), private :: atmos_refstate_update_dt = -1.0_dp
67 
68  real(DP), private :: last_updated
69 
70  real(RP), private, allocatable :: atmos_refstate1d_pres(:)
71  real(RP), private, allocatable :: atmos_refstate1d_temp(:)
72  real(RP), private, allocatable :: atmos_refstate1d_dens(:)
73  real(RP), private, allocatable :: atmos_refstate1d_pott(:)
74  real(RP), private, allocatable :: atmos_refstate1d_qv (:)
75 
76  !-----------------------------------------------------------------------------
77 contains
78  !-----------------------------------------------------------------------------
80  subroutine atmos_refstate_setup( &
81  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
82  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
83  use scale_const, only: &
84  undef => const_undef
85  use scale_prc, only: &
86  prc_abort
87  implicit none
88  integer, intent(in) :: KA, KS, KE
89  integer, intent(in) :: IA, IS, IE
90  integer, intent(in) :: JA, JS, JE
91 
92  real(RP), intent(in) :: CZ ( ka)
93  real(RP), intent(in) :: FZ (0:ka)
94  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
95  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
96  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
97 
98  namelist / param_atmos_refstate / &
99  atmos_refstate_in_basename, &
100  atmos_refstate_out_basename, &
101  atmos_refstate_out_title, &
102  atmos_refstate_out_dtype, &
103  atmos_refstate_type, &
104  atmos_refstate_temp_sfc, &
105  atmos_refstate_rh, &
106  atmos_refstate_pott_uniform, &
107  atmos_refstate_update_dt
108 
109  integer :: ierr
110  !---------------------------------------------------------------------------
111 
112  log_newline
113  log_info("ATMOS_REFSTATE_setup",*) 'Setup'
114 
115  !--- read namelist
116  rewind(io_fid_conf)
117  read(io_fid_conf,nml=param_atmos_refstate,iostat=ierr)
118  if( ierr < 0 ) then !--- missing
119  log_info("ATMOS_REFSTATE_setup",*) 'Not found namelist. Default used.'
120  elseif( ierr > 0 ) then !--- fatal error
121  log_error("ATMOS_REFSTATE_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_REFSTATE. Check!'
122  call prc_abort
123  endif
124  log_nml(param_atmos_refstate)
125 
126  allocate( atmos_refstate_pres(ka,ia,ja) )
127  allocate( atmos_refstate_temp(ka,ia,ja) )
128  allocate( atmos_refstate_dens(ka,ia,ja) )
129  allocate( atmos_refstate_pott(ka,ia,ja) )
130  allocate( atmos_refstate_qv(ka,ia,ja) )
131  atmos_refstate_pres(:,:,:) = undef
132  atmos_refstate_temp(:,:,:) = undef
133  atmos_refstate_dens(:,:,:) = undef
134  atmos_refstate_pott(:,:,:) = undef
135  atmos_refstate_qv(:,:,:) = undef
136 
137  allocate( atmos_refstate1d_pres(ka) )
138  allocate( atmos_refstate1d_temp(ka) )
139  allocate( atmos_refstate1d_dens(ka) )
140  allocate( atmos_refstate1d_pott(ka) )
141  allocate( atmos_refstate1d_qv(ka) )
142  atmos_refstate1d_pres(:) = undef
143  atmos_refstate1d_temp(:) = undef
144  atmos_refstate1d_dens(:) = undef
145  atmos_refstate1d_pott(:) = undef
146  atmos_refstate1d_qv(:) = undef
147 
148  log_newline
149  log_info("ATMOS_REFSTATE_setup",*) 'Reference state settings'
150  if ( atmos_refstate_in_basename /= '' ) then
151  log_info_cont(*) 'Input file of reference state : ', trim(atmos_refstate_in_basename)
152  else
153  log_info_cont(*) 'Input file of reference state : Nothing, generate internally'
154  endif
155 
156  if ( atmos_refstate_out_basename /= '' ) then
157  log_info_cont(*) 'Output file of reference state : ', trim(atmos_refstate_out_basename)
158  else
159  log_info_cont(*) 'Output file of reference state : No output'
160  endif
161 
162  ! input or generate reference profile
163  if ( atmos_refstate_in_basename /= '' ) then
164 
165  call atmos_refstate_read( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
166  cz(:), fz(:), & ! [IN]
167  real_cz(:,:,:), real_fz(:,:,:), & ! [IN]
168  real_phi(:,:,:) ) ! [IN]
169 
170  else
171  if ( atmos_refstate_type == 'ISA' ) then
172 
173  log_info_cont(*) 'Reference type : ISA'
174  log_info_cont(*) 'Surface temperature [K] : ', atmos_refstate_temp_sfc
175  log_info_cont(*) 'Surface & environment RH [%] : ', atmos_refstate_rh
176 
177  call atmos_refstate_generate_isa( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
178  cz(:), fz(:), & ! [IN]
179  real_cz(:,:,:), real_fz(:,:,:), & ! [IN]
180  real_phi(:,:,:) ) ! [IN]
181 
182  elseif ( atmos_refstate_type == 'UNIFORM' ) then
183 
184  log_info_cont(*) 'Reference type : UNIFORM POTT'
185  log_info_cont(*) 'Potential temperature : ', atmos_refstate_pott_uniform
186 
187  call atmos_refstate_generate_uniform( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
188  cz(:), fz(:), & ! [IN]
189  real_cz(:,:,:), real_fz(:,:,:), & ! [IN]
190  real_phi(:,:,:) ) ! [IN]
191 
192  elseif ( atmos_refstate_type == 'ZERO' ) then
193 
194  log_info_cont(*) 'Reference type : ZERO'
195 
196  call atmos_refstate_generate_zero( ka, ia, ja )
197 
198  elseif ( atmos_refstate_type == 'INIT' ) then
199 
200  if ( atmos_refstate_update_dt >= 0.0_rp ) then
202  endif
203 
204  log_info_cont(*) 'Reference type : Generate from initial data'
205  log_info_cont(*) 'Update state? : ', atmos_refstate_update_flag
206  log_info_cont(*) 'Update interval [sec] : ', atmos_refstate_update_dt
207 
208  else
209  log_error("ATMOS_REFSTATE_setup",*) 'ATMOS_REFSTATE_TYPE must be "ISA", "UNIFORM", "ZERO", or "INIT". Check! : ', trim(atmos_refstate_type)
210  call prc_abort
211  endif
212 
213  endif
214 
215  return
216  end subroutine atmos_refstate_setup
217 
218  !-----------------------------------------------------------------------------
220  subroutine atmos_refstate_read( &
221  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
222  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
223  use scale_file_cartesc, only: &
225  file_cartesc_check_coordinates, &
226  file_cartesc_read, &
228  use scale_prc, only: &
229  prc_abort
230  implicit none
231 
232  integer, intent(in) :: KA, KS, KE
233  integer, intent(in) :: IA, IS, IE
234  integer, intent(in) :: JA, JS, JE
235  real(RP), intent(in) :: CZ ( ka)
236  real(RP), intent(in) :: FZ (0:ka)
237  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
238  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
239  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
240 
241  integer :: fid
242  !---------------------------------------------------------------------------
243 
244  log_newline
245  log_info("ATMOS_REFSTATE_read",*) 'Input reference state profile '
246 
247  if ( atmos_refstate_in_basename /= '' ) then
248 
249  call file_cartesc_open( atmos_refstate_in_basename, fid )
250 
251  if ( atmos_refstate_in_check_coordinates ) then
252  call file_cartesc_check_coordinates( fid, atmos=.true. )
253  end if
254 
255  ! 1D
256  call file_cartesc_read( fid, 'PRES_ref', 'Z', atmos_refstate1d_pres(:) )
257  call file_cartesc_read( fid, 'TEMP_ref', 'Z', atmos_refstate1d_temp(:) )
258  call file_cartesc_read( fid, 'DENS_ref', 'Z', atmos_refstate1d_dens(:) )
259  call file_cartesc_read( fid, 'POTT_ref', 'Z', atmos_refstate1d_pott(:) )
260  call file_cartesc_read( fid, 'QV_ref', 'Z', atmos_refstate1d_qv(:) )
261 
262  ! 3D
263  call file_cartesc_read( fid, 'PRES_ref3D', 'ZXY', atmos_refstate_pres(:,:,:) )
264  call file_cartesc_read( fid, 'TEMP_ref3D', 'ZXY', atmos_refstate_temp(:,:,:) )
265  call file_cartesc_read( fid, 'DENS_ref3D', 'ZXY', atmos_refstate_dens(:,:,:) )
266  call file_cartesc_read( fid, 'POTT_ref3D', 'ZXY', atmos_refstate_pott(:,:,:) )
267  call file_cartesc_read( fid, 'QV_ref3D', 'ZXY', atmos_refstate_qv(:,:,:) )
268 
269  else
270  log_error("ATMOS_REFSTATE_read",*) 'refstate file is not specified.'
271  call prc_abort
272  endif
273 
274  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
275  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
276 
277  return
278  end subroutine atmos_refstate_read
279 
280  !-----------------------------------------------------------------------------
282  subroutine atmos_refstate_write
283  use scale_file_cartesc, only: &
284  file_cartesc_write
285  implicit none
286 
287  logical, save :: first = .true.
288  !---------------------------------------------------------------------------
289 
290  if ( .not. first ) return
291  first = .false.
292 
293  if ( atmos_refstate_out_basename /= '' ) then
294 
295  log_newline
296  log_info("ATMOS_REFSTATE_write",*) 'Output reference state profile '
297 
298  ! 1D
299  call file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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 file_cartesc_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  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
331  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
332  use scale_const, only: &
333  epsvap => const_epsvap, &
334  pstd => const_pstd
335  use scale_atmos_profile, only: &
336  profile_isa => atmos_profile_isa
337  use scale_atmos_hydrostatic, only: &
338  hydrostatic_buildrho => atmos_hydrostatic_buildrho
339  use scale_atmos_saturation, only: &
340  saturation_psat_all => atmos_saturation_psat_all, &
341  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
342  implicit none
343  integer, intent(in) :: KA, KS, KE
344  integer, intent(in) :: IA, IS, IE
345  integer, intent(in) :: JA, JS, JE
346 
347  real(RP), intent(in) :: CZ ( ka)
348  real(RP), intent(in) :: FZ (0:ka)
349  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
350  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
351  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
352 
353  real(RP) :: temp(ka)
354  real(RP) :: pres(ka)
355  real(RP) :: dens(ka)
356  real(RP) :: pott(ka)
357  real(RP) :: qv (ka)
358  real(RP) :: qc (ka)
359 
360  real(RP) :: temp_sfc
361  real(RP) :: pres_sfc
362  real(RP) :: pott_sfc
363  real(RP) :: qv_sfc
364  real(RP) :: qc_sfc
365 
366  real(RP) :: qsat(ka)
367  real(RP) :: psat_sfc
368 
369  integer :: k
370  !---------------------------------------------------------------------------
371 
372  pott_sfc = atmos_refstate_temp_sfc
373  pres_sfc = pstd
374 
375  call profile_isa( ka, ks, ke, & ! [IN]
376  pott_sfc, pres_sfc, & ! [IN]
377  cz(:), & ! [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( ka, ks, ke, &
387  pott(:), qv(:), qc(:), & ! [IN]
388  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
389  cz(:), fz(:), & ! [IN]
390  dens(:), temp(:), pres(:), & ! [OUT]
391  temp_sfc ) ! [OUT]
392 
393  ! calc QV from RH
394  call saturation_psat_all( temp_sfc, psat_sfc )
395  call saturation_dens2qsat_all( ka, ks, ke, &
396  temp(:), dens(:), & ! [IN]
397  qsat(:) ) ! [OUT]
398 
399  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
400  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp-epsvap) * psat_sfc )
401  do k = ks, ke
402  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
403  enddo
404 
405  ! make density & pressure profile in moist condition
406  call hydrostatic_buildrho( ka, ks, ke, &
407  pott(:), qv(:), qc(:), & ! [IN]
408  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
409  cz(:), fz(:), & ! [IN]
410  dens(:), temp(:), pres(:), & ! [OUT]
411  temp_sfc ) ! [OUT]
412 
413  atmos_refstate1d_pres(:) = pres(:)
414  atmos_refstate1d_temp(:) = temp(:)
415  atmos_refstate1d_dens(:) = dens(:)
416  atmos_refstate1d_pott(:) = pott(:)
417  atmos_refstate1d_qv(:) = qv(:)
418 
419  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
420  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
421 
422  return
423  end subroutine atmos_refstate_generate_isa
424 
425  !-----------------------------------------------------------------------------
427  subroutine atmos_refstate_generate_uniform( &
428  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
429  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
430  use scale_const, only: &
431  epsvap => const_epsvap, &
432  pstd => const_pstd
433  use scale_atmos_hydrostatic, only: &
434  hydrostatic_buildrho => atmos_hydrostatic_buildrho
435  use scale_atmos_saturation, only: &
436  saturation_psat_all => atmos_saturation_psat_all, &
437  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
438  implicit none
439  integer, intent(in) :: KA, KS, KE
440  integer, intent(in) :: IA, IS, IE
441  integer, intent(in) :: JA, JS, JE
442 
443  real(RP), intent(in) :: CZ ( ka)
444  real(RP), intent(in) :: FZ (0:ka)
445  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
446  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
447  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
448 
449  real(RP) :: temp(ka)
450  real(RP) :: pres(ka)
451  real(RP) :: dens(ka)
452  real(RP) :: pott(ka)
453  real(RP) :: qv (ka)
454  real(RP) :: qc (ka)
455 
456  real(RP) :: temp_sfc
457  real(RP) :: pres_sfc
458  real(RP) :: pott_sfc
459  real(RP) :: qv_sfc
460  real(RP) :: qc_sfc
461 
462  real(RP) :: qsat(ka)
463  real(RP) :: psat_sfc
464 
465  integer :: k
466  !---------------------------------------------------------------------------
467 
468  pres_sfc = pstd
469  pott_sfc = atmos_refstate_temp_sfc
470  qv_sfc = 0.0_rp
471  qc_sfc = 0.0_rp
472 
473  do k = 1, ka
474  pott(k) = atmos_refstate_pott_uniform
475  qv(k) = 0.0_rp
476  qc(k) = 0.0_rp
477  enddo
478 
479  ! make density & pressure profile in dry condition
480  call hydrostatic_buildrho( ka, ks, ke, &
481  pott(:), qv(:), qc(:), & ! [IN]
482  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
483  cz(:), fz(:), & ! [IN]
484  dens(:), temp(:), pres(:), & ! [OUT]
485  temp_sfc ) ! [OUT]
486 
487  ! calc QV from RH
488  call saturation_psat_all( temp_sfc, psat_sfc )
489  call saturation_dens2qsat_all( ka, ks, ke, &
490  temp(:), dens(:), & ! [IN]
491  qsat(:) ) ! [OUT]
492 
493  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
494  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp - epsvap) * psat_sfc )
495  do k = ks, ke
496  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
497  enddo
498 
499  ! make density & pressure profile in moist condition
500  call hydrostatic_buildrho( ka, ks, ke, &
501  pott(:), qv(:), qc(:), & ! [IN]
502  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
503  cz(:), fz(:), & ! [IN]
504  dens(:), temp(:), pres(:), & ! [OUT]
505  temp_sfc ) ! [OUT]
506 
507  atmos_refstate1d_pres(:) = pres(:)
508  atmos_refstate1d_temp(:) = temp(:)
509  atmos_refstate1d_dens(:) = dens(:)
510  atmos_refstate1d_pott(:) = pott(:)
511  atmos_refstate1d_qv(:) = qv(:)
512 
513  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
514  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
515 
516  return
517  end subroutine atmos_refstate_generate_uniform
518 
519  !-----------------------------------------------------------------------------
521  subroutine atmos_refstate_generate_zero( KA, IA, JA )
522  implicit none
523  integer, intent(in) :: KA, IA, JA
524 
525  integer :: k, i, j
526  !---------------------------------------------------------------------------
527 
528  do k = 1, ka
529  atmos_refstate1d_pres(k) = 0.0_rp
530  atmos_refstate1d_temp(k) = 0.0_rp
531  atmos_refstate1d_dens(k) = 0.0_rp
532  atmos_refstate1d_pott(k) = 0.0_rp
533  atmos_refstate1d_qv(k) = 0.0_rp
534  end do
535 
536  do j = 1, ja
537  do i = 1, ia
538  do k = 1, ka
539  atmos_refstate_dens(k,i,j) = 0.0_rp
540  atmos_refstate_temp(k,i,j) = 0.0_rp
541  atmos_refstate_pres(k,i,j) = 0.0_rp
542  atmos_refstate_pott(k,i,j) = 0.0_rp
543  atmos_refstate_qv(k,i,j) = 0.0_rp
544  enddo
545  enddo
546  enddo
547 
548  return
549  end subroutine atmos_refstate_generate_zero
550 
551  !-----------------------------------------------------------------------------
553  subroutine atmos_refstate_update( &
554  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
555  DENS, POTT, TEMP, PRES, QV, &
556  CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA, &
557  nowsec, &
558  force )
559  use scale_statistics, only: &
560  statistics_horizontal_mean
561  use scale_interp_vert, only: &
563  implicit none
564 
565  integer, intent(in) :: KA, KS, KE
566  integer, intent(in) :: IA, IS, IE
567  integer, intent(in) :: JA, JS, JE
568  real(RP), intent(in) :: DENS (ka,ia,ja)
569  real(RP), intent(in) :: POTT (ka,ia,ja)
570  real(RP), intent(in) :: TEMP (ka,ia,ja)
571  real(RP), intent(in) :: PRES (ka,ia,ja)
572  real(RP), intent(in) :: QV (ka,ia,ja)
573  real(RP), intent(in) :: CZ ( ka)
574  real(RP), intent(in) :: FZ (0:ka)
575  real(RP), intent(in) :: FDZ ( ka-1)
576  real(RP), intent(in) :: RCDZ ( ka)
577  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
578  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
579  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
580  real(RP), intent(in) :: AREA ( ia,ja)
581  real(DP), intent(in) :: nowsec
582  logical, intent(in), optional :: force
583 
584  real(RP) :: work(ka,ia,ja)
585  logical :: force_
586 
587  integer :: k
588  !---------------------------------------------------------------------------
589 
590  if ( present(force) ) then
591  force_ = force
592  else
593  force_ = .false.
594  end if
595 
596  if ( force_ .or. ( nowsec - last_updated >= atmos_refstate_update_dt ) ) then
597 
598  log_info("ATMOS_REFSTATE_update",*) 'update reference state'
599 
600  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
601  cz(:), real_cz(:,:,:), temp(:,:,:), & ! [IN]
602  work(:,:,:) ) ! [OUT]
603  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
604  work(:,:,:), area(:,:), & ! [IN]
605  atmos_refstate1d_temp(:) ) ! [IN]
606 
607  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
608  cz(:), real_cz(:,:,:), pres(:,:,:), & ! [IN]
609  work(:,:,:) ) ! [OUT]
610  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
611  work(:,:,:), area(:,:), & ! [IN]
612  atmos_refstate1d_pres(:) ) ! [OUT]
613 
614  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
615  cz(:), real_cz(:,:,:), dens(:,:,:), & ! [IN]
616  work(:,:,:) ) ! [OUT]
617  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
618  work(:,:,:), area(:,:), & ! [IN]
619  atmos_refstate1d_dens(:) ) ! [OUT]
620 
621  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
622  cz(:), real_cz(:,:,:), pott(:,:,:), & ! [IN]
623  work(:,:,:) ) ! [OUT]
624  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
625  work(:,:,:), area(:,:), & ! [IN]
626  atmos_refstate1d_pott(:) ) ! [OUT]
627 
628  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
629  cz(:), real_cz(:,:,:), qv(:,:,:), & ! [IN]
630  work(:,:,:) ) ! [OUT]
631  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
632  work(:,:,:), area(:,:), & ! [IN]
633  atmos_refstate1d_qv(:) ) ! [OUT]
634 
635  do k = ke-1, ks, -1 ! fill undefined value
636  if( atmos_refstate1d_dens(k) <= 0.0_rp ) atmos_refstate1d_dens(k) = atmos_refstate1d_dens(k+1)
637  if( atmos_refstate1d_temp(k) <= 0.0_rp ) atmos_refstate1d_temp(k) = atmos_refstate1d_temp(k+1)
638  if( atmos_refstate1d_pres(k) <= 0.0_rp ) atmos_refstate1d_pres(k) = atmos_refstate1d_pres(k+1)
639  if( atmos_refstate1d_pott(k) <= 0.0_rp ) atmos_refstate1d_pott(k) = atmos_refstate1d_pott(k+1)
640  if( atmos_refstate1d_qv(k) <= 0.0_rp ) atmos_refstate1d_qv(k) = atmos_refstate1d_qv(k+1)
641  enddo
642 
643  call atmos_refstate_smoothing( ka, ks, ke, fdz(:), rcdz(:), atmos_refstate1d_pott(:) )
644  call atmos_refstate_smoothing( ka, ks, ke, fdz(:), rcdz(:), atmos_refstate1d_qv(:) )
645 
646  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
647  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
648 
649  last_updated = nowsec
650 
651 
652  log_newline
653  log_info("ATMOS_REFSTATE_update",*) 'Generated reference state of atmosphere:'
654  log_info_cont(*) '============================================================================='
655  log_info_cont(*) ' z*-coord.: pressure: temperature: density: pot.temp.: water vapor'
656  do k = ks, ke
657  log_info_cont('(6F13.5)') cz(k), &
658  atmos_refstate1d_pres(k), &
659  atmos_refstate1d_temp(k), &
660  atmos_refstate1d_dens(k), &
661  atmos_refstate1d_pott(k), &
662  atmos_refstate1d_qv(k)
663  enddo
664  log_info_cont(*) '============================================================================='
665 
666  ! output reference profile
668 
669  endif
670 
671  return
672  end subroutine atmos_refstate_update
673 
674  !-----------------------------------------------------------------------------
676  subroutine atmos_refstate_calc3d( &
677  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
678  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
679  use scale_const, only: &
680  undef => const_undef, &
681  rdry => const_rdry, &
682  cpdry => const_cpdry, &
683  p00 => const_pre00
684  use scale_comm_cartesc, only: &
685  comm_vars8, &
686  comm_wait
687  use scale_interp_vert, only: &
689  use scale_atmos_hydrostatic, only: &
690  hydrostatic_buildrho_atmos => atmos_hydrostatic_buildrho_atmos, &
691  hydrostatic_buildrho_atmos_rev => atmos_hydrostatic_buildrho_atmos_rev
692  implicit none
693  integer, intent(in) :: KA, KS, KE
694  integer, intent(in) :: IA, IS, IE
695  integer, intent(in) :: JA, JS, JE
696 
697  real(RP), intent(in) :: CZ ( ka)
698  real(RP), intent(in) :: FZ (0:ka)
699  real(RP), intent(in) :: REAL_CZ ( ka,ia,ja)
700  real(RP), intent(in) :: REAL_FZ (0:ka,ia,ja)
701  real(RP), intent(in) :: REAL_PHI( ka,ia,ja)
702 
703  real(RP) :: dens(ka,ia,ja)
704  real(RP) :: temp(ka,ia,ja)
705  real(RP) :: pres(ka,ia,ja)
706  real(RP) :: pott(ka,ia,ja)
707  real(RP) :: qv (ka,ia,ja)
708  real(RP) :: qc (ka,ia,ja)
709  real(RP) :: dz (ka,ia,ja)
710 
711  real(RP) :: dens_toa_1D
712  real(RP) :: temp_toa_1D
713  real(RP) :: pres_toa_1D
714  real(RP) :: qc_1D
715  real(RP) :: dz_1D
716 
717  real(RP) :: work(ka,ia,ja)
718  real(RP) :: RovCP
719  integer :: k, i, j
720  !---------------------------------------------------------------------------
721 
722  rovcp = rdry / cpdry
723 
724  !--- potential temperature
725  do j = js, je
726  do i = is, ie
727  work(:,i,j) = atmos_refstate1d_pott(:)
728  enddo
729  enddo
730 
731  call interp_vert_z2xi( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
732  real_cz(:,:,:), cz(:), & ! [IN]
733  work(:,:,:), & ! [IN]
734  pott(:,:,:) ) ! [OUT]
735 
736  !--- water vapor
737  do j = js, je
738  do i = is, ie
739  work(:,i,j) = atmos_refstate1d_qv(:)
740  enddo
741  enddo
742 
743  call interp_vert_z2xi( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
744  real_cz(:,:,:), cz(:), & ! [IN]
745  work(:,:,:), & ! [IN]
746  qv(:,:,:) ) ! [OUT]
747 
748 
749  !--- build up density to TOA (1D)
750  qc_1d = 0.0_rp
751  dz_1d = fz(ke) - cz(ke)
752 
753  call hydrostatic_buildrho_atmos( atmos_refstate1d_pott(ke), & ! [IN]
754  atmos_refstate1d_qv(ke), & ! [IN]
755  qc_1d, & ! [IN]
756  atmos_refstate1d_dens(ke), & ! [IN]
757  atmos_refstate1d_pott(ke), & ! [IN]
758  atmos_refstate1d_qv(ke), & ! [IN]
759  qc_1d, & ! [IN]
760  dz_1d, & ! [IN]
761  ke+1, & ! [IN]
762  dens_toa_1d, & ! [OUT]
763  temp_toa_1d, & ! [OUT]
764  pres_toa_1d ) ! [OUT]
765 
766  ! build down density from TOA (3D)
767  do j = js, je
768  do i = is, ie
769  dz(ks-1,i,j) = real_cz(ks,i,j) - real_fz(ks-1,i,j) ! distance from surface to cell center
770  do k = ks, ke-1
771  dz(k,i,j) = real_cz(k+1,i,j) - real_cz(k,i,j) ! distance from cell center to cell center
772  enddo
773  dz(ke,i,j) = real_fz(ke,i,j) - real_cz(ke,i,j) ! distance from cell center to TOA
774  enddo
775  enddo
776 
777  do j = js, je
778  do i = is, ie
779  dens(ke+1,i,j) = dens_toa_1d
780  temp(ke+1,i,j) = temp_toa_1d
781  pres(ke+1,i,j) = pres_toa_1d
782  pott(ke+1,i,j) = pott(ke,i,j)
783  qv(ke+1,i,j) = qv(ke,i,j)
784  enddo
785  enddo
786 
787  do j = js, je
788  do i = is, ie
789  pott(ks-1,i,j) = pott(ks,i,j)
790  qv(ks-1,i,j) = qv(ks,i,j)
791  enddo
792  enddo
793 
794  qc(:,:,:) = 0.0_rp
795 
796  call hydrostatic_buildrho_atmos_rev( ia, is, ie, ja, js, je, &
797  pott(ke,:,:), qv(ke,:,:), qc(ke,:,:), & ! [IN]
798  dens(ke+1,:,:), pott(ke+1,:,:), qv(ke+1,:,:), qc(ke+1,:,:), & ! [IN]
799  dz(ke,:,:), ke+1, & ! [IN]
800  dens(ke ,:,:), temp(ke ,:,:), pres(ke ,:,:) ) ! [OUT]
801 
802  call hydrostatic_buildrho_atmos_rev( ka, ks, ke, ia, is, ie, ja, js, je, &
803  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
804  dz(:,:,:), & ! [IN]
805  dens(:,:,:), & ! [INOUT]
806  temp(:,:,:), pres(:,:,:) ) ! [OUT]
807 
808  do j = js, je
809  do i = is, ie
810  do k = ks, ke
811  atmos_refstate_dens(k,i,j) = dens(k,i,j)
812  atmos_refstate_temp(k,i,j) = temp(k,i,j)
813  atmos_refstate_pres(k,i,j) = pres(k,i,j)
814  atmos_refstate_pott(k,i,j) = pott(k,i,j)
815  atmos_refstate_qv(k,i,j) = qv(k,i,j)
816  enddo
817  enddo
818  enddo
819 
820  ! boundary condition
821  do j = js, je
822  do i = is, ie
823 
824  atmos_refstate_temp(1:ks-1, i,j) = temp(ks,i,j)
825  atmos_refstate_temp(ke+1:ka,i,j) = temp_toa_1d
826 
827  atmos_refstate_qv(1:ks-1, i,j) = qv(ks,i,j)
828  atmos_refstate_qv(ke+1:ka,i,j) = qv(ke,i,j)
829 
830  atmos_refstate_pres(1:ks-2, i,j) = undef
831  atmos_refstate_pres(ks-1, i,j) = atmos_refstate_pres(ks+1,i,j) &
832  - atmos_refstate_dens(ks ,i,j) * ( real_phi(ks-1,i,j) - real_phi(ks+1,i,j) )
833  atmos_refstate_pres(ke+1, i,j) = atmos_refstate_pres(ke-1,i,j) &
834  - atmos_refstate_dens(ke ,i,j) * ( real_phi(ke+1,i,j) - real_phi(ke-1,i,j) )
835  atmos_refstate_pres(ke+2:ka,i,j) = undef
836 
837  atmos_refstate_dens(1:ks-2, i,j) = undef
838  atmos_refstate_dens(ks-1, i,j) = atmos_refstate_pres(ks-1,i,j) / ( atmos_refstate_temp(ks-1,i,j) * rdry )
839  atmos_refstate_dens(ke+1, i,j) = atmos_refstate_pres(ke+1,i,j) / ( atmos_refstate_temp(ke+1,i,j) * rdry )
840  atmos_refstate_dens(ke+2:ka,i,j) = undef
841 
842  atmos_refstate_pott(1:ks-2, i,j) = undef
843  atmos_refstate_pott(ks-1, i,j) = atmos_refstate_temp(ks-1,i,j) * ( p00 / atmos_refstate_pres(ks-1,i,j) )**rovcp
844  atmos_refstate_pott(ke+1, i,j) = atmos_refstate_temp(ke+1,i,j) * ( p00 / atmos_refstate_pres(ke+1,i,j) )**rovcp
845  atmos_refstate_pott(ke+2:ka,i,j) = undef
846  enddo
847  enddo
848 
849  call comm_vars8( atmos_refstate_dens(:,:,:), 1 )
850  call comm_vars8( atmos_refstate_temp(:,:,:), 2 )
851  call comm_vars8( atmos_refstate_pres(:,:,:), 3 )
852  call comm_vars8( atmos_refstate_pott(:,:,:), 4 )
853  call comm_vars8( atmos_refstate_qv(:,:,:), 5 )
854  call comm_wait ( atmos_refstate_dens(:,:,:), 1, .false. )
855  call comm_wait ( atmos_refstate_temp(:,:,:), 2, .false. )
856  call comm_wait ( atmos_refstate_pres(:,:,:), 3, .false. )
857  call comm_wait ( atmos_refstate_pott(:,:,:), 4, .false. )
858  call comm_wait ( atmos_refstate_qv(:,:,:), 5, .false. )
859 
860  return
861  end subroutine atmos_refstate_calc3d
862 
863  !-----------------------------------------------------------------------------
864  subroutine atmos_refstate_smoothing( &
865  KA, KS, KE, &
866  FDZ, RCDZ, &
867  phi )
868  use scale_const, only: &
869  eps => const_eps
870  implicit none
871  integer, intent(in) :: KA, KS, KE
872 
873  real(RP), intent(in) :: FDZ (ka-1)
874  real(RP), intent(in) :: RCDZ(ka)
875 
876  real(RP), intent(inout) :: phi(ka)
877 
878  real(RP) :: dev (ka)
879  real(RP) :: correction(ka)
880  real(RP) :: fact(ka)
881 
882  integer, parameter :: iter_max = 100
883  real(RP) :: sig0, sig1, zerosw
884  logical :: updated
885  integer :: k, iter
886  !---------------------------------------------------------------------------
887 
888  dev(ks) = 0.0_rp
889  dev(ke) = 0.0_rp
890 
891  correction(ks-1:ks+1) = 0.0_rp
892  correction(ke-1:ke+1) = 0.0_rp
893 
894  fact(ks-1:ks+1) = 0.0_rp
895  fact(ke-1:ke+1) = 0.0_rp
896 
897  do iter = 1, iter_max
898  updated = .false.
899 
900  do k = ks+1, ke-1
901  dev(k) = phi(k) - ( fdz(k-1)*phi(k+1) + fdz(k)*phi(k-1) ) / ( fdz(k) + fdz(k-1) )
902  enddo
903 
904  do k = ks+2, ke-2
905  sig0 = dev(k) * dev(k-1)
906  sig1 = dev(k) * dev(k+1)
907  if ( sig0 < -eps .and. sig1 < -eps ) then
908  correction(k) = dev(k) &
909  / ( 2.0_rp*rcdz(k) + ( fdz(k-1)*rcdz(k+1) + fdz(k)*rcdz(k-1) ) / ( fdz(k) + fdz(k-1) ) )
910  updated = .true.
911  else
912  correction(k) = 0.0_rp
913  end if
914  enddo
915 
916  sig1 = dev(ks+1) * dev(ks+2)
917  if ( sig1 < -eps ) then
918  correction(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  updated = .true.
921  else
922  correction(ks+1) = 0.0_rp
923  end if
924 
925  sig0 = dev(ke-1) * dev(ke-2)
926  if ( sig0 < -eps ) then
927  correction(ke-1) = dev(ke-1) &
928  / ( 2.0_rp*rcdz(ke-1) + (fdz(ke-2)*rcdz(ke)+fdz(ke-1)*rcdz(ke-2))/(fdz(ke-1)+fdz(ke-2)) )
929  updated = .true.
930  else
931  correction(ke-1) = 0.0_rp
932  end if
933 
934  if ( .NOT. updated ) exit
935 
936  do k = ks+1, ke-1
937  zerosw = 0.5_rp - sign( 0.5_rp, abs(correction(k))-eps ) ! if correction(k) == 0 then fact(k) = 0.0
938  fact(k) = correction(k) / ( correction(k) - correction(k+1) - correction(k-1) + zerosw )
939  enddo
940 
941  do k = ks, ke
942  phi(k) = phi(k) + ( correction(k+1) * fact(k+1) &
943  - correction(k ) * fact(k ) * 2.0_rp &
944  + correction(k-1) * fact(k-1) ) * rcdz(k)
945  enddo
946 
947  if ( iter == iter_max ) then
948  log_info("ATMOS_REFSTATE_smoothing",*) "iteration not converged!", phi
949  endif
950  enddo
951 
952  return
953  end subroutine atmos_refstate_smoothing
954 
955 end module scale_atmos_refstate
subroutine, public atmos_refstate_update(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, POTT, TEMP, PRES, QV, CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA, nowsec, force)
Update reference state profile (Horizontal average)
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_temp
refernce temperature [K]
subroutine, public interp_vert_z2xi(KA, KS, KE, IA, IS, IE, JA, JS, JE, Z, Xi, var, var_Xi)
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
module atmosphere / saturation
module atmosphere / reference state
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pott
refernce potential temperature [K]
module INTERPOLATION vertical
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
subroutine, public interp_vert_xi2z(KA, KS, KE, IA, IS, IE, JA, JS, JE, Xi, Z, var, var_Z)
subroutine, public atmos_refstate_write
Write reference state profile.
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
real(rp), public const_undef
Definition: scale_const.F90:41
module COMMUNICATION
module atmosphere / vertical profile
logical, public atmos_refstate_update_flag
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_dens
refernce density [kg/m3]
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
module PROCESS
Definition: scale_prc.F90:11
module atmosphere / hydrostatic barance
subroutine, public atmos_refstate_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI)
Setup.
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:69
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pres
refernce pressure [Pa]
module CONSTANT
Definition: scale_const.F90:11
module profiler
Definition: scale_prof.F90:11
real(rp), public const_eps
small number
Definition: scale_const.F90:33
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_qv
refernce vapor [kg/kg]
subroutine, public atmos_refstate_read(KA, KS, KE, IA, IS, IE, JA, JS, JE, CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI)
Read reference state profile.
module PRECISION
module file / cartesianC
module Statistics
module STDIO
Definition: scale_io.F90:10
subroutine, public file_cartesc_open(basename, fid, aggregate)
open a netCDF file for read
subroutine atmos_refstate_calc3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI)
apply 1D reference to 3D (terrain-following) with re-calc hydrostatic balance
subroutine, public file_cartesc_close(fid)
Close a netCDF file.
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:87
subroutine atmos_refstate_smoothing(KA, KS, KE, FDZ, RCDZ, phi)