SCALE-RM
scale_atmos_hydrostatic.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_hydrostatic
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_const, only: &
21  grav => const_grav, &
22  rdry => const_rdry, &
23  rvap => const_rvap, &
24  cvdry => const_cvdry, &
25  cvvap => const_cvvap, &
26  cpdry => const_cpdry, &
27  cpvap => const_cpvap, &
28  laps => const_laps, &
29  lapsdry => const_lapsdry, &
30  p00 => const_pre00, &
31  epstvap => const_epstvap
32  !-----------------------------------------------------------------------------
33  implicit none
34  private
35  !-----------------------------------------------------------------------------
36  !
37  !++ Public procedure
38  !
39  public :: atmos_hydrostatic_setup
40  public :: atmos_hydrostatic_buildrho
41  public :: atmos_hydrostatic_buildrho_real
42  public :: atmos_hydrostatic_buildrho_atmos
43  public :: atmos_hydrostatic_buildrho_atmos_rev
44  public :: atmos_hydrostatic_buildrho_bytemp
45  public :: atmos_hydrostatic_buildrho_bytemp_atmos
46 
47  public :: atmos_hydrostatic_barometric_law_mslp
48  public :: atmos_hydrostatic_barometric_law_pres
49 
50  interface atmos_hydrostatic_buildrho
51  module procedure atmos_hydrostatic_buildrho_1d
52  module procedure atmos_hydrostatic_buildrho_3d
53  end interface atmos_hydrostatic_buildrho
54 
55  interface atmos_hydrostatic_buildrho_real
56  module procedure atmos_hydrostatic_buildrho_real_3d
57  end interface atmos_hydrostatic_buildrho_real
58 
59  interface atmos_hydrostatic_buildrho_atmos
60  module procedure atmos_hydrostatic_buildrho_atmos_0d
61  module procedure atmos_hydrostatic_buildrho_atmos_1d
62  module procedure atmos_hydrostatic_buildrho_atmos_3d
63  end interface atmos_hydrostatic_buildrho_atmos
64 
65  interface atmos_hydrostatic_buildrho_atmos_rev
66  module procedure atmos_hydrostatic_buildrho_atmos_rev_1d
67  module procedure atmos_hydrostatic_buildrho_atmos_rev_2d
68  module procedure atmos_hydrostatic_buildrho_atmos_rev_3d
69  end interface atmos_hydrostatic_buildrho_atmos_rev
70 
71  interface atmos_hydrostatic_buildrho_bytemp
72  module procedure atmos_hydrostatic_buildrho_bytemp_1d
73  module procedure atmos_hydrostatic_buildrho_bytemp_3d
74  end interface atmos_hydrostatic_buildrho_bytemp
75 
76  interface atmos_hydrostatic_buildrho_bytemp_atmos
77  module procedure atmos_hydrostatic_buildrho_bytemp_atmos_1d
78  module procedure atmos_hydrostatic_buildrho_bytemp_atmos_3d
79  end interface atmos_hydrostatic_buildrho_bytemp_atmos
80 
81  interface atmos_hydrostatic_barometric_law_mslp
82  module procedure atmos_hydrostatic_barometric_law_mslp_0d
83  module procedure atmos_hydrostatic_barometric_law_mslp_2d
84  end interface atmos_hydrostatic_barometric_law_mslp
85 
86  interface atmos_hydrostatic_barometric_law_pres
87  module procedure atmos_hydrostatic_barometric_law_pres_0d
88  module procedure atmos_hydrostatic_barometric_law_pres_2d
89  end interface atmos_hydrostatic_barometric_law_pres
90 
91  !-----------------------------------------------------------------------------
92  !
93  !++ Public parameters & variables
94  !
95  !-----------------------------------------------------------------------------
96  !
97  !++ Private procedure
98  !
99  private :: atmos_hydrostatic_buildrho_atmos_2d
100 
101  !-----------------------------------------------------------------------------
102  !
103  !++ Private parameters & variables
104  !
105  integer, private, parameter :: itelim = 100
106  real(RP), private :: criteria
107  logical, private :: hydrostatic_uselapserate = .false.
108  integer, private :: hydrostatic_buildrho_real_kref = 1
109  integer, private :: hydrostatic_barometric_law_mslp_kref = 1
110 
111  !-----------------------------------------------------------------------------
112 contains
113  !-----------------------------------------------------------------------------
115  subroutine atmos_hydrostatic_setup
116  use scale_prc, only: &
117  prc_abort
118  use scale_const, only: &
119  const_eps
120  implicit none
121 
122  namelist / param_atmos_hydrostatic / &
123  hydrostatic_uselapserate, &
124  hydrostatic_buildrho_real_kref, &
125  hydrostatic_barometric_law_mslp_kref
126 
127  integer :: ierr
128  !---------------------------------------------------------------------------
129 
130  log_newline
131  log_info("ATMOS_HYDROSTATIC_setup",*) 'Setup'
132 
133  !--- read namelist
134  rewind(io_fid_conf)
135  read(io_fid_conf,nml=param_atmos_hydrostatic,iostat=ierr)
136  if( ierr < 0 ) then !--- missing
137  log_info("ATMOS_HYDROSTATIC_setup",*) 'Not found namelist. Default used.'
138  elseif( ierr > 0 ) then !--- fatal error
139  log_error("ATMOS_HYDROSTATIC_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_HYDROSTATIC. Check!'
140  call prc_abort
141  endif
142  log_nml(param_atmos_hydrostatic)
143 
144  criteria = sqrt( const_eps )
145 
146  log_newline
147  log_info("ATMOS_HYDROSTATIC_setup",*) 'Use lapse rate for estimation of surface temperature? : ', hydrostatic_uselapserate
148  log_info("ATMOS_HYDROSTATIC_setup",*) 'Buildrho conversion criteria : ', criteria
149 
150  return
151  end subroutine atmos_hydrostatic_setup
152 
153  !-----------------------------------------------------------------------------
155 !OCL SERIAL
156  subroutine atmos_hydrostatic_buildrho_1d( &
157  KA, KS, KE, &
158  pott, qv, qc, &
159  pres_sfc, pott_sfc, qv_sfc, qc_sfc, &
160  cz, fz, &
161  dens, temp, pres, &
162  temp_sfc )
163  use scale_prc, only: &
164  prc_abort
165  use scale_atmos_hydrometeor, only: &
166  cv_vapor, &
167  cv_water, &
168  cp_vapor, &
169  cp_water
170  implicit none
171 
172  integer, intent(in) :: KA, KS, KE
173  real(RP), intent(in) :: pott(ka)
174  real(RP), intent(in) :: qv (ka)
175  real(RP), intent(in) :: qc (ka)
176  real(RP), intent(in) :: pres_sfc
177  real(RP), intent(in) :: pott_sfc
178  real(RP), intent(in) :: qv_sfc
179  real(RP), intent(in) :: qc_sfc
180  real(RP), intent(in) :: cz (ka)
181  real(RP), intent(in) :: fz (0:ka)
182  real(RP), intent(out) :: dens(ka)
183  real(RP), intent(out) :: temp(ka)
184  real(RP), intent(out) :: pres(ka)
185  real(RP), intent(out) :: temp_sfc
186 
187  real(RP) :: dens_sfc
188  real(RP) :: Rtot_sfc
189  real(RP) :: CVtot_sfc
190  real(RP) :: CPtot_sfc
191  real(RP) :: CPovCV_sfc
192  real(RP) :: Rtot
193  real(RP) :: CVtot
194  real(RP) :: CPtot
195  real(RP) :: CPovCV
196  real(RP) :: dz(ka)
197 
198  integer :: k
199  !---------------------------------------------------------------------------
200 
201  !--- from surface to lowermost atmosphere
202 
203  rtot_sfc = rdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
204  + rvap * qv_sfc
205  cvtot_sfc = cvdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
206  + cv_vapor * qv_sfc &
207  + cv_water * qc_sfc
208  cptot_sfc = cpdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
209  + cp_vapor * qv_sfc &
210  + cp_water * qc_sfc
211  cpovcv_sfc = cptot_sfc / cvtot_sfc
212 
213  rtot = rdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
214  + rvap * qv(ks)
215  cvtot = cvdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
216  + cv_vapor * qv(ks) &
217  + cv_water * qc(ks)
218  cptot = cpdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
219  + cp_vapor * qv(ks) &
220  + cp_water * qc(ks)
221  cpovcv = cptot / cvtot
222 
223  ! density at surface
224  dens_sfc = p00 / rtot_sfc / pott_sfc * ( pres_sfc/p00 )**(cvtot_sfc/cptot_sfc)
225  temp_sfc = pres_sfc / ( dens_sfc * rtot_sfc )
226 
227  dz(ks-1) = cz(ks) - fz(ks-1)
228  do k = ks, ke-1
229  dz(k) = cz(k+1) - cz(k)
230  end do
231 
232  ! make density at lowermost cell center
233  if ( hydrostatic_uselapserate ) then
234 
235  temp(ks) = pott_sfc - lapsdry * dz(ks-1) ! use dry lapse rate
236  pres(ks) = p00 * ( temp(ks)/pott(ks) )**(cptot/rtot)
237  dens(ks) = p00 / rtot / pott(ks) * ( pres(ks)/p00 )**(cvtot/cptot)
238 
239  else ! use itelation
240 
241  call atmos_hydrostatic_buildrho_atmos_0d( pott(ks), qv(ks), qc(ks), & ! [IN]
242  dens_sfc, pott_sfc, qv_sfc, qc_sfc, & ! [IN]
243  dz(ks-1), ks-1, & ! [IN]
244  dens(ks), temp(ks), pres(ks) ) ! [OUT]
245 
246  endif
247 
248  !--- from lowermost atmosphere to top of atmosphere
249  call atmos_hydrostatic_buildrho_atmos_1d( ka, ks, ke, &
250  pott(:), qv(:), qc(:), & ! [IN]
251  dz(:), & ! [IN]
252  dens(:), & ! [INOUT]
253  temp(:), pres(:) ) ! [OUT]
254  ! fill dummy
255  dens( 1:ks-1) = 0.0_rp
256  dens(ke+1:ka ) = 0.0_rp
257 
258  return
259  end subroutine atmos_hydrostatic_buildrho_1d
260 
261  !-----------------------------------------------------------------------------
263  subroutine atmos_hydrostatic_buildrho_3d( &
264  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
265  pott, qv, qc, &
266  pres_sfc, pott_sfc, &
267  qv_sfc, qc_sfc, &
268  cz, fz, area, &
269  dens, temp, pres, &
270  temp_sfc )
271  use scale_statistics, only: &
272  statistics_horizontal_mean
273  implicit none
274 
275  integer, intent(in) :: KA, KS, KE
276  integer, intent(in) :: IA, IS, IE
277  integer, intent(in) :: JA, JS, JE
278 
279  real(RP), intent(in) :: pott(ka,ia,ja)
280  real(RP), intent(in) :: qv (ka,ia,ja)
281  real(RP), intent(in) :: qc (ka,ia,ja)
282  real(RP), intent(in) :: pres_sfc(ia,ja)
283  real(RP), intent(in) :: pott_sfc(ia,ja)
284  real(RP), intent(in) :: qv_sfc (ia,ja)
285  real(RP), intent(in) :: qc_sfc (ia,ja)
286  real(RP), intent(in) :: cz( ka,ia,ja)
287  real(RP), intent(in) :: fz(0:ka,ia,ja)
288  real(RP), intent(in) :: area(ia,ja)
289 
290  real(RP), intent(out) :: dens(ka,ia,ja)
291  real(RP), intent(out) :: temp(ka,ia,ja)
292  real(RP), intent(out) :: pres(ka,ia,ja)
293  real(RP), intent(out) :: temp_sfc(ia,ja)
294 
295  real(RP) :: dz(ka,ia,ja), dz_top(ia,ja)
296 
297  ! TOA
298  real(RP) :: pott_toa(ia,ja)
299  real(RP) :: qv_toa (ia,ja)
300  real(RP) :: qc_toa (ia,ja)
301  real(RP) :: dens_toa(ia,ja)
302  real(RP) :: temp_toa(ia,ja)
303  real(RP) :: pres_toa(ia,ja)
304 
305  real(RP) :: dens_mean
306 
307  integer :: k, i, j
308  !---------------------------------------------------------------------------
309 
310  !--- from surface to lowermost atmosphere
311  !$omp parallel do OMP_SCHEDULE_ collapse(2)
312  do j = js, je
313  do i = is, ie
314  call atmos_hydrostatic_buildrho_1d( ka, ks, ke, &
315  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
316  pres_sfc(i,j), pott_sfc(i,j), qv_sfc(i,j), qc_sfc(i,j), & ! [IN]
317  cz(:,i,j), fz(:,i,j), & ! [IN]
318  dens(:,i,j), temp(:,i,j), pres(:,i,j), & ! [OUT]
319  temp_sfc(i,j) ) ! [OUT]
320  end do
321  end do
322 
323  !$omp parallel do OMP_SCHEDULE_ collapse(2)
324  do j = js, je
325  do i = is, ie
326  do k = ks, ke-1
327  dz(k,i,j) = cz(k+1,i,j) - cz(k,i,j)
328  end do
329  dz_top(i,j) = fz(ke,i,j) - cz(ke,i,j) ! distance from cell center to TOA
330 
331  ! value at TOA
332  pott_toa(i,j) = pott(ke,i,j)
333  qv_toa(i,j) = qv(ke,i,j)
334  qc_toa(i,j) = qc(ke,i,j)
335  end do
336  end do
337 
338  call atmos_hydrostatic_buildrho_atmos_2d( ia, is, ie, ja, js, je, &
339  pott_toa(:,:), qv_toa(:,:), qc_toa(:,:), & ! [IN]
340  dens(ke,:,:), pott(ke,:,:), qv(ke,:,:), qc(ke,:,:), & ! [IN]
341  dz_top(:,:), ke+1, & ! [IN]
342  dens_toa(:,:), temp_toa(:,:), pres_toa(:,:) ) ! [OUT]
343 
344  call statistics_horizontal_mean( ia, is, ie, ja, js, je, &
345  dens_toa(:,:), area(:,:), & ! [IN]
346  dens_mean ) ! [OUT]
347 
348  !$omp parallel do OMP_SCHEDULE_ collapse(2)
349  do j = js, je
350  do i = is, ie
351  dens_toa(i,j) = dens_mean
352  enddo
353  enddo
354 
355  call atmos_hydrostatic_buildrho_atmos_rev_2d( ia, is, ie, ja, js, je, &
356  pott(ke,:,:), qv(ke,:,:), qc(ke,:,:), & ! [IN]
357  dens_toa(:,:), pott_toa(:,:), qv_toa(:,:), qc_toa(:,:), & ! [IN]
358  dz_top(:,:), ke+1, & ! [IN]
359  dens(ke,:,:), temp(ke,:,:), pres(ke,:,:) ) ! [OUT]
360 
361  !--- from top of atmosphere to lowermost atmosphere
362  call atmos_hydrostatic_buildrho_atmos_rev_3d( ka, ks, ke, ia, is, ie, ja, js, je, &
363  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
364  dz(:,:,:), & ! [IN]
365  dens(:,:,:), & ! [INOUT]
366  temp(:,:,:), pres(:,:,:) ) ! [OUT]
367 
368  return
369  end subroutine atmos_hydrostatic_buildrho_3d
370 
371  !-----------------------------------------------------------------------------
373  subroutine atmos_hydrostatic_buildrho_real_3d( &
374  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
375  pott, qv, qc, &
376  cz, &
377  pres, &
378  dens, temp )
379  use scale_atmos_hydrometeor, only: &
380  cv_vapor, &
381  cv_water, &
382  cp_vapor, &
383  cp_water
384  implicit none
385  integer, intent(in) :: KA, KS, KE
386  integer, intent(in) :: IA, IS, IE
387  integer, intent(in) :: JA, JS, JE
388 
389  real(RP), intent(in) :: pott(ka,ia,ja)
390  real(RP), intent(in) :: qv (ka,ia,ja)
391  real(RP), intent(in) :: qc (ka,ia,ja)
392  real(RP), intent(in) :: cz (ka,ia,ja)
393 
394  real(RP), intent(inout) :: pres(ka,ia,ja)
395 
396  real(RP), intent(out) :: dens(ka,ia,ja)
397  real(RP), intent(out) :: temp(ka,ia,ja)
398 
399  real(RP) :: dz(ka,ia,ja)
400 
401  real(RP) :: pott_toa(ia,ja)
402  real(RP) :: qv_toa (ia,ja)
403  real(RP) :: qc_toa (ia,ja)
404 
405  real(RP) :: Rtot
406  real(RP) :: CVtot
407  real(RP) :: CPtot
408 
409  integer :: kref
410  integer :: k, i, j
411  !---------------------------------------------------------------------------
412 
413  !--- from surface to lowermost atmosphere
414 
415  !$omp parallel do OMP_SCHEDULE_ collapse(2)
416  do j = js, je
417  do i = is, ie
418  do k = ks, ke-1
419  dz(k,i,j) = cz(k+1,i,j) - cz(k,i,j) ! distance from cell center to cell center
420  enddo
421  enddo
422  enddo
423 
424  !$omp parallel do OMP_SCHEDULE_ collapse(2)
425  do j = js, je
426  do i = is, ie
427  pott_toa(i,j) = pott(ke,i,j)
428  qv_toa(i,j) = qv(ke,i,j)
429  qc_toa(i,j) = qc(ke,i,j)
430  enddo
431  enddo
432 
433  kref = hydrostatic_buildrho_real_kref + ks - 1
434 
435  ! calc density at reference level
436  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
437  !$omp private(Rtot,CVtot,CPtot)
438  do j = js, je
439  do i = is, ie
440  rtot = rdry * ( 1.0_rp - qv(kref,i,j) - qc(kref,i,j) ) &
441  + rvap * qv(kref,i,j)
442  cvtot = cvdry * ( 1.0_rp - qv(kref,i,j) - qc(kref,i,j) ) &
443  + cv_vapor * qv(kref,i,j) &
444  + cv_water * qc(kref,i,j)
445  cptot = cpdry * ( 1.0_rp - qv(kref,i,j) - qc(kref,i,j) ) &
446  + cp_vapor * qv(kref,i,j) &
447  + cp_water * qc(kref,i,j)
448  dens(kref,i,j) = p00 / ( rtot * pott(kref,i,j) ) * ( pres(kref,i,j)/p00 )**(cvtot/cptot)
449  enddo
450  enddo
451 
452  !--- from lowermost atmosphere to top of atmosphere
453  call atmos_hydrostatic_buildrho_atmos_3d( ka, ks, ke, ia, is, ie, ja, js, je, &
454  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
455  dz(:,:,:), & ! [IN]
456  dens(:,:,:), & ! [INOUT]
457  temp(:,:,:), pres(:,:,:), & ! [OUT]
458  kref = kref ) ! [IN]
459  call atmos_hydrostatic_buildrho_atmos_rev_3d( ka, ks, ke, ia, is, ie, ja, js, je, &
460  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
461  dz(:,:,:), & ! [IN]
462  dens(:,:,:), & ! [INOUT]
463  temp(:,:,:), pres(:,:,:), & ! [OUT]
464  kref = kref ) ! [IN]
465 
466 !!$ call ATMOS_HYDROSTATIC_buildrho_atmos_2D( IA, IS, IE, JA, JS, JE, &
467 !!$ pott_toa(:,:), qv_toa(:,:), qc_toa(:,:), & ! [IN]
468 !!$ dens(KE,:,:), & ! [IN]
469 !!$ pott(KE,:,:), qv(KE,:,:), qc(KE,:,:), & ! [IN]
470 !!$ dz(KE+1,:,:), KE+1, & ! [IN]
471 !!$ dens_toa(:,:), temp_toa(:,:), pres_toa(:,:) ) ! [OUT]
472 
473  ! density at TOA
474  dens( 1:ks-1,:,:) = 0.0_rp ! fill dummy
475 !!$ dens(KE+2:KA ,:,:) = 0.0_RP ! fill dummy
476  dens(ke+1:ka ,:,:) = 0.0_rp ! fill dummy
477 
478  return
479  end subroutine atmos_hydrostatic_buildrho_real_3d
480 
481  !-----------------------------------------------------------------------------
483 !OCL SERIAL
484 !OCL NOSIMD
485  subroutine atmos_hydrostatic_buildrho_atmos_0d( &
486  pott_L2, qv_L2, qc_L2, &
487  dens_L1, pott_L1, qv_L1, qc_L1, &
488  dz, k, &
489  dens_L2, temp_L2, pres_L2 )
490  use scale_prc, only: &
491  prc_abort
492  use scale_atmos_hydrometeor, only: &
493  cv_vapor, &
494  cv_water, &
495  cp_vapor, &
496  cp_water
497  implicit none
498  real(RP), intent(in) :: pott_L2
499  real(RP), intent(in) :: qv_L2
500  real(RP), intent(in) :: qc_L2
501  real(RP), intent(in) :: dens_L1
502  real(RP), intent(in) :: pott_L1
503  real(RP), intent(in) :: qv_L1
504  real(RP), intent(in) :: qc_L1
505  real(RP), intent(in) :: dz
506  integer, intent(in) :: k
507 
508  real(RP), intent(out) :: dens_L2
509  real(RP), intent(out) :: temp_L2
510  real(RP), intent(out) :: pres_L2
511 
512  real(RP) :: Rtot_L1 , Rtot_L2
513  real(RP) :: CVtot_L1 , CVtot_L2
514  real(RP) :: CPtot_L1 , CPtot_L2
515  real(RP) :: CPovCV_L1, CPovCV_L2
516 
517  real(RP) :: pres_L1
518  real(RP) :: dens_s, dhyd, dgrd
519  integer :: ite
520  logical :: converged
521  !---------------------------------------------------------------------------
522 
523  rtot_l1 = rdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
524  + rvap * qv_l1
525  cvtot_l1 = cvdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
526  + cv_vapor * qv_l1 &
527  + cv_water * qc_l1
528  cptot_l1 = cpdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
529  + cp_vapor * qv_l1 &
530  + cp_water * qc_l1
531  cpovcv_l1 = cptot_l1 / cvtot_l1
532 
533  rtot_l2 = rdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
534  + rvap * qv_l2
535  cvtot_l2 = cvdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
536  + cv_vapor * qv_l2 &
537  + cv_water * qc_l2
538  cptot_l2 = cpdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
539  + cp_vapor * qv_l2 &
540  + cp_water * qc_l2
541  cpovcv_l2 = cptot_l2 / cvtot_l2
542 
543  dens_s = 0.0_rp
544  dens_l2 = dens_l1 ! first guess
545 
546  pres_l1 = p00 * ( dens_l1 * rtot_l1 * pott_l1 / p00 )**cpovcv_l1
547 
548  converged = .false.
549  do ite = 1, itelim
550  if ( abs(dens_l2-dens_s) <= criteria ) then
551  converged = .true.
552  exit
553  endif
554 
555  dens_s = dens_l2
556  pres_l2 = p00 * ( dens_s * rtot_l2 * pott_l2 / p00 )**cpovcv_l2
557 
558  dhyd = + ( pres_l1 - pres_l2 ) / dz & ! dp/dz
559  - grav * 0.5_rp * ( dens_l1 + dens_s ) ! rho*g
560 
561  dgrd = - pres_l2 * cpovcv_l2 / dens_s / dz &
562  - grav * 0.5_rp
563 
564  dens_l2 = dens_s - dhyd/dgrd
565 
566  if ( dens_l2*0.0_rp /= 0.0_rp ) exit
567  enddo
568 
569  if ( .NOT. converged ) then
570  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_0D",*) 'iteration not converged!', &
571  k,dens_l2,ite,dens_s,dhyd,dgrd
572  call prc_abort
573  endif
574 
575  pres_l2 = p00 * ( dens_l2 * rtot_l2 * pott_l2 / p00 )**cpovcv_l2
576  temp_l2 = pres_l2 / ( dens_l2 * rtot_l2 )
577 
578  return
579  end subroutine atmos_hydrostatic_buildrho_atmos_0d
580 
581  !-----------------------------------------------------------------------------
583 !OCL SERIAL
584  subroutine atmos_hydrostatic_buildrho_atmos_1d( &
585  KA, KS, KE, &
586  pott, qv, qc, &
587  dz, &
588  dens, &
589  temp, pres, &
590  kref )
591  use scale_prc, only: &
592  prc_abort
593  use scale_atmos_hydrometeor, only: &
594  cv_vapor, &
595  cv_water, &
596  cp_vapor, &
597  cp_water
598  implicit none
599  integer, intent(in) :: KA, KS, KE
600 
601  real(RP), intent(in) :: pott(ka)
602  real(RP), intent(in) :: qv (ka)
603  real(RP), intent(in) :: qc (ka)
604  real(RP), intent(in) :: dz (ka)
605 
606  real(RP), intent(inout) :: dens(ka)
607 
608  real(RP), intent(out) :: temp(ka)
609  real(RP), intent(out) :: pres(ka)
610 
611  integer, intent(in), optional :: kref
612 
613  real(RP) :: Rtot (ka)
614  real(RP) :: CPovCV(ka)
615  real(RP) :: CVtot
616  real(RP) :: CPtot
617 
618  real(RP) :: dens_s, dhyd, dgrd
619  integer :: ite
620  logical :: converged
621 
622  integer :: kref_
623  integer :: k
624  !---------------------------------------------------------------------------
625 
626  if ( present(kref) ) then
627  kref_ = kref
628  else
629  kref_ = ks
630  end if
631 
632  do k = kref_, ke
633  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
634  + rvap * qv(k)
635  cvtot = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
636  + cv_vapor * qv(k) &
637  + cv_water * qc(k)
638  cptot = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
639  + cp_vapor * qv(k) &
640  + cp_water * qc(k)
641  cpovcv(k) = cptot / cvtot
642  enddo
643 
644  pres(kref_) = p00 * ( dens(kref_) * rtot(kref_) * pott(kref_) / p00 )**cpovcv(kref_)
645 
646  do k = kref_+1, ke
647 
648  dens_s = 0.0_rp
649  dens(k) = dens(k-1) ! first guess
650 
651  converged = .false.
652  do ite = 1, itelim
653  if ( abs(dens(k)-dens_s) <= criteria ) then
654  converged = .true.
655  exit
656  endif
657 
658  dens_s = dens(k)
659  pres(k) = p00 * ( dens_s * rtot(k) * pott(k) / p00 )**cpovcv(k)
660 
661  dhyd = + ( pres(k-1) - pres(k) ) / dz(k-1) & ! dpdz
662  - grav * 0.5_rp * ( dens(k-1) + dens_s ) ! rho*g
663 
664  dgrd = - cpovcv(k) * pres(k) / dens_s / dz(k-1) &
665  - grav * 0.5_rp
666 
667  dens(k) = dens_s - dhyd/dgrd
668 
669  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
670  enddo
671 
672  if ( .NOT. converged ) then
673  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_1D",*) 'iteration not converged!', &
674  k,dens(k),ite,dens_s,dhyd,dgrd
675  call prc_abort
676  endif
677 
678  pres(k) = p00 * ( dens(k) * rtot(k) * pott(k) / p00 )**cpovcv(k)
679 
680  enddo
681 
682  do k = kref_, ke
683  temp(k) = pres(k) / ( dens(k) * rtot(k) )
684  enddo
685 
686  dens( 1:ks-1) = dens(ks)
687  dens(ke+1:ka ) = dens(ke)
688  pres( 1:ks-1) = pres(ks)
689  pres(ke+1:ka ) = pres(ke)
690  temp( 1:ks-1) = temp(ks)
691  temp(ke+1:ka ) = temp(ke)
692 
693  return
694  end subroutine atmos_hydrostatic_buildrho_atmos_1d
695 
696  !-----------------------------------------------------------------------------
698 !OCL SERIAL
699  subroutine atmos_hydrostatic_buildrho_atmos_rev_1d( &
700  KA, KS, KE, &
701  pott, qv, qc, &
702  dz, &
703  dens, &
704  temp, pres, &
705  kref )
706  use scale_prc, only: &
707  prc_abort
708  use scale_atmos_hydrometeor, only: &
709  cv_vapor, &
710  cv_water, &
711  cp_vapor, &
712  cp_water
713  implicit none
714  integer, intent(in) :: KA, KS, KE
715 
716  real(RP), intent(in) :: pott(ka)
717  real(RP), intent(in) :: qv (ka)
718  real(RP), intent(in) :: qc (ka)
719  real(RP), intent(in) :: dz (ka)
720 
721  real(RP), intent(inout) :: dens(ka)
722 
723  real(RP), intent(out) :: temp(ka)
724  real(RP), intent(out) :: pres(ka)
725 
726  integer, intent(in), optional :: kref
727 
728  real(RP) :: Rtot (ka)
729  real(RP) :: CPovCV(ka)
730  real(RP) :: CVtot
731  real(RP) :: CPtot
732 
733  real(RP) :: dens_s, dhyd, dgrd
734  integer :: ite
735  logical :: converged
736 
737  integer :: kref_
738  integer :: k
739  !---------------------------------------------------------------------------
740 
741  if ( present(kref) ) then
742  kref_ = kref
743  else
744  kref_ = ke
745  end if
746 
747  do k = ks, kref_
748  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
749  + rvap * qv(k)
750  cvtot = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
751  + cv_vapor * qv(k) &
752  + cv_water * qc(k)
753  cptot = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
754  + cp_vapor * qv(k) &
755  + cp_water * qc(k)
756  cpovcv(k) = cptot / cvtot
757  enddo
758 
759  pres(kref_) = p00 * ( dens(kref_) * rtot(kref_) * pott(kref_) / p00 )**cpovcv(kref_)
760 
761  do k = kref_-1, ks, -1
762 
763  dens_s = 0.0_rp
764  dens(k) = dens(k+1) ! first guess
765 
766  converged = .false.
767  do ite = 1, itelim
768  if ( abs(dens(k)-dens_s) <= criteria ) then
769  converged = .true.
770  exit
771  endif
772 
773  dens_s = dens(k)
774  pres(k) = p00 * ( dens_s * rtot(k) * pott(k) / p00 )**cpovcv(k)
775 
776  dhyd = - ( pres(k+1) - pres(k) ) / dz(k) & ! dpdz
777  - grav * 0.5_rp * ( dens(k+1) + dens_s ) ! rho*g
778 
779  dgrd = + cpovcv(k) * pres(k) / dens_s / dz(k) &
780  - grav * 0.5_rp
781 
782  dens(k) = dens_s - dhyd/dgrd
783 
784  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
785  enddo
786 
787  if ( .NOT. converged ) then
788  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_rev_1D",*) 'iteration not converged!', &
789  k,dens(k),ite,dens_s,dhyd,dgrd
790  call prc_abort
791  endif
792 
793  pres(k) = p00 * ( dens(k) * rtot(k) * pott(k) / p00 )**cpovcv(k)
794 
795  enddo
796 
797  do k = ks, kref_
798  temp(k) = pres(k) / ( dens(k) * rtot(k) )
799  enddo
800 
801  dens( 1:ks-1) = dens(ks)
802  dens(ke+1:ka ) = dens(ke)
803  pres( 1:ks-1) = pres(ks)
804  pres(ke+1:ka ) = pres(ke)
805  temp( 1:ks-1) = temp(ks)
806  temp(ke+1:ka ) = temp(ke)
807 
808  return
809  end subroutine atmos_hydrostatic_buildrho_atmos_rev_1d
810 
811  !-----------------------------------------------------------------------------
813  subroutine atmos_hydrostatic_buildrho_atmos_2d( &
814  IA, IS, IE, JA, JS, JE, &
815  pott_L2, qv_L2, qc_L2, &
816  dens_L1, pott_L1, qv_L1, qc_L1, &
817  dz, k, &
818  dens_L2, temp_L2, pres_L2 )
819  implicit none
820  integer, intent(in) :: IA, IS, IE
821  integer, intent(in) :: JA, JS, JE
822 
823  real(RP), intent(in) :: pott_L2(ia,ja)
824  real(RP), intent(in) :: qv_L2 (ia,ja)
825  real(RP), intent(in) :: qc_L2 (ia,ja)
826  real(RP), intent(in) :: dens_L1(ia,ja)
827  real(RP), intent(in) :: pott_L1(ia,ja)
828  real(RP), intent(in) :: qv_L1 (ia,ja)
829  real(RP), intent(in) :: qc_L1 (ia,ja)
830  real(RP), intent(in) :: dz (ia,ja)
831  integer, intent(in) :: k
832 
833  real(RP), intent(out) :: dens_L2(ia,ja)
834  real(RP), intent(out) :: temp_L2(ia,ja)
835  real(RP), intent(out) :: pres_L2(ia,ja)
836 
837  integer :: i, j
838  !---------------------------------------------------------------------------
839 
840  !$omp parallel do OMP_SCHEDULE_ collapse(2)
841  do j = js, je
842  do i = is, ie
843  call atmos_hydrostatic_buildrho_atmos_0d( &
844  pott_l2(i,j), qv_l2(i,j), qc_l2(i,j), & ! [IN]
845  dens_l1(i,j), pott_l1(i,j), qv_l1(i,j), qc_l1(i,j), & ! [IN]
846  dz(i,j), k, & ! [IN]
847  dens_l2(i,j), temp_l2(i,j), pres_l2(i,j) ) ! [OUT]
848  enddo
849  enddo
850 
851  return
852  end subroutine atmos_hydrostatic_buildrho_atmos_2d
853 
854  !-----------------------------------------------------------------------------
856  subroutine atmos_hydrostatic_buildrho_atmos_rev_2d( &
857  IA, IS, IE, JA, JS, JE, &
858  pott_L1, qv_L1, qc_L1, &
859  dens_L2, pott_L2, qv_L2, qc_L2, &
860  dz, k, &
861  dens_L1, temp_L1, pres_L1 )
862  implicit none
863  integer, intent(in) :: IA, IS, IE
864  integer, intent(in) :: JA, JS, JE
865 
866  real(RP), intent(in) :: pott_L1(ia,ja)
867  real(RP), intent(in) :: qv_L1 (ia,ja)
868  real(RP), intent(in) :: qc_L1 (ia,ja)
869  real(RP), intent(in) :: dens_L2(ia,ja)
870  real(RP), intent(in) :: pott_L2(ia,ja)
871  real(RP), intent(in) :: qv_L2 (ia,ja)
872  real(RP), intent(in) :: qc_L2 (ia,ja)
873  real(RP), intent(in) :: dz (ia,ja)
874  integer, intent(in) :: k
875 
876  real(RP), intent(out) :: dens_L1(ia,ja)
877  real(RP), intent(out) :: temp_L1(ia,ja)
878  real(RP), intent(out) :: pres_L1(ia,ja)
879 
880  integer :: i, j
881  !---------------------------------------------------------------------------
882 
883  !$omp parallel do OMP_SCHEDULE_ collapse(2)
884  do j = js, je
885  do i = is, ie
886  call atmos_hydrostatic_buildrho_atmos_0d( &
887  pott_l1(i,j), qv_l1(i,j), qc_l1(i,j), & ! [IN]
888  dens_l2(i,j), pott_l2(i,j), qv_l2(i,j), qc_l2(i,j), & ! [IN]
889  -dz(i,j), k, & ! [IN]
890  dens_l1(i,j), temp_l1(i,j), pres_l1(i,j) ) ! [OUT]
891  enddo
892  enddo
893 
894  return
895  end subroutine atmos_hydrostatic_buildrho_atmos_rev_2d
896 
897  !-----------------------------------------------------------------------------
899  subroutine atmos_hydrostatic_buildrho_atmos_3d( &
900  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
901  pott, qv, qc, &
902  dz, &
903  dens, &
904  temp, pres, &
905  kref )
906  implicit none
907  integer, intent(in) :: KA, KS, KE
908  integer, intent(in) :: IA, IS, IE
909  integer, intent(in) :: JA, JS, JE
910 
911  real(RP), intent(in) :: pott(ka,ia,ja)
912  real(RP), intent(in) :: qv (ka,ia,ja)
913  real(RP), intent(in) :: qc (ka,ia,ja)
914  real(RP), intent(in) :: dz (ka,ia,ja)
915 
916  real(RP), intent(inout) :: dens(ka,ia,ja)
917 
918  real(RP), intent(out) :: temp(ka,ia,ja)
919  real(RP), intent(out) :: pres(ka,ia,ja)
920 
921  integer, intent(in), optional :: kref
922 
923  integer :: i, j
924  !---------------------------------------------------------------------------
925 
926  !$omp parallel do OMP_SCHEDULE_ collapse(2)
927  do j = js, je
928  do i = is, ie
929  call atmos_hydrostatic_buildrho_atmos_1d( ka, ks, ke, &
930  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
931  dz(:,i,j), & ! [IN]
932  dens(:,i,j), & ! [INOUT]
933  temp(:,i,j), pres(:,i,j), & ! [OUT]
934  kref=kref ) ! [IN]
935  enddo
936  enddo
937 
938  return
939  end subroutine atmos_hydrostatic_buildrho_atmos_3d
940 
941  !-----------------------------------------------------------------------------
943  subroutine atmos_hydrostatic_buildrho_atmos_rev_3d( &
944  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
945  pott, qv, qc, &
946  dz, &
947  dens, &
948  temp, pres, &
949  kref )
950  implicit none
951 
952  integer, intent(in) :: KA, KS, KE
953  integer, intent(in) :: IA, IS, IE
954  integer, intent(in) :: JA, JS, JE
955  real(RP), intent(inout) :: dens(ka,ia,ja)
956  real(RP), intent(out) :: temp(ka,ia,ja)
957  real(RP), intent(out) :: pres(ka,ia,ja)
958  real(RP), intent(in) :: pott(ka,ia,ja)
959  real(RP), intent(in) :: qv (ka,ia,ja)
960  real(RP), intent(in) :: qc (ka,ia,ja)
961  real(RP), intent(in) :: dz (ka,ia,ja)
962  integer, intent(in), optional :: kref
963 
964  integer :: i, j
965  !---------------------------------------------------------------------------
966 
967  !$omp parallel do OMP_SCHEDULE_ collapse(2)
968  do j = js, je
969  do i = is, ie
970  call atmos_hydrostatic_buildrho_atmos_rev_1d( ka, ks, ke, &
971  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
972  dz(:,i,j), & ! [IN]
973  dens(:,i,j), & ! [INOUT]
974  temp(:,i,j), pres(:,i,j), & ! [OUT]
975  kref=kref ) ! [IN]
976  enddo
977  enddo
978 
979  return
980  end subroutine atmos_hydrostatic_buildrho_atmos_rev_3d
981 
982  !-----------------------------------------------------------------------------
984 !OCL SERIAL
985  subroutine atmos_hydrostatic_buildrho_bytemp_1d( &
986  KA, KS, KE, &
987  temp, qv, qc, &
988  pres_sfc, temp_sfc, qv_sfc, qc_sfc, &
989  cz, fz, &
990  dens, pott, pres, pott_sfc )
991  use scale_prc, only: &
992  prc_abort
993  use scale_atmos_hydrometeor, only: &
994  cv_vapor, &
995  cv_water, &
996  cp_vapor, &
997  cp_water
998  implicit none
999  integer, intent(in) :: KA, KS, KE
1000 
1001  real(RP), intent(in) :: temp(ka)
1002  real(RP), intent(in) :: qv (ka)
1003  real(RP), intent(in) :: qc (ka)
1004  real(RP), intent(in) :: pres_sfc
1005  real(RP), intent(in) :: temp_sfc
1006  real(RP), intent(in) :: qv_sfc
1007  real(RP), intent(in) :: qc_sfc
1008  real(RP), intent(in) :: cz( ka)
1009  real(RP), intent(in) :: fz(0:ka)
1010 
1011  real(RP), intent(out) :: dens(ka)
1012  real(RP), intent(out) :: pott(ka)
1013  real(RP), intent(out) :: pres(ka)
1014  real(RP), intent(out) :: pott_sfc
1015 
1016  real(RP) :: dens_sfc
1017 
1018  real(RP) :: Rtot_sfc
1019  real(RP) :: CVtot_sfc
1020  real(RP) :: CPtot_sfc
1021  real(RP) :: Rtot
1022  real(RP) :: CVtot
1023  real(RP) :: CPtot
1024 
1025  real(RP) :: RovCP_sfc
1026  real(RP) :: dens_s, dhyd, dgrd
1027  integer :: ite
1028  logical :: converged
1029 
1030  real(RP) :: dz(ka)
1031 
1032  integer :: k
1033  !---------------------------------------------------------------------------
1034 
1035  !--- from surface to lowermost atmosphere
1036 
1037  rtot_sfc = rdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1038  + rvap * qv_sfc
1039  cvtot_sfc = cvdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1040  + cv_vapor * qv_sfc &
1041  + cv_water * qc_sfc
1042  cptot_sfc = cpdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1043  + cp_vapor * qv_sfc &
1044  + cp_water * qc_sfc
1045 
1046  rtot = rdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1047  + rvap * qv(ks)
1048  cvtot = cvdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1049  + cv_vapor * qv(ks) &
1050  + cv_water * qc(ks)
1051  cptot = cpdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1052  + cp_vapor * qv(ks) &
1053  + cp_water * qc(ks)
1054 
1055  ! density at surface
1056  rovcp_sfc = rtot_sfc / cptot_sfc
1057  dens_sfc = pres_sfc / ( rtot_sfc * temp_sfc )
1058  pott_sfc = temp_sfc * ( p00/pres_sfc )**rovcp_sfc
1059 
1060  ! make density at lowermost cell center
1061  dens_s = 0.0_rp
1062  dens(ks) = dens_sfc ! first guess
1063 
1064  dz(ks-1) = cz(ks) - fz(ks-1)
1065  do k = ks, ke-1
1066  dz(k) = cz(k+1) - cz(k)
1067  end do
1068 
1069  converged = .false.
1070  do ite = 1, itelim
1071  if ( abs(dens(ks)-dens_s) <= criteria ) then
1072  converged = .true.
1073  exit
1074  endif
1075 
1076  dens_s = dens(ks)
1077 
1078  dhyd = + ( pres_sfc - dens_s * rtot * temp(ks) ) / dz(ks-1) & ! dp/dz
1079  - grav * 0.5_rp * ( dens_sfc + dens_s ) ! rho*g
1080 
1081  dgrd = - rtot * temp(ks) / dz(ks-1) &
1082  - grav * 0.5_rp
1083 
1084  dens(ks) = dens_s - dhyd/dgrd
1085 
1086  if ( dens(ks)*0.0_rp /= 0.0_rp ) exit
1087  enddo
1088 
1089  if ( .NOT. converged ) then
1090  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_1D",*) 'iteration not converged!', &
1091  dens(ks),ite,dens_s,dhyd,dgrd
1092  call prc_abort
1093  endif
1094 
1095  !--- from lowermost atmosphere to top of atmosphere
1096  call atmos_hydrostatic_buildrho_bytemp_atmos_1d( ka, ks, ke, &
1097  temp(:), qv(:), qc(:), & ! [IN]
1098  dz(:), & ! [IN]
1099  dens(:), & ! [INOUT]
1100  pott(:), pres(:) ) ! [OUT]
1101 
1102  return
1103  end subroutine atmos_hydrostatic_buildrho_bytemp_1d
1104 
1105  !-----------------------------------------------------------------------------
1107  subroutine atmos_hydrostatic_buildrho_bytemp_3d( &
1108  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1109  temp, qv, qc, &
1110  pres_sfc, temp_sfc, qv_sfc, qc_sfc, &
1111  cz, fz, &
1112  dens, pott, pres, &
1113  pott_sfc )
1114  implicit none
1115  integer, intent(in) :: KA, KS, KE
1116  integer, intent(in) :: IA, IS, IE
1117  integer, intent(in) :: JA, JS, JE
1118 
1119  real(RP), intent(in) :: temp(ka,ia,ja)
1120  real(RP), intent(in) :: qv (ka,ia,ja)
1121  real(RP), intent(in) :: qc (ka,ia,ja)
1122  real(RP), intent(in) :: pres_sfc(ia,ja)
1123  real(RP), intent(in) :: temp_sfc(ia,ja)
1124  real(RP), intent(in) :: qv_sfc (ia,ja)
1125  real(RP), intent(in) :: qc_sfc (ia,ja)
1126  real(RP), intent(in) :: cz(ka,ia,ja)
1127  real(RP), intent(in) :: fz(ka,ia,ja)
1128 
1129  real(RP), intent(out) :: dens(ka,ia,ja)
1130  real(RP), intent(out) :: pott(ka,ia,ja)
1131  real(RP), intent(out) :: pres(ka,ia,ja)
1132  real(RP), intent(out) :: pott_sfc(ia,ja)
1133 
1134  integer :: i, j
1135  !---------------------------------------------------------------------------
1136 
1137  !--- from surface to lowermost atmosphere
1138  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1139  do j = js, je
1140  do i = is, ie
1141  call atmos_hydrostatic_buildrho_bytemp_1d( ka, ks, ke, &
1142  temp(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
1143  pres_sfc(i,j), temp_sfc(i,j), qv_sfc(i,j), qc_sfc(i,j), & ! [IN]
1144  cz(:,i,j), fz(:,i,j), & ! [IN]
1145  dens(:,i,j), pott(:,i,j), pres(:,i,j), & ! [OUT]
1146  pott_sfc(i,j) ) ! [OUT]
1147  enddo
1148  enddo
1149 
1150  return
1151  end subroutine atmos_hydrostatic_buildrho_bytemp_3d
1152 
1153  !-----------------------------------------------------------------------------
1155 !OCL SERIAL
1156  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_1d( &
1157  KA, KS, KE, &
1158  temp, qv, qc, &
1159  dz, &
1160  dens, &
1161  pott, pres )
1162  use scale_prc, only: &
1163  prc_abort
1164  use scale_atmos_hydrometeor, only: &
1165  cv_vapor, &
1166  cv_water, &
1167  cp_vapor, &
1168  cp_water
1169  implicit none
1170  integer, intent(in) :: KA, KS, KE
1171 
1172  real(RP), intent(in) :: temp(ka)
1173  real(RP), intent(in) :: qv (ka)
1174  real(RP), intent(in) :: qc (ka)
1175  real(RP), intent(in) :: dz (ka)
1176 
1177  real(RP), intent(inout) :: dens(ka)
1178 
1179  real(RP), intent(out) :: pott(ka)
1180  real(RP), intent(out) :: pres(ka)
1181 
1182  real(RP) :: Rtot (ka)
1183  real(RP) :: CVtot (ka)
1184  real(RP) :: CPtot (ka)
1185 
1186  real(RP) :: RovCP
1187  real(RP) :: dens_s, dhyd, dgrd
1188  integer :: ite
1189  logical :: converged
1190 
1191  integer :: k
1192  !---------------------------------------------------------------------------
1193 
1194  do k = ks, ke
1195  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
1196  + rvap * qv(k)
1197  cvtot(k) = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
1198  + cv_vapor * qv(k) &
1199  + cv_water * qc(k)
1200  cptot(k) = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
1201  + cp_vapor * qv(k) &
1202  + cp_water * qc(k)
1203  enddo
1204 
1205  pres(ks) = dens(ks) * rtot(ks) * temp(ks)
1206 
1207  do k = ks+1, ke
1208 
1209  dens_s = 0.0_rp
1210  dens(k) = dens(k-1) ! first guess
1211 
1212  converged = .false.
1213  do ite = 1, itelim
1214  if ( abs(dens(k)-dens_s) <= criteria ) then
1215  converged = .true.
1216  exit
1217  endif
1218 
1219  dens_s = dens(k)
1220 
1221  dhyd = + ( pres(k-1) - dens_s * rtot(k) * temp(k) ) / dz(k-1) & ! dp/dz
1222  - grav * 0.5_rp * ( dens(k-1) + dens_s ) ! rho*g
1223 
1224  dgrd = - rtot(k) * temp(k) / dz(k-1) &
1225  - 0.5_rp * grav
1226 
1227  dens(k) = dens_s - dhyd/dgrd
1228 
1229  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
1230  enddo
1231 
1232  if ( .NOT. converged ) then
1233  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_atmos_1D",*) 'iteration not converged!', &
1234  k,dens(k),ite,dens_s,dhyd,dgrd
1235  call prc_abort
1236  endif
1237 
1238  pres(k) = dens(k) * rtot(k) * temp(k)
1239 
1240  enddo
1241 
1242  do k = ks, ke
1243  rovcp = rtot(k) / cptot(k)
1244  pott(k) = temp(k) * ( p00 / pres(k) )**rovcp
1245  enddo
1246 
1247  return
1248  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_1d
1249 
1250  !-----------------------------------------------------------------------------
1252 !OCL SERIAL
1253  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d( &
1254  KA, KS, KE, &
1255  temp, qv, qc, &
1256  dz, &
1257  dens, &
1258  pott, pres )
1259  use scale_prc, only: &
1260  prc_abort
1261  use scale_atmos_hydrometeor, only: &
1262  cv_vapor, &
1263  cv_water, &
1264  cp_vapor, &
1265  cp_water
1266  implicit none
1267  integer, intent(in) :: KA, KS, KE
1268 
1269  real(RP), intent(in) :: temp(ka)
1270  real(RP), intent(in) :: qv (ka)
1271  real(RP), intent(in) :: qc (ka)
1272  real(RP), intent(in) :: dz (ka)
1273 
1274  real(RP), intent(inout) :: dens(ka)
1275 
1276  real(RP), intent(out) :: pott(ka)
1277  real(RP), intent(out) :: pres(ka)
1278 
1279  real(RP) :: Rtot (ka)
1280  real(RP) :: CVtot (ka)
1281  real(RP) :: CPtot (ka)
1282 
1283  real(RP) :: RovCP
1284  real(RP) :: dens_s, dhyd, dgrd
1285  integer :: ite
1286  logical :: converged
1287 
1288  integer :: k
1289  !---------------------------------------------------------------------------
1290 
1291  do k = ks, ke
1292  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
1293  + rvap * qv(k)
1294  cvtot(k) = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
1295  + cv_vapor * qv(k) &
1296  + cv_water * qc(k)
1297  cptot(k) = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
1298  + cp_vapor * qv(k) &
1299  + cp_water * qc(k)
1300  enddo
1301 
1302  pres(ke) = dens(ke) * rtot(ke) * temp(ke)
1303 
1304  do k = ke-1, ke, -1
1305 
1306  dens_s = 0.0_rp
1307  dens(k) = dens(k+1) ! first guess
1308 
1309  converged = .false.
1310  do ite = 1, itelim
1311  if ( abs(dens(k)-dens_s) <= criteria ) then
1312  converged = .true.
1313  exit
1314  endif
1315 
1316  dens_s = dens(k)
1317 
1318  dhyd = - ( pres(k+1) - dens_s * rtot(k) * temp(k) ) / dz(k) & ! dp/dz
1319  - grav * 0.5_rp * ( dens(k+1) + dens_s ) ! rho*g
1320 
1321  dgrd = + rtot(k) * temp(k) / dz(k) &
1322  - 0.5_rp * grav
1323 
1324  dens(k) = dens_s - dhyd/dgrd
1325 
1326  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
1327  enddo
1328 
1329  if ( .NOT. converged ) then
1330  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_atmos_rev_1D",*) 'iteration not converged!', &
1331  k,dens(k),ite,dens_s,dhyd,dgrd
1332  call prc_abort
1333  endif
1334 
1335  pres(k) = dens(k) * rtot(k) * temp(k)
1336 
1337  enddo
1338 
1339  do k = ks, ke
1340  rovcp = rtot(k) / cptot(k)
1341  pott(k) = temp(k) * ( p00 / pres(k) )**rovcp
1342  enddo
1343 
1344  return
1345  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d
1346 
1347  !-----------------------------------------------------------------------------
1349  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_3d( &
1350  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1351  temp, qv, qc, &
1352  dz, &
1353  dens, &
1354  pott, pres )
1355  implicit none
1356  integer, intent(in) :: KA, KS, KE
1357  integer, intent(in) :: IA, IS, IE
1358  integer, intent(in) :: JA, JS, JE
1359 
1360  real(RP), intent(in) :: temp(ka,ia,ja)
1361  real(RP), intent(in) :: qv (ka,ia,ja)
1362  real(RP), intent(in) :: qc (ka,ia,ja)
1363  real(RP), intent(in) :: dz (ka,ia,ja)
1364 
1365  real(RP), intent(inout) :: dens(ka,ia,ja)
1366  real(RP), intent(out) :: pott(ka,ia,ja)
1367  real(RP), intent(out) :: pres(ka,ia,ja)
1368 
1369  integer :: i, j
1370  !---------------------------------------------------------------------------
1371 
1372  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1373  do j = js, je
1374  do i = is, ie
1375  call atmos_hydrostatic_buildrho_bytemp_atmos_1d( ka, ks, ke, &
1376  temp(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
1377  dz(:,i,j), & ! [IN]
1378  dens(:,i,j), & ! [INOUT]
1379  pott(:,i,j), pres(:,i,j) ) ! [OUT]
1380  enddo
1381  enddo
1382 
1383  return
1384  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_3d
1385 
1386  !-----------------------------------------------------------------------------
1388 !OCL SERIAL
1389 !OCL NOSIMD
1390  subroutine atmos_hydrostatic_barometric_law_mslp_0d( &
1391  KA, KS, KE, &
1392  pres, temp, qv, &
1393  cz, &
1394  mslp )
1395  implicit none
1396  integer, intent(in) :: KA, KS, KE
1397 
1398  real(RP), intent(in) :: pres(ka)
1399  real(RP), intent(in) :: temp(ka)
1400  real(RP), intent(in) :: qv (ka)
1401  real(RP), intent(in) :: cz (ka)
1402 
1403  real(RP), intent(out) :: mslp
1404 
1405  ! work
1406  integer :: kref
1407  real(RP) :: vtemp
1408  !---------------------------------------------------------------------------
1409 
1410  kref = hydrostatic_barometric_law_mslp_kref + ks - 1
1411 
1412  ! virtual temperature
1413  vtemp = temp(kref) * (1.0_rp + epstvap * qv(kref) )
1414 
1415  ! barometric law assuming constant lapse rate
1416  mslp = pres(kref) * ( ( vtemp + laps * cz(kref) ) / vtemp ) ** ( grav / ( rdry * laps ) )
1417 
1418  return
1419  end subroutine atmos_hydrostatic_barometric_law_mslp_0d
1420 
1421  !-----------------------------------------------------------------------------
1423  subroutine atmos_hydrostatic_barometric_law_mslp_2d( &
1424  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1425  pres, temp, qv, &
1426  cz, &
1427  mslp )
1428  implicit none
1429  integer, intent(in) :: KA, KS, KE
1430  integer, intent(in) :: IA, IS, IE
1431  integer, intent(in) :: JA, JS, JE
1432 
1433  real(RP), intent(in) :: pres(ka,ia,ja)
1434  real(RP), intent(in) :: temp(ka,ia,ja)
1435  real(RP), intent(in) :: qv (ka,ia,ja)
1436  real(RP), intent(in) :: cz (ka,ia,ja)
1437 
1438  real(RP), intent(out) :: mslp(ia,ja)
1439 
1440  integer :: i, j
1441  !---------------------------------------------------------------------------
1442 
1443  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1444  do j = js, je
1445  do i = is, ie
1446  call atmos_hydrostatic_barometric_law_mslp_0d( ka, ks, ke, &
1447  pres(:,i,j), temp(:,i,j), qv(:,i,j), & ! [IN]
1448  cz(:,i,j), & ! [IN]
1449  mslp(i,j) ) ! [OUT]
1450  enddo
1451  enddo
1452 
1453  return
1454  end subroutine atmos_hydrostatic_barometric_law_mslp_2d
1455 
1456  !-----------------------------------------------------------------------------
1458 !OCL SERIAL
1459 !OCL NOSIMD
1460  subroutine atmos_hydrostatic_barometric_law_pres_0d( &
1461  mslp, temp, &
1462  dz, &
1463  pres )
1464  implicit none
1465  real(RP), intent(in) :: mslp
1466  real(RP), intent(in) :: temp
1467  real(RP), intent(in) :: dz
1468 
1469  real(RP), intent(out) :: pres
1470 
1471  ! work
1472  real(RP) :: TM
1473  !---------------------------------------------------------------------------
1474 
1475  tm = temp + laps * dz * 0.5_rp ! column-mean air temperature
1476 
1477  ! barometric law
1478  pres = mslp / exp( grav * dz / ( rdry * tm ) )
1479 
1480  return
1481  end subroutine atmos_hydrostatic_barometric_law_pres_0d
1482 
1483  !-----------------------------------------------------------------------------
1485  subroutine atmos_hydrostatic_barometric_law_pres_2d( &
1486  IA, IS, IE, JA, JS, JE, &
1487  mslp, temp, &
1488  dz, &
1489  pres )
1490  implicit none
1491  integer, intent(in) :: IA, IS, IE
1492  integer, intent(in) :: JA, JS, JE
1493 
1494  real(RP), intent(in) :: mslp(ia,ja)
1495  real(RP), intent(in) :: temp(ia,ja)
1496  real(RP), intent(in) :: dz (ia,ja)
1497 
1498  real(RP), intent(out) :: pres(ia,ja)
1499 
1500  integer :: i, j
1501  !---------------------------------------------------------------------------
1502 
1503  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1504  do j = js, je
1505  do i = is, ie
1506  call atmos_hydrostatic_barometric_law_pres_0d( mslp(i,j), temp(i,j), dz(i,j), & ! [IN]
1507  pres(i,j) ) ! [OUT]
1508  enddo
1509  enddo
1510 
1511  return
1512  end subroutine atmos_hydrostatic_barometric_law_pres_2d
1513 
1514 end module scale_atmos_hydrostatic
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:57
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d
subroutine atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d(KA, KS, KE, temp, qv, qc, dz, dens, pott, pres)
Build up density from upermost atmosphere (1D)
Definition: scale_atmos_hydrostatic.F90:1259
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_const::const_cvvap
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
Definition: scale_const.F90:65
scale_const::const_laps
real(rp), public const_laps
lapse rate of ISA [K/m]
Definition: scale_const.F90:58
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_real_3d
subroutine atmos_hydrostatic_buildrho_real_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, pott, qv, qc, cz, pres, dens, temp)
Build up density from surface (3D), not to reverse from TOA.
Definition: scale_atmos_hydrostatic.F90:379
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:130
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_bytemp_atmos_3d
subroutine atmos_hydrostatic_buildrho_bytemp_atmos_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, temp, qv, qc, dz, dens, pott, pres)
Build up density from lowermost atmosphere (3D)
Definition: scale_atmos_hydrostatic.F90:1355
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_atmos_0d
subroutine atmos_hydrostatic_buildrho_atmos_0d(pott_L2, qv_L2, qc_L2, dens_L1, pott_L1, qv_L1, qc_L1, dz, k, dens_L2, temp_L2, pres_L2)
Build up density (0D)
Definition: scale_atmos_hydrostatic.F90:490
scale_const::const_lapsdry
real(rp), public const_lapsdry
dry adiabatic lapse rate [K/m]
Definition: scale_const.F90:59
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_hydrostatic
module atmosphere / hydrostatic barance
Definition: scale_atmos_hydrostatic.F90:12
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_3d
subroutine atmos_hydrostatic_buildrho_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, pott, qv, qc, pres_sfc, pott_sfc, qv_sfc, qc_sfc, cz, fz, area, dens, temp, pres, temp_sfc)
Build up density from surface (3D)
Definition: scale_atmos_hydrostatic.F90:271
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_bytemp_3d
subroutine atmos_hydrostatic_buildrho_bytemp_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, temp, qv, qc, pres_sfc, temp_sfc, qv_sfc, qc_sfc, cz, fz, dens, pott, pres, pott_sfc)
Build up density from surface (3D)
Definition: scale_atmos_hydrostatic.F90:1114
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_atmos_1d
subroutine atmos_hydrostatic_buildrho_atmos_1d(KA, KS, KE, pott, qv, qc, dz, dens, temp, pres, kref)
Build up density from lowermost atmosphere (1D)
Definition: scale_atmos_hydrostatic.F90:591
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_atmos_rev_1d
subroutine atmos_hydrostatic_buildrho_atmos_rev_1d(KA, KS, KE, pott, qv, qc, dz, dens, temp, pres, kref)
Build up density from upermost atmosphere (1D)
Definition: scale_atmos_hydrostatic.F90:706
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_1d
subroutine atmos_hydrostatic_buildrho_1d(KA, KS, KE, pott, qv, qc, pres_sfc, pott_sfc, qv_sfc, qc_sfc, cz, fz, dens, temp, pres, temp_sfc)
Build up density from surface (1D)
Definition: scale_atmos_hydrostatic.F90:163
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:64
scale_atmos_hydrometeor::cp_vapor
real(rp), public cp_vapor
CP for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:131
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:133
scale_atmos_hydrostatic::atmos_hydrostatic_setup
subroutine, public atmos_hydrostatic_setup
Setup.
Definition: scale_atmos_hydrostatic.F90:116
scale_atmos_hydrometeor::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:132
Generated by   doxygen 1.8.13