SCALE-RM
scale_atmos_refstate.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_refstate
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_finalize
29  public :: atmos_refstate_read
30  public :: atmos_refstate_write
31  public :: atmos_refstate_update
32 
33  !-----------------------------------------------------------------------------
34  !
35  !++ Public parameters & variables
36  !
37  real(rp), public, allocatable :: atmos_refstate_pres(:,:,:)
38  real(rp), public, allocatable :: atmos_refstate_temp(:,:,:)
39  real(rp), public, allocatable :: atmos_refstate_dens(:,:,:)
40  real(rp), public, allocatable :: atmos_refstate_pott(:,:,:)
41  real(rp), public, allocatable :: atmos_refstate_qv (:,:,:)
42 
43  !-----------------------------------------------------------------------------
44  !
45  !++ Private procedure
46  !
47  private :: atmos_refstate_generate_isa
48  private :: atmos_refstate_generate_uniform
49  private :: atmos_refstate_generate_zero
50 
51  !-----------------------------------------------------------------------------
52  !
53  !++ Private parameters & variables
54  !
55  character(len=H_LONG), private :: atmos_refstate_in_basename = ''
56  logical, private :: atmos_refstate_in_check_coordinates = .true.
57  character(len=H_LONG), private :: atmos_refstate_out_basename = ''
58  character(len=H_MID), private :: atmos_refstate_out_title = 'SCALE-RM RefState'
59  character(len=H_SHORT), private :: atmos_refstate_out_dtype = 'DEFAULT'
60 
61  character(len=H_SHORT), private :: atmos_refstate_type = 'UNIFORM'
62  real(rp), private :: atmos_refstate_temp_sfc = 300.0_rp
63  real(rp), private :: atmos_refstate_rh = 0.0_rp
64  real(rp), private :: atmos_refstate_pott_uniform = 300.0_rp
65  real(dp), private :: atmos_refstate_update_dt = -1.0_dp
66  logical, private :: atmos_refstate_update_flag = .false.
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  !$acc enter data create(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
137 
138  allocate( atmos_refstate1d_pres(ka) )
139  allocate( atmos_refstate1d_temp(ka) )
140  allocate( atmos_refstate1d_dens(ka) )
141  allocate( atmos_refstate1d_pott(ka) )
142  allocate( atmos_refstate1d_qv(ka) )
143  atmos_refstate1d_pres(:) = undef
144  atmos_refstate1d_temp(:) = undef
145  atmos_refstate1d_dens(:) = undef
146  atmos_refstate1d_pott(:) = undef
147  atmos_refstate1d_qv(:) = undef
148  !$acc enter data create(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
149 
150  log_newline
151  log_info("ATMOS_REFSTATE_setup",*) 'Reference state settings'
152  if ( atmos_refstate_in_basename /= '' ) then
153  log_info_cont(*) 'Input file of reference state : ', trim(atmos_refstate_in_basename)
154  else
155  log_info_cont(*) 'Input file of reference state : Nothing, generate internally'
156  endif
157 
158  if ( atmos_refstate_out_basename /= '' ) then
159  log_info_cont(*) 'Output file of reference state : ', trim(atmos_refstate_out_basename)
160  else
161  log_info_cont(*) 'Output file of reference state : No output'
162  endif
163 
164  atmos_refstate_update_flag = .false.
165 
166  ! input or generate reference profile
167  if ( atmos_refstate_in_basename /= '' ) then
168 
169  call atmos_refstate_read( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
170  real_phi(:,:,:) ) ! [IN]
171 
172  else
173  if ( atmos_refstate_type == 'ISA' ) then
174 
175  log_info_cont(*) 'Reference type : ISA'
176  log_info_cont(*) 'Surface temperature [K] : ', atmos_refstate_temp_sfc
177  log_info_cont(*) 'Surface & environment RH [%] : ', atmos_refstate_rh
178 
179  call atmos_refstate_generate_isa( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
180  cz(:), fz(:), & ! [IN]
181  real_cz(:,:,:), real_fz(:,:,:), & ! [IN]
182  real_phi(:,:,:) ) ! [IN]
183 
184  elseif ( atmos_refstate_type == 'UNIFORM' ) then
185 
186  log_info_cont(*) 'Reference type : UNIFORM POTT'
187  log_info_cont(*) 'Potential temperature : ', atmos_refstate_pott_uniform
188 
189  call atmos_refstate_generate_uniform( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
190  cz(:), fz(:), & ! [IN]
191  real_cz(:,:,:), real_fz(:,:,:), & ! [IN]
192  real_phi(:,:,:) ) ! [IN]
193 
194  elseif ( atmos_refstate_type == 'ZERO' ) then
195 
196  log_info_cont(*) 'Reference type : ZERO'
197 
198  call atmos_refstate_generate_zero( ka, ia, ja )
199 
200  elseif ( atmos_refstate_type == 'INIT' ) then
201 
202  atmos_refstate_update_flag = .true.
203 
204  log_info_cont(*) 'Reference type : Generate from initial data'
205  log_info_cont(*) 'Update interval [sec] : ', atmos_refstate_update_dt
206 
207  else
208  log_error("ATMOS_REFSTATE_setup",*) 'ATMOS_REFSTATE_TYPE must be "ISA", "UNIFORM", "ZERO", or "INIT". Check! : ', trim(atmos_refstate_type)
209  call prc_abort
210  endif
211 
212  endif
213 
214  last_updated = -1.0_rp
215 
216  return
217  end subroutine atmos_refstate_setup
218 
219  subroutine atmos_refstate_finalize
220 
221  !$acc exit data delete(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
222  !$acc exit data delete(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
223 
224  deallocate( atmos_refstate_pres )
225  deallocate( atmos_refstate_temp )
226  deallocate( atmos_refstate_dens )
227  deallocate( atmos_refstate_pott )
228  deallocate( atmos_refstate_qv )
229 
230  deallocate( atmos_refstate1d_pres )
231  deallocate( atmos_refstate1d_temp )
232  deallocate( atmos_refstate1d_dens )
233  deallocate( atmos_refstate1d_pott )
234  deallocate( atmos_refstate1d_qv )
235 
236  return
237  end subroutine atmos_refstate_finalize
238 
239  !-----------------------------------------------------------------------------
241  subroutine atmos_refstate_read( &
242  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
243  REAL_PHI )
244  use scale_file_cartesc, only: &
245  file_cartesc_open, &
246  file_cartesc_check_coordinates, &
247  file_cartesc_read, &
248  file_cartesc_close
249  use scale_prc, only: &
250  prc_abort
251  implicit none
252  integer, intent(in) :: ka, ks, ke
253  integer, intent(in) :: ia, is, ie
254  integer, intent(in) :: ja, js, je
255  real(rp), intent(in) :: real_phi(ka,ia,ja)
256 
257  integer :: fid
258  !---------------------------------------------------------------------------
259 
260  log_newline
261  log_info("ATMOS_REFSTATE_read",*) 'Input reference state profile '
262 
263  if ( atmos_refstate_in_basename /= '' ) then
264 
265  ! 1D
266  call file_cartesc_open( atmos_refstate_in_basename, fid, single=.true. )
267  call file_cartesc_read( fid, 'PRES_ref', 'Z', atmos_refstate1d_pres(:) )
268  call file_cartesc_read( fid, 'TEMP_ref', 'Z', atmos_refstate1d_temp(:) )
269  call file_cartesc_read( fid, 'DENS_ref', 'Z', atmos_refstate1d_dens(:) )
270  call file_cartesc_read( fid, 'POTT_ref', 'Z', atmos_refstate1d_pott(:) )
271  call file_cartesc_read( fid, 'QV_ref', 'Z', atmos_refstate1d_qv(:) )
272  call file_cartesc_close( fid )
273 
274  ! 3D
275  call file_cartesc_open( atmos_refstate_in_basename, fid )
276  if ( atmos_refstate_in_check_coordinates ) then
277  call file_cartesc_check_coordinates( fid, atmos=.true. )
278  end if
279  call file_cartesc_read( fid, 'PRES_ref3D', 'ZXY', atmos_refstate_pres(:,:,:) )
280  call file_cartesc_read( fid, 'TEMP_ref3D', 'ZXY', atmos_refstate_temp(:,:,:) )
281  call file_cartesc_read( fid, 'DENS_ref3D', 'ZXY', atmos_refstate_dens(:,:,:) )
282  call file_cartesc_read( fid, 'POTT_ref3D', 'ZXY', atmos_refstate_pott(:,:,:) )
283  call file_cartesc_read( fid, 'QV_ref3D', 'ZXY', atmos_refstate_qv(:,:,:) )
284  call file_cartesc_close( fid )
285 
286  else
287  log_error("ATMOS_REFSTATE_read",*) 'refstate file is not specified.'
288  call prc_abort
289  endif
290 
291  !$acc update device(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
292  !$acc update device(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
293 
294  call atmos_refstate_fillhalo( ka, ks, ke, ia, is, ie, ja, js, je, real_phi(:,:,:) )
295 
296  return
297  end subroutine atmos_refstate_read
298 
299  !-----------------------------------------------------------------------------
301  subroutine atmos_refstate_write
302  use scale_file_cartesc, only: &
303  file_cartesc_write
304  implicit none
305 
306  logical, save :: first = .true.
307  !---------------------------------------------------------------------------
308 
309  if ( .not. first ) return
310  first = .false.
311 
312  if ( atmos_refstate_out_basename /= '' ) then
313 
314  !$acc update host(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
315  !$acc update host(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
316 
317  log_newline
318  log_info("ATMOS_REFSTATE_write",*) 'Output reference state profile '
319 
320  ! 1D
321  call file_cartesc_write( atmos_refstate1d_pres(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
322  'PRES_ref', 'Reference profile of pres.', 'Pa', 'Z', atmos_refstate_out_dtype ) ! [IN]
323  call file_cartesc_write( atmos_refstate1d_temp(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
324  'TEMP_ref', 'Reference profile of temp.', 'K', 'Z', atmos_refstate_out_dtype ) ! [IN]
325  call file_cartesc_write( atmos_refstate1d_dens(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
326  'DENS_ref', 'Reference profile of rho', 'kg/m3', 'Z', atmos_refstate_out_dtype ) ! [IN]
327  call file_cartesc_write( atmos_refstate1d_pott(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
328  'POTT_ref', 'Reference profile of theta', 'K', 'Z', atmos_refstate_out_dtype ) ! [IN]
329  call file_cartesc_write( atmos_refstate1d_qv(:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
330  'QV_ref', 'Reference profile of qv', 'kg/kg', 'Z', atmos_refstate_out_dtype ) ! [IN]
331 
332  ! 3D
333  call file_cartesc_write( atmos_refstate_pres(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
334  'PRES_ref3D', 'Reference profile of pres.', 'Pa', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
335  call file_cartesc_write( atmos_refstate_temp(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
336  'TEMP_ref3D', 'Reference profile of temp.', 'K', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
337  call file_cartesc_write( atmos_refstate_dens(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
338  'DENS_ref3D', 'Reference profile of rho', 'kg/m3', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
339  call file_cartesc_write( atmos_refstate_pott(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
340  'POTT_ref3D', 'Reference profile of theta', 'K', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
341  call file_cartesc_write( atmos_refstate_qv(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, & ! [IN]
342  'QV_ref3D', 'Reference profile of qv', 'kg/kg', 'ZXY', atmos_refstate_out_dtype ) ! [IN]
343 
344  endif
345 
346  return
347  end subroutine atmos_refstate_write
348 
349  !-----------------------------------------------------------------------------
351  subroutine atmos_refstate_generate_isa( &
352  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
353  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
354  use scale_const, only: &
355  epsvap => const_epsvap, &
356  pstd => const_pstd
357  use scale_prc, only: &
358  prc_abort
359  use scale_atmos_profile, only: &
360  profile_isa => atmos_profile_isa
361  use scale_atmos_hydrostatic, only: &
362  hydrostatic_buildrho => atmos_hydrostatic_buildrho
363  use scale_atmos_saturation, only: &
364  saturation_psat_all => atmos_saturation_psat_all, &
365  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
366  implicit none
367  integer, intent(in) :: ka, ks, ke
368  integer, intent(in) :: ia, is, ie
369  integer, intent(in) :: ja, js, je
370 
371  real(rp), intent(in) :: cz ( ka)
372  real(rp), intent(in) :: fz (0:ka)
373  real(rp), intent(in) :: real_cz ( ka,ia,ja)
374  real(rp), intent(in) :: real_fz (0:ka,ia,ja)
375  real(rp), intent(in) :: real_phi( ka,ia,ja)
376 
377  real(rp) :: temp(ka)
378  real(rp) :: pres(ka)
379  real(rp) :: dens(ka)
380  real(rp) :: pott(ka)
381  real(rp) :: qv (ka)
382  real(rp) :: qc (ka)
383 
384  real(rp) :: temp_sfc
385  real(rp) :: pres_sfc
386  real(rp) :: pott_sfc
387  real(rp) :: qv_sfc
388  real(rp) :: qc_sfc
389 
390  real(rp) :: qsat(ka)
391  real(rp) :: psat_sfc
392 
393 #ifdef _OPENACC
394  real(rp) :: work1(ka)
395  real(rp) :: work2(ka)
396  real(rp) :: work3(ka)
397 #endif
398 
399  logical :: converged
400  integer :: k
401  !---------------------------------------------------------------------------
402 
403  pott_sfc = atmos_refstate_temp_sfc
404  pres_sfc = pstd
405 
406  call profile_isa( ka, ks, ke, & ! [IN]
407  pott_sfc, pres_sfc, & ! [IN]
408  cz(:), & ! [IN]
409  pott(:) ) ! [OUT]
410 
411  qv(:) = 0.0_rp
412  qc(:) = 0.0_rp
413  qv_sfc = 0.0_rp
414  qc_sfc = 0.0_rp
415 
416  ! make density & pressure profile in dry condition
417  call hydrostatic_buildrho( ka, ks, ke, &
418  pott(:), qv(:), qc(:), & ! [IN]
419  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
420  cz(:), fz(:), & ! [IN]
421 #ifdef _OPENACC
422  work1(:), work2(:), work3(:), & ! [WORK]
423 #endif
424  dens(:), temp(:), pres(:), & ! [OUT]
425  temp_sfc, & ! [OUT]
426  converged ) ! [OUT]
427 
428  if ( .not. converged ) then
429  log_error("ATMOS_REFSTATE_generate_isa",*) "not converged"
430  call prc_abort
431  end if
432 
433  ! calc QV from RH
434  call saturation_psat_all( temp_sfc, psat_sfc )
435  call saturation_dens2qsat_all( ka, ks, ke, &
436  temp(:), dens(:), & ! [IN]
437  qsat(:) ) ! [OUT]
438 
439  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
440  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp-epsvap) * psat_sfc )
441  do k = ks, ke
442  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
443  enddo
444 
445  ! make density & pressure profile in moist condition
446  call hydrostatic_buildrho( ka, ks, ke, &
447  pott(:), qv(:), qc(:), & ! [IN]
448  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
449  cz(:), fz(:), & ! [IN]
450 #ifdef _OPENACC
451  work1(:), work2(:), work3(:), & ! [WORK]
452 #endif
453  dens(:), temp(:), pres(:), & ! [OUT]
454  temp_sfc, & ! [OUT]
455  converged )
456 
457  if ( .not. converged ) then
458  log_error("ATMOS_REFSTATE_generate_isa",*) "not converged"
459  call prc_abort
460  end if
461 
462  atmos_refstate1d_pres(:) = pres(:)
463  atmos_refstate1d_temp(:) = temp(:)
464  atmos_refstate1d_dens(:) = dens(:)
465  atmos_refstate1d_pott(:) = pott(:)
466  atmos_refstate1d_qv(:) = qv(:)
467 
468  !$acc update device(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
469 
470  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
471  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
472 
473  return
474  end subroutine atmos_refstate_generate_isa
475 
476  !-----------------------------------------------------------------------------
478  subroutine atmos_refstate_generate_uniform( &
479  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
480  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
481  use scale_const, only: &
482  epsvap => const_epsvap, &
483  pstd => const_pstd
484  use scale_prc, only: &
485  prc_abort
486  use scale_atmos_hydrostatic, only: &
487  hydrostatic_buildrho => atmos_hydrostatic_buildrho
488  use scale_atmos_saturation, only: &
489  saturation_psat_all => atmos_saturation_psat_all, &
490  saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
491  implicit none
492  integer, intent(in) :: ka, ks, ke
493  integer, intent(in) :: ia, is, ie
494  integer, intent(in) :: ja, js, je
495 
496  real(rp), intent(in) :: cz ( ka)
497  real(rp), intent(in) :: fz (0:ka)
498  real(rp), intent(in) :: real_cz ( ka,ia,ja)
499  real(rp), intent(in) :: real_fz (0:ka,ia,ja)
500  real(rp), intent(in) :: real_phi( ka,ia,ja)
501 
502  real(rp) :: temp(ka)
503  real(rp) :: pres(ka)
504  real(rp) :: dens(ka)
505  real(rp) :: pott(ka)
506  real(rp) :: qv (ka)
507  real(rp) :: qc (ka)
508 
509  real(rp) :: temp_sfc
510  real(rp) :: pres_sfc
511  real(rp) :: pott_sfc
512  real(rp) :: qv_sfc
513  real(rp) :: qc_sfc
514 
515  real(rp) :: qsat(ka)
516  real(rp) :: psat_sfc
517 
518 #ifdef _OPENACC
519  real(rp) :: work1(ka)
520  real(rp) :: work2(ka)
521  real(rp) :: work3(ka)
522 #endif
523 
524  logical :: converged
525  integer :: k
526  !---------------------------------------------------------------------------
527 
528  pres_sfc = pstd
529  pott_sfc = atmos_refstate_temp_sfc
530  qv_sfc = 0.0_rp
531  qc_sfc = 0.0_rp
532 
533  do k = 1, ka
534  pott(k) = atmos_refstate_pott_uniform
535  qv(k) = 0.0_rp
536  qc(k) = 0.0_rp
537  enddo
538 
539  ! make density & pressure profile in dry condition
540  call hydrostatic_buildrho( ka, ks, ke, &
541  pott(:), qv(:), qc(:), & ! [IN]
542  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
543  cz(:), fz(:), & ! [IN]
544 #ifdef _OPENACC
545  work1(:), work2(:), work3(:), & ! [WORK]
546 #endif
547  dens(:), temp(:), pres(:), & ! [OUT]
548  temp_sfc, & ! [OUT]
549  converged ) ! [OUT]
550  if ( .not. converged ) then
551  log_error("ATMOS_REFSTATE_generate_uniform",*) "not converged"
552  call prc_abort
553  end if
554 
555  ! calc QV from RH
556  call saturation_psat_all( temp_sfc, psat_sfc )
557  call saturation_dens2qsat_all( ka, ks, ke, &
558  temp(:), dens(:), & ! [IN]
559  qsat(:) ) ! [OUT]
560 
561  psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc ! rh * e
562  qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp - epsvap) * psat_sfc )
563  do k = ks, ke
564  qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
565  enddo
566 
567  ! make density & pressure profile in moist condition
568  call hydrostatic_buildrho( ka, ks, ke, &
569  pott(:), qv(:), qc(:), & ! [IN]
570  pres_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
571  cz(:), fz(:), & ! [IN]
572 #ifdef _OPENACC
573  work1(:), work2(:), work3(:), & ! [WORK]
574 #endif
575  dens(:), temp(:), pres(:), & ! [OUT]
576  temp_sfc, & ! [OUT]
577  converged ) ! [OUT]
578  if ( .not. converged ) then
579  log_error("ATMOS_REFSTATE_generate_uniform",*) "not converged"
580  call prc_abort
581  end if
582 
583  do k = ks, ke
584  atmos_refstate1d_pres(k) = pres(k)
585  atmos_refstate1d_temp(k) = temp(k)
586  atmos_refstate1d_dens(k) = dens(k)
587  atmos_refstate1d_pott(k) = pott(k)
588  atmos_refstate1d_qv(k) = qv(k)
589  end do
590 
591  !$acc update device(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
592 
593  call atmos_refstate_calc3d( ka, ks, ke, ia, is, ie, ja, js, je, &
594  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
595 
596  return
597  end subroutine atmos_refstate_generate_uniform
598 
599  !-----------------------------------------------------------------------------
601  subroutine atmos_refstate_generate_zero( KA, IA, JA )
602  implicit none
603  integer, intent(in) :: ka, ia, ja
604 
605  integer :: k, i, j
606  !---------------------------------------------------------------------------
607 
608  !$omp parallel do
609  do k = 1, ka
610  atmos_refstate1d_pres(k) = 0.0_rp
611  atmos_refstate1d_temp(k) = 0.0_rp
612  atmos_refstate1d_dens(k) = 0.0_rp
613  atmos_refstate1d_pott(k) = 0.0_rp
614  atmos_refstate1d_qv(k) = 0.0_rp
615  end do
616 
617  !$acc update device(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
618 
619  !$omp parallel do collapse(3)
620  do j = 1, ja
621  do i = 1, ia
622  do k = 1, ka
623  atmos_refstate_dens(k,i,j) = 0.0_rp
624  atmos_refstate_temp(k,i,j) = 0.0_rp
625  atmos_refstate_pres(k,i,j) = 0.0_rp
626  atmos_refstate_pott(k,i,j) = 0.0_rp
627  atmos_refstate_qv(k,i,j) = 0.0_rp
628  enddo
629  enddo
630  enddo
631 
632  !$acc update device(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
633 
634  return
635  end subroutine atmos_refstate_generate_zero
636 
637  !-----------------------------------------------------------------------------
639  subroutine atmos_refstate_update( &
640  KA, KS, KE, IA, IS, IE, ISB, IEB, JA, JS, JE, JSB, JEB, &
641  DENS, POTT, TEMP, PRES, QV, &
642  CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA, &
643  nowsec )
644  use scale_statistics, only: &
645  statistics_horizontal_mean
646  use scale_interp_vert, only: &
647  interp_vert_xi2z
648  implicit none
649 
650  integer, intent(in) :: ka, ks, ke
651  integer, intent(in) :: ia, is, ie, isb, ieb
652  integer, intent(in) :: ja, js, je, jsb, jeb
653  real(rp), intent(in) :: dens (ka,ia,ja)
654  real(rp), intent(in) :: pott (ka,ia,ja)
655  real(rp), intent(in) :: temp (ka,ia,ja)
656  real(rp), intent(in) :: pres (ka,ia,ja)
657  real(rp), intent(in) :: qv (ka,ia,ja)
658  real(rp), intent(in) :: cz ( ka)
659  real(rp), intent(in) :: fz (0:ka)
660  real(rp), intent(in) :: fdz ( ka-1)
661  real(rp), intent(in) :: rcdz ( ka)
662  real(rp), intent(in) :: real_cz ( ka,ia,ja)
663  real(rp), intent(in) :: real_fz (0:ka,ia,ja)
664  real(rp), intent(in) :: real_phi( ka,ia,ja)
665  real(rp), intent(in) :: area ( ia,ja)
666  real(dp), intent(in) :: nowsec
667 
668  real(rp) :: work(ka,ia,ja)
669 
670  integer :: k
671  !---------------------------------------------------------------------------
672 
673  if ( .not. atmos_refstate_update_flag ) return
674 
675  if ( last_updated < 0.0_rp .or. ( nowsec - last_updated >= atmos_refstate_update_dt ) ) then
676 
677  !$acc data copyin(DENS, POTT, TEMP, PRES, QV, CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA) create(work)
678 
679  log_info("ATMOS_REFSTATE_update",*) 'update reference state'
680 
681  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
682  cz(:), real_cz(:,:,:), temp(:,:,:), & ! [IN]
683  work(:,:,:) ) ! [OUT]
684  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
685  work(:,:,:), area(:,:), & ! [IN]
686  atmos_refstate1d_temp(:) ) ! [OUT]
687 
688  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
689  cz(:), real_cz(:,:,:), pres(:,:,:), & ! [IN]
690  work(:,:,:) ) ! [OUT]
691  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
692  work(:,:,:), area(:,:), & ! [IN]
693  atmos_refstate1d_pres(:) ) ! [OUT]
694 
695  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
696  cz(:), real_cz(:,:,:), dens(:,:,:), & ! [IN]
697  work(:,:,:) ) ! [OUT]
698  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
699  work(:,:,:), area(:,:), & ! [IN]
700  atmos_refstate1d_dens(:) ) ! [OUT]
701 
702  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
703  cz(:), real_cz(:,:,:), pott(:,:,:), & ! [IN]
704  work(:,:,:) ) ! [OUT]
705  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
706  work(:,:,:), area(:,:), & ! [IN]
707  atmos_refstate1d_pott(:) ) ! [OUT]
708 
709  call interp_vert_xi2z( ka, ks, ke, ia, is, ie, ja, js, je, &
710  cz(:), real_cz(:,:,:), qv(:,:,:), & ! [IN]
711  work(:,:,:) ) ! [OUT]
712  call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
713  work(:,:,:), area(:,:), & ! [IN]
714  atmos_refstate1d_qv(:) ) ! [OUT]
715 
716  !$acc update host(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
717 
718  do k = ke-1, ks, -1 ! fill undefined value
719  if( atmos_refstate1d_dens(k) <= 0.0_rp ) atmos_refstate1d_dens(k) = atmos_refstate1d_dens(k+1)
720  if( atmos_refstate1d_temp(k) <= 0.0_rp ) atmos_refstate1d_temp(k) = atmos_refstate1d_temp(k+1)
721  if( atmos_refstate1d_pres(k) <= 0.0_rp ) atmos_refstate1d_pres(k) = atmos_refstate1d_pres(k+1)
722  if( atmos_refstate1d_pott(k) <= 0.0_rp ) atmos_refstate1d_pott(k) = atmos_refstate1d_pott(k+1)
723  if( atmos_refstate1d_qv(k) <= 0.0_rp ) atmos_refstate1d_qv(k) = atmos_refstate1d_qv(k+1)
724  enddo
725 
726  call atmos_refstate_smoothing( ka, ks, ke, fdz(:), rcdz(:), atmos_refstate1d_pott(:) )
727  call atmos_refstate_smoothing( ka, ks, ke, fdz(:), rcdz(:), atmos_refstate1d_qv(:) )
728 
729  !$acc update device(ATMOS_REFSTATE1D_pres, ATMOS_REFSTATE1D_temp, ATMOS_REFSTATE1D_dens, ATMOS_REFSTATE1D_pott, ATMOS_REFSTATE1D_qv)
730 
731  call atmos_refstate_calc3d( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
732  cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
733 
734  last_updated = nowsec
735 
736 
737  log_newline
738  log_info("ATMOS_REFSTATE_update",*) 'Generated reference state of atmosphere:'
739  log_info_cont(*) '============================================================================='
740  log_info_cont(*) ' z*-coord.: pressure: temperature: density: pot.temp.: water vapor'
741  do k = ks, ke
742  log_info_cont('(6F13.5)') cz(k), &
743  atmos_refstate1d_pres(k), &
744  atmos_refstate1d_temp(k), &
745  atmos_refstate1d_dens(k), &
746  atmos_refstate1d_pott(k), &
747  atmos_refstate1d_qv(k)
748  enddo
749  log_info_cont(*) '============================================================================='
750 
751  ! output reference profile
752  call atmos_refstate_write
753 
754  !$acc end data
755 
756  endif
757 
758  if ( atmos_refstate_update_dt < 0.0_rp ) atmos_refstate_update_flag = .false.
759 
760  return
761  end subroutine atmos_refstate_update
762 
763  !-----------------------------------------------------------------------------
765  subroutine atmos_refstate_calc3d( &
766  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
767  CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
768  use scale_prc, only: &
769  prc_abort
770  use scale_interp_vert, only: &
771  interp_vert_z2xi
772  use scale_atmos_hydrostatic, only: &
773  hydrostatic_buildrho_atmos => atmos_hydrostatic_buildrho_atmos, &
774  hydrostatic_buildrho_atmos_rev => atmos_hydrostatic_buildrho_atmos_rev
775  implicit none
776  integer, intent(in) :: KA, KS, KE
777  integer, intent(in) :: IA, IS, IE
778  integer, intent(in) :: JA, JS, JE
779 
780  real(RP), intent(in) :: CZ ( KA)
781  real(RP), intent(in) :: FZ (0:KA)
782  real(RP), intent(in) :: REAL_CZ ( KA,IA,JA)
783  real(RP), intent(in) :: REAL_FZ (0:KA,IA,JA)
784  real(RP), intent(in) :: REAL_PHI( KA,IA,JA)
785 
786  real(RP) :: dens(KA,IA,JA)
787  real(RP) :: temp(KA,IA,JA)
788  real(RP) :: pres(KA,IA,JA)
789  real(RP) :: pott(KA,IA,JA)
790  real(RP) :: qv (KA,IA,JA)
791  real(RP) :: qc (KA,IA,JA)
792  real(RP) :: dz (KA,IA,JA)
793 
794  real(RP) :: pott1(IA,JA)
795  real(RP) :: pott2(IA,JA)
796  real(RP) :: dens1(IA,JA)
797  real(RP) :: dens2(IA,JA)
798  real(RP) :: temp1(IA,JA)
799  real(RP) :: pres1(IA,JA)
800  real(RP) :: qv1 (IA,JA)
801  real(RP) :: qv2 (IA,JA)
802  real(RP) :: qc1 (IA,JA)
803  real(RP) :: qc2 (IA,JA)
804  real(RP) :: dz1 (IA,JA)
805 
806  real(RP) :: dens_toa_1D
807  real(RP) :: temp_toa_1D
808  real(RP) :: pres_toa_1D
809  real(RP) :: qc_1D
810  real(RP) :: dz_1D
811 
812  real(RP) :: work(KA,IA,JA)
813 
814  logical :: converged
815  integer :: k, i, j
816  !---------------------------------------------------------------------------
817 
818  !$acc data copyin(CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI) &
819  !$acc create(dens, temp, pres, pott, qv, qc, dz, work, pott1, pott2, dens1, dens2, temp1, pres1, qv1, qv2, qc1, qc2, dz1)
820 
821  !--- potential temperature
822  !$acc kernels
823  do j = js, je
824  do i = is, ie
825  work(:,i,j) = atmos_refstate1d_pott(:)
826  enddo
827  enddo
828  !$acc end kernels
829 
830  call interp_vert_z2xi( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
831  real_cz(:,:,:), cz(:), & ! [IN]
832  work(:,:,:), & ! [IN]
833  pott(:,:,:) ) ! [OUT]
834 
835  !--- water vapor
836  !$acc kernels
837  do j = js, je
838  do i = is, ie
839  work(:,i,j) = atmos_refstate1d_qv(:)
840  enddo
841  enddo
842  !$acc end kernels
843 
844  call interp_vert_z2xi( ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
845  real_cz(:,:,:), cz(:), & ! [IN]
846  work(:,:,:), & ! [IN]
847  qv(:,:,:) ) ! [OUT]
848 
849  !--- build up density to TOA (1D)
850  qc_1d = 0.0_rp
851  dz_1d = fz(ke) - cz(ke)
852 
853  call hydrostatic_buildrho_atmos( atmos_refstate1d_pott(ke), & ! [IN]
854  atmos_refstate1d_qv(ke), & ! [IN]
855  qc_1d, & ! [IN]
856  atmos_refstate1d_dens(ke), & ! [IN]
857  atmos_refstate1d_pott(ke), & ! [IN]
858  atmos_refstate1d_qv(ke), & ! [IN]
859  qc_1d, & ! [IN]
860  dz_1d, & ! [IN]
861  ke+1, & ! [IN]
862  dens_toa_1d, & ! [OUT]
863  temp_toa_1d, & ! [OUT]
864  pres_toa_1d, & ! [OUT]
865  converged ) ! [OUT]
866 
867  if ( .not. converged ) then
868  log_error("ATMOS_REFSTATE_calc3D",*) "not converged"
869  call prc_abort
870  end if
871 
872  ! build down density from TOA (3D)
873  !$acc kernels
874  do j = js, je
875  do i = is, ie
876  dz(ks-1,i,j) = real_cz(ks,i,j) - real_fz(ks-1,i,j) ! distance from surface to cell center
877  do k = ks, ke-1
878  dz(k,i,j) = real_cz(k+1,i,j) - real_cz(k,i,j) ! distance from cell center to cell center
879  enddo
880  dz(ke,i,j) = real_fz(ke,i,j) - real_cz(ke,i,j) ! distance from cell center to TOA
881  enddo
882  enddo
883  !$acc end kernels
884 
885  !$acc kernels
886  do j = js, je
887  do i = is, ie
888  dens(ke+1,i,j) = dens_toa_1d
889  temp(ke+1,i,j) = temp_toa_1d
890  pres(ke+1,i,j) = pres_toa_1d
891  pott(ke+1,i,j) = pott(ke,i,j)
892  qv(ke+1,i,j) = qv(ke,i,j)
893  enddo
894  enddo
895  !$acc end kernels
896 
897  !$acc kernels
898  do j = js, je
899  do i = is, ie
900  pott(ks-1,i,j) = pott(ks,i,j)
901  qv(ks-1,i,j) = qv(ks,i,j)
902  enddo
903  enddo
904  !$acc end kernels
905 
906  !$acc kernels
907  qc(:,:,:) = 0.0_rp
908  !$acc end kernels
909 
910 
911  !$acc kernels
912  do j = js, je
913  do i = is, ie
914  pott1(i,j) = pott(ke,i,j)
915  qv1(i,j) = qv(ke,i,j)
916  qc1(i,j) = qc(ke,i,j)
917  dens2(i,j) = dens(ke+1,i,j)
918  pott2(i,j) = pott(ke+1,i,j)
919  qv2(i,j) = qv(ke+1,i,j)
920  qc2(i,j) = qc(ke+1,i,j)
921  dz1(i,j) = dz(ke,i,j)
922  end do
923  end do
924  !$acc end kernels
925  call hydrostatic_buildrho_atmos_rev( ia, is, ie, ja, js, je, &
926  pott1(:,:), qv1(:,:), qc1(:,:), & ! [IN]
927  dens2(:,:), pott2(:,:), qv2(:,:), qc2(:,:), & ! [IN]
928  dz1(:,:), ke+1, & ! [IN]
929  dens1(:,:), temp1(:,:), pres1(:,:) ) ! [OUT]
930  !$acc kernels
931  do j = js, je
932  do i = is, ie
933  dens(ke,i,j) = dens1(i,j)
934  temp(ke,i,j) = temp1(i,j)
935  pres(ke,i,j) = pres1(i,j)
936  end do
937  end do
938  !$acc end kernels
939 
940  call hydrostatic_buildrho_atmos_rev( ka, ks, ke, ia, is, ie, ja, js, je, &
941  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
942  dz(:,:,:), & ! [IN]
943  dens(:,:,:), & ! [INOUT]
944  temp(:,:,:), pres(:,:,:) ) ! [OUT]
945 
946  !$acc kernels
947  do j = js, je
948  do i = is, ie
949  do k = ks, ke
950  atmos_refstate_dens(k,i,j) = dens(k,i,j)
951  atmos_refstate_temp(k,i,j) = temp(k,i,j)
952  atmos_refstate_pres(k,i,j) = pres(k,i,j)
953  atmos_refstate_pott(k,i,j) = pott(k,i,j)
954  atmos_refstate_qv(k,i,j) = qv(k,i,j)
955  enddo
956  enddo
957  enddo
958  !$acc end kernels
959 
960  !$acc end data
961 
962  call atmos_refstate_fillhalo( ka, ks, ke, ia, is, ie, ja, js, je, real_phi(:,:,:) )
963 
964  return
965  end subroutine atmos_refstate_calc3d
966 
967  !-----------------------------------------------------------------------------
968  subroutine atmos_refstate_fillhalo( &
969  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
970  REAL_PHI )
971  use scale_const, only: &
972  rdry => const_rdry, &
973  cpdry => const_cpdry, &
974  p00 => const_pre00, &
975  undef => const_undef
976  use scale_comm_cartesc, only: &
977  comm_vars8, &
978  comm_wait
979  integer, intent(in) :: KA, KS, KE
980  integer, intent(in) :: IA, IS, IE
981  integer, intent(in) :: JA, JS, JE
982  real(RP), intent(in) :: REAL_PHI(KA,IA,JA)
983 
984  real(RP) :: RovCP
985  integer :: k, i, j
986 
987  rovcp = rdry / cpdry
988 
989  !$acc kernels
990  !$acc loop collapse(2)
991  do j = js, je
992  do i = is, ie
993  !$acc loop seq
994  do k = 1, ks-1
995  atmos_refstate_temp(k, i,j) = atmos_refstate_temp(ks,i,j)
996  end do
997  !$acc loop seq
998  do k = ke+1, ka
999  atmos_refstate_temp(k,i,j) = atmos_refstate_temp(ke,i,j)
1000  end do
1001 
1002  !$acc loop seq
1003  do k = 1, ks-1
1004  atmos_refstate_qv(k, i,j) = atmos_refstate_qv(ks,i,j)
1005  end do
1006  !$acc loop seq
1007  do k = ke+1, ka
1008  atmos_refstate_qv(k,i,j) = atmos_refstate_qv(ke,i,j)
1009  end do
1010 
1011  !$acc loop seq
1012  do k = 1, ks-2
1013  atmos_refstate_pres(k, i,j) = undef
1014  end do
1015  atmos_refstate_pres(ks-1, i,j) = atmos_refstate_pres(ks+1,i,j) &
1016  - atmos_refstate_dens(ks ,i,j) * ( real_phi(ks-1,i,j) - real_phi(ks+1,i,j) )
1017  atmos_refstate_pres(ke+1, i,j) = max( 1e-20_rp, &
1018  atmos_refstate_pres(ke-1,i,j) &
1019  - atmos_refstate_dens(ke ,i,j) * ( real_phi(ke+1,i,j) - real_phi(ke-1,i,j) ) )
1020  !$acc loop seq
1021  do k = ke+2, ka
1022  atmos_refstate_pres(k,i,j) = undef
1023  end do
1024 
1025  !$acc loop seq
1026  do k = 1, ks-2
1027  atmos_refstate_dens(k, i,j) = undef
1028  end do
1029  atmos_refstate_dens(ks-1, i,j) = atmos_refstate_pres(ks-1,i,j) / ( atmos_refstate_temp(ks-1,i,j) * rdry )
1030  atmos_refstate_dens(ke+1, i,j) = atmos_refstate_pres(ke+1,i,j) / ( atmos_refstate_temp(ke+1,i,j) * rdry )
1031  !$acc loop seq
1032  do k = ke+2, ka
1033  atmos_refstate_dens(k,i,j) = undef
1034  end do
1035 
1036  !$acc loop seq
1037  do k = 1, ks-2
1038  atmos_refstate_pott(k, i,j) = undef
1039  end do
1040  atmos_refstate_pott(ks-1, i,j) = atmos_refstate_temp(ks-1,i,j) * ( p00 / atmos_refstate_pres(ks-1,i,j) )**rovcp
1041  atmos_refstate_pott(ke+1, i,j) = atmos_refstate_temp(ke+1,i,j) * ( p00 / atmos_refstate_pres(ke+1,i,j) )**rovcp
1042  !$acc loop seq
1043  do k = ke+2, ka
1044  atmos_refstate_pott(k,i,j) = undef
1045  end do
1046  enddo
1047  enddo
1048  !$acc end kernels
1049 
1050  call comm_vars8( atmos_refstate_dens(:,:,:), 1 )
1051  call comm_vars8( atmos_refstate_temp(:,:,:), 2 )
1052  call comm_vars8( atmos_refstate_pres(:,:,:), 3 )
1053  call comm_vars8( atmos_refstate_pott(:,:,:), 4 )
1054  call comm_vars8( atmos_refstate_qv(:,:,:), 5 )
1055  call comm_wait ( atmos_refstate_dens(:,:,:), 1, .false. )
1056  call comm_wait ( atmos_refstate_temp(:,:,:), 2, .false. )
1057  call comm_wait ( atmos_refstate_pres(:,:,:), 3, .false. )
1058  call comm_wait ( atmos_refstate_pott(:,:,:), 4, .false. )
1059  call comm_wait ( atmos_refstate_qv(:,:,:), 5, .false. )
1060 
1061  !$acc update host(ATMOS_REFSTATE_pres, ATMOS_REFSTATE_temp, ATMOS_REFSTATE_dens, ATMOS_REFSTATE_pott, ATMOS_REFSTATE_qv)
1062 
1063  return
1064  end subroutine atmos_refstate_fillhalo
1065 
1066  !-----------------------------------------------------------------------------
1067 !OCL SERIAL
1068  subroutine atmos_refstate_smoothing( &
1069  KA, KS, KE, &
1070  FDZ, RCDZ, &
1071  phi )
1072  use scale_const, only: &
1073  eps => const_eps
1074  implicit none
1075  integer, intent(in) :: KA, KS, KE
1076 
1077  real(RP), intent(in) :: FDZ (KA-1)
1078  real(RP), intent(in) :: RCDZ(KA)
1079 
1080  real(RP), intent(inout) :: phi(KA)
1081 
1082  real(RP) :: dev (KA)
1083  real(RP) :: correction(KA)
1084  real(RP) :: fact(KA)
1085 
1086  integer, parameter :: iter_max = 100
1087  real(RP) :: sig0, sig1, zerosw
1088  logical :: updated
1089  integer :: k, iter
1090  !---------------------------------------------------------------------------
1091 
1092  dev(ks) = 0.0_rp
1093  dev(ke) = 0.0_rp
1094 
1095  correction(ks-1:ks+1) = 0.0_rp
1096  correction(ke-1:ke+1) = 0.0_rp
1097 
1098  fact(ks-1:ks+1) = 0.0_rp
1099  fact(ke-1:ke+1) = 0.0_rp
1100 
1101  do iter = 1, iter_max
1102  updated = .false.
1103 
1104  do k = ks+1, ke-1
1105  dev(k) = phi(k) - ( fdz(k-1)*phi(k+1) + fdz(k)*phi(k-1) ) / ( fdz(k) + fdz(k-1) )
1106  enddo
1107 
1108  do k = ks+2, ke-2
1109  sig0 = dev(k) * dev(k-1)
1110  sig1 = dev(k) * dev(k+1)
1111  if ( sig0 < -eps .and. sig1 < -eps ) then
1112  correction(k) = dev(k) &
1113  / ( 2.0_rp*rcdz(k) + ( fdz(k-1)*rcdz(k+1) + fdz(k)*rcdz(k-1) ) / ( fdz(k) + fdz(k-1) ) )
1114  updated = .true.
1115  else
1116  correction(k) = 0.0_rp
1117  end if
1118  enddo
1119 
1120  sig1 = dev(ks+1) * dev(ks+2)
1121  if ( sig1 < -eps ) then
1122  correction(ks+1) = dev(ks+1) &
1123  / ( 2.0_rp*rcdz(ks+1) + (fdz(ks)*rcdz(ks+2)+fdz(ks+1)*rcdz(ks))/(fdz(ks+1)+fdz(ks)) )
1124  updated = .true.
1125  else
1126  correction(ks+1) = 0.0_rp
1127  end if
1128 
1129  sig0 = dev(ke-1) * dev(ke-2)
1130  if ( sig0 < -eps ) then
1131  correction(ke-1) = dev(ke-1) &
1132  / ( 2.0_rp*rcdz(ke-1) + (fdz(ke-2)*rcdz(ke)+fdz(ke-1)*rcdz(ke-2))/(fdz(ke-1)+fdz(ke-2)) )
1133  updated = .true.
1134  else
1135  correction(ke-1) = 0.0_rp
1136  end if
1137 
1138  if ( .NOT. updated ) exit
1139 
1140  do k = ks+1, ke-1
1141  zerosw = 0.5_rp - sign( 0.5_rp, abs(correction(k))-eps ) ! if correction(k) == 0 then fact(k) = 0.0
1142  fact(k) = correction(k) / ( correction(k) - correction(k+1) - correction(k-1) + zerosw )
1143  enddo
1144 
1145  do k = ks, ke
1146  phi(k) = phi(k) + ( correction(k+1) * fact(k+1) &
1147  - correction(k ) * fact(k ) * 2.0_rp &
1148  + correction(k-1) * fact(k-1) ) * rcdz(k)
1149  enddo
1150 
1151  if ( iter == iter_max ) then
1152  log_info("ATMOS_REFSTATE_smoothing",*) "iteration not converged!", phi
1153  endif
1154  enddo
1155 
1156  return
1157  end subroutine atmos_refstate_smoothing
1158 
1159 end module scale_atmos_refstate
scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:75
scale_atmos_refstate::atmos_refstate_dens
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_dens
refernce density [kg/m3]
Definition: scale_atmos_refstate.F90:39
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_profile
module atmosphere / vertical profile
Definition: scale_atmos_profile.F90:12
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_atmos_refstate::atmos_refstate_pres
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pres
refernce pressure [Pa]
Definition: scale_atmos_refstate.F90:37
scale_atmos_refstate
module atmosphere / reference state
Definition: scale_atmos_refstate.F90:12
scale_atmos_refstate::atmos_refstate_pott
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_pott
refernce potential temperature [K]
Definition: scale_atmos_refstate.F90:40
scale_atmos_refstate::atmos_refstate_fillhalo
subroutine atmos_refstate_fillhalo(KA, KS, KE, IA, IS, IE, JA, JS, JE, REAL_PHI)
Definition: scale_atmos_refstate.F90:971
scale_atmos_refstate::atmos_refstate_qv
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_qv
refernce vapor [kg/kg]
Definition: scale_atmos_refstate.F90:41
scale_atmos_refstate::atmos_refstate_update
subroutine, public atmos_refstate_update(KA, KS, KE, IA, IS, IE, ISB, IEB, JA, JS, JE, JSB, JEB, DENS, POTT, TEMP, PRES, QV, CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA, nowsec)
Update reference state profile (Horizontal average)
Definition: scale_atmos_refstate.F90:644
scale_interp_vert::interp_vert_z2xi
subroutine, public interp_vert_z2xi(KA, KS, KE, IA, IS, IE, JA, JS, JE, Z, Xi, var, var_Xi)
Definition: scale_interp_vert.F90:265
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const::const_pstd
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:96
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_refstate::atmos_refstate_temp
real(rp), dimension(:,:,:), allocatable, public atmos_refstate_temp
refernce temperature [K]
Definition: scale_atmos_refstate.F90:38
scale_file_cartesc::file_cartesc_close
subroutine, public file_cartesc_close(fid)
Close a netCDF file.
Definition: scale_file_cartesC.F90:1044
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:60
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_refstate::atmos_refstate_write
subroutine, public atmos_refstate_write
Write reference state profile.
Definition: scale_atmos_refstate.F90:302
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_interp_vert::interp_vert_xi2z
subroutine, public interp_vert_xi2z(KA, KS, KE, IA, IS, IE, JA, JS, JE, Xi, Z, var, var_Z)
Definition: scale_interp_vert.F90:212
scale_atmos_refstate::atmos_refstate_read
subroutine, public atmos_refstate_read(KA, KS, KE, IA, IS, IE, JA, JS, JE, REAL_PHI)
Read reference state profile.
Definition: scale_atmos_refstate.F90:244
scale_atmos_hydrostatic
module atmosphere / hydrostatic barance
Definition: scale_atmos_hydrostatic.F90:12
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_atmos_refstate::atmos_refstate_finalize
subroutine, public atmos_refstate_finalize
Definition: scale_atmos_refstate.F90:220
scale_file_cartesc::file_cartesc_open
subroutine, public file_cartesc_open(basename, fid, single, aggregate)
open a netCDF file for read
Definition: scale_file_cartesC.F90:760
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:59
scale_atmos_refstate::atmos_refstate_smoothing
subroutine atmos_refstate_smoothing(KA, KS, KE, FDZ, RCDZ, phi)
Definition: scale_atmos_refstate.F90:1072
scale_atmos_refstate::atmos_refstate_setup
subroutine, public atmos_refstate_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI)
Setup.
Definition: scale_atmos_refstate.F90:83
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_interp_vert
module INTERPOLATION vertical
Definition: scale_interp_vert.F90:11
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:97
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_file_cartesc
module file / cartesianC
Definition: scale_file_cartesC.F90:11
scale_atmos_refstate::atmos_refstate_calc3d
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
Definition: scale_atmos_refstate.F90:768
Generated by   doxygen 1.8.17