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  use scale_const, only: &
385  undef => const_undef
386  implicit none
387  integer, intent(in) :: KA, KS, KE
388  integer, intent(in) :: IA, IS, IE
389  integer, intent(in) :: JA, JS, JE
390 
391  real(RP), intent(in) :: pott(KA,IA,JA)
392  real(RP), intent(in) :: qv (KA,IA,JA)
393  real(RP), intent(in) :: qc (KA,IA,JA)
394  real(RP), intent(in) :: cz (KA,IA,JA)
395 
396  real(RP), intent(inout) :: pres(KA,IA,JA)
397 
398  real(RP), intent(out) :: dens(KA,IA,JA)
399  real(RP), intent(out) :: temp(KA,IA,JA)
400 
401  real(RP) :: dz(KA,IA,JA)
402 
403  real(RP) :: pott_toa(IA,JA)
404  real(RP) :: qv_toa (IA,JA)
405  real(RP) :: qc_toa (IA,JA)
406 
407  real(RP) :: Rtot
408  real(RP) :: CVtot
409  real(RP) :: CPtot
410 
411  integer :: kref(IA,JA)
412  integer :: k, i, j
413  !---------------------------------------------------------------------------
414 
415  !--- from surface to lowermost atmosphere
416 
417  !$omp parallel do OMP_SCHEDULE_ collapse(2)
418  do j = js, je
419  do i = is, ie
420  do k = ks, ke-1
421  dz(k,i,j) = cz(k+1,i,j) - cz(k,i,j) ! distance from cell center to cell center
422  enddo
423  enddo
424  enddo
425 
426  !$omp parallel do OMP_SCHEDULE_ collapse(2)
427  do j = js, je
428  do i = is, ie
429  pott_toa(i,j) = pott(ke,i,j)
430  qv_toa(i,j) = qv(ke,i,j)
431  qc_toa(i,j) = qc(ke,i,j)
432  enddo
433  enddo
434 
435  !$omp parallel do OMP_SCHEDULE_ collapse(2)
436  do j = js, je
437  do i = is, ie
438  kref(i,j) = hydrostatic_buildrho_real_kref + ks - 1
439  do k = kref(i,j), ke
440  if ( pres(k,i,j) .ne. undef ) then
441  kref(i,j) = k
442  exit
443  end if
444  end do
445  end do
446  end do
447 
448  ! calc density at reference level
449  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
450  !$omp private(Rtot,CVtot,CPtot,k)
451  do j = js, je
452  do i = is, ie
453  k = kref(i,j)
454  rtot = rdry * ( 1.0_rp - qv(k,i,j) - qc(k,i,j) ) &
455  + rvap * qv(k,i,j)
456  cvtot = cvdry * ( 1.0_rp - qv(k,i,j) - qc(k,i,j) ) &
457  + cv_vapor * qv(k,i,j) &
458  + cv_water * qc(k,i,j)
459  cptot = cpdry * ( 1.0_rp - qv(k,i,j) - qc(k,i,j) ) &
460  + cp_vapor * qv(k,i,j) &
461  + cp_water * qc(k,i,j)
462  dens(k,i,j) = p00 / ( rtot * pott(k,i,j) ) * ( pres(k,i,j)/p00 )**(cvtot/cptot)
463  enddo
464  enddo
465 
466  !--- from lowermost atmosphere to top of atmosphere
467  call atmos_hydrostatic_buildrho_atmos_3d( ka, ks, ke, ia, is, ie, ja, js, je, &
468  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
469  dz(:,:,:), & ! [IN]
470  dens(:,:,:), & ! [INOUT]
471  temp(:,:,:), pres(:,:,:), & ! [OUT]
472  kref = kref ) ! [IN]
473  call atmos_hydrostatic_buildrho_atmos_rev_3d( ka, ks, ke, ia, is, ie, ja, js, je, &
474  pott(:,:,:), qv(:,:,:), qc(:,:,:), & ! [IN]
475  dz(:,:,:), & ! [IN]
476  dens(:,:,:), & ! [INOUT]
477  temp(:,:,:), pres(:,:,:), & ! [OUT]
478  kref = kref ) ! [IN]
479 
480 !!$ call ATMOS_HYDROSTATIC_buildrho_atmos_2D( IA, IS, IE, JA, JS, JE, &
481 !!$ pott_toa(:,:), qv_toa(:,:), qc_toa(:,:), & ! [IN]
482 !!$ dens(KE,:,:), & ! [IN]
483 !!$ pott(KE,:,:), qv(KE,:,:), qc(KE,:,:), & ! [IN]
484 !!$ dz(KE+1,:,:), KE+1, & ! [IN]
485 !!$ dens_toa(:,:), temp_toa(:,:), pres_toa(:,:) ) ! [OUT]
486 
487  ! density at TOA
488  dens( 1:ks-1,:,:) = 0.0_rp ! fill dummy
489 !!$ dens(KE+2:KA ,:,:) = 0.0_RP ! fill dummy
490  dens(ke+1:ka ,:,:) = 0.0_rp ! fill dummy
491 
492  return
493  end subroutine atmos_hydrostatic_buildrho_real_3d
494 
495  !-----------------------------------------------------------------------------
497 !OCL SERIAL
498 !OCL NOSIMD
499  subroutine atmos_hydrostatic_buildrho_atmos_0d( &
500  pott_L2, qv_L2, qc_L2, &
501  dens_L1, pott_L1, qv_L1, qc_L1, &
502  dz, k, &
503  dens_L2, temp_L2, pres_L2 )
504  use scale_prc, only: &
505  prc_abort
506  use scale_atmos_hydrometeor, only: &
507  cv_vapor, &
508  cv_water, &
509  cp_vapor, &
510  cp_water
511  implicit none
512  real(RP), intent(in) :: pott_L2
513  real(RP), intent(in) :: qv_L2
514  real(RP), intent(in) :: qc_L2
515  real(RP), intent(in) :: dens_L1
516  real(RP), intent(in) :: pott_L1
517  real(RP), intent(in) :: qv_L1
518  real(RP), intent(in) :: qc_L1
519  real(RP), intent(in) :: dz
520  integer, intent(in) :: k
521 
522  real(RP), intent(out) :: dens_L2
523  real(RP), intent(out) :: temp_L2
524  real(RP), intent(out) :: pres_L2
525 
526  real(RP) :: Rtot_L1 , Rtot_L2
527  real(RP) :: CVtot_L1 , CVtot_L2
528  real(RP) :: CPtot_L1 , CPtot_L2
529  real(RP) :: CPovCV_L1, CPovCV_L2
530 
531  real(RP) :: pres_L1
532  real(RP) :: dens_s, dhyd, dgrd
533  integer :: ite
534  logical :: converged
535  !---------------------------------------------------------------------------
536 
537  rtot_l1 = rdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
538  + rvap * qv_l1
539  cvtot_l1 = cvdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
540  + cv_vapor * qv_l1 &
541  + cv_water * qc_l1
542  cptot_l1 = cpdry * ( 1.0_rp - qv_l1 - qc_l1 ) &
543  + cp_vapor * qv_l1 &
544  + cp_water * qc_l1
545  cpovcv_l1 = cptot_l1 / cvtot_l1
546 
547  rtot_l2 = rdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
548  + rvap * qv_l2
549  cvtot_l2 = cvdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
550  + cv_vapor * qv_l2 &
551  + cv_water * qc_l2
552  cptot_l2 = cpdry * ( 1.0_rp - qv_l2 - qc_l2 ) &
553  + cp_vapor * qv_l2 &
554  + cp_water * qc_l2
555  cpovcv_l2 = cptot_l2 / cvtot_l2
556 
557  dens_s = 0.0_rp
558  dens_l2 = dens_l1 ! first guess
559 
560  pres_l1 = p00 * ( dens_l1 * rtot_l1 * pott_l1 / p00 )**cpovcv_l1
561 
562  converged = .false.
563  do ite = 1, itelim
564  if ( abs(dens_l2-dens_s) <= criteria ) then
565  converged = .true.
566  exit
567  endif
568 
569  dens_s = dens_l2
570  pres_l2 = p00 * ( dens_s * rtot_l2 * pott_l2 / p00 )**cpovcv_l2
571 
572  dhyd = + ( pres_l1 - pres_l2 ) / dz & ! dp/dz
573  - grav * 0.5_rp * ( dens_l1 + dens_s ) ! rho*g
574 
575  dgrd = - pres_l2 * cpovcv_l2 / dens_s / dz &
576  - grav * 0.5_rp
577 
578  dens_l2 = dens_s - dhyd/dgrd
579 
580  if ( dens_l2*0.0_rp /= 0.0_rp ) exit
581  enddo
582 
583  if ( .NOT. converged ) then
584  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_0D",*) 'iteration not converged!', &
585  k,dens_l2,ite,dens_s,dhyd,dgrd
586  call prc_abort
587  endif
588 
589  pres_l2 = p00 * ( dens_l2 * rtot_l2 * pott_l2 / p00 )**cpovcv_l2
590  temp_l2 = pres_l2 / ( dens_l2 * rtot_l2 )
591 
592  return
593  end subroutine atmos_hydrostatic_buildrho_atmos_0d
594 
595  !-----------------------------------------------------------------------------
597 !OCL SERIAL
598  subroutine atmos_hydrostatic_buildrho_atmos_1d( &
599  KA, KS, KE, &
600  pott, qv, qc, &
601  dz, &
602  dens, &
603  temp, pres )
604  use scale_const, only: &
605  undef => const_undef
606  use scale_prc, only: &
607  prc_abort
608  use scale_atmos_hydrometeor, only: &
609  cv_vapor, &
610  cv_water, &
611  cp_vapor, &
612  cp_water
613  implicit none
614  integer, intent(in) :: KA, KS, KE
615 
616  real(RP), intent(in) :: pott(KA)
617  real(RP), intent(in) :: qv (KA)
618  real(RP), intent(in) :: qc (KA)
619  real(RP), intent(in) :: dz (KA)
620 
621  real(RP), intent(inout) :: dens(KA)
622 
623  real(RP), intent(out) :: temp(KA)
624  real(RP), intent(out) :: pres(KA)
625 
626  real(RP) :: Rtot (KA)
627  real(RP) :: CPovCV(KA)
628  real(RP) :: CVtot
629  real(RP) :: CPtot
630 
631  real(RP) :: dens_s, dhyd, dgrd
632  integer :: ite
633  logical :: converged
634 
635  integer :: k
636  !---------------------------------------------------------------------------
637 
638  do k = ks, ke
639  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
640  + rvap * qv(k)
641  cvtot = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
642  + cv_vapor * qv(k) &
643  + cv_water * qc(k)
644  cptot = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
645  + cp_vapor * qv(k) &
646  + cp_water * qc(k)
647  cpovcv(k) = cptot / cvtot
648  enddo
649 
650  pres(ks) = p00 * ( dens(ks) * rtot(ks) * pott(ks) / p00 )**cpovcv(ks)
651 
652  do k = ks+1, ke
653 
654  dens_s = 0.0_rp
655  dens(k) = dens(k-1) ! first guess
656 
657  converged = .false.
658  do ite = 1, itelim
659  if ( abs(dens(k)-dens_s) <= criteria ) then
660  converged = .true.
661  exit
662  endif
663 
664  dens_s = dens(k)
665  pres(k) = p00 * ( dens_s * rtot(k) * pott(k) / p00 )**cpovcv(k)
666 
667  dhyd = + ( pres(k-1) - pres(k) ) / dz(k-1) & ! dpdz
668  - grav * 0.5_rp * ( dens(k-1) + dens_s ) ! rho*g
669 
670  dgrd = - cpovcv(k) * pres(k) / dens_s / dz(k-1) &
671  - grav * 0.5_rp
672 
673  dens(k) = dens_s - dhyd/dgrd
674 
675  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
676  enddo
677 
678  if ( .NOT. converged ) then
679  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_1D",*) 'iteration not converged!', &
680  k,dens(k),ite,dens_s,dhyd,dgrd
681  call prc_abort
682  endif
683 
684  pres(k) = p00 * ( dens(k) * rtot(k) * pott(k) / p00 )**cpovcv(k)
685 
686  enddo
687 
688  do k = ks, ke
689  temp(k) = pres(k) / ( dens(k) * rtot(k) )
690  enddo
691 
692  dens(ke+1:ka ) = dens(ke)
693  pres(ke+1:ka ) = pres(ke)
694  temp(ke+1:ka ) = temp(ke)
695 
696  return
697  end subroutine atmos_hydrostatic_buildrho_atmos_1d
698 
699  !-----------------------------------------------------------------------------
701 !OCL SERIAL
702  subroutine atmos_hydrostatic_buildrho_atmos_rev_1d( &
703  KA, KS, KE, &
704  pott, qv, qc, &
705  dz, &
706  dens, &
707  temp, pres )
708  use scale_const, only: &
709  undef => const_undef
710  use scale_prc, only: &
711  prc_abort
712  use scale_atmos_hydrometeor, only: &
713  cv_vapor, &
714  cv_water, &
715  cp_vapor, &
716  cp_water
717  implicit none
718  integer, intent(in) :: KA, KS, KE
719 
720  real(RP), intent(in) :: pott(KA)
721  real(RP), intent(in) :: qv (KA)
722  real(RP), intent(in) :: qc (KA)
723  real(RP), intent(in) :: dz (KA)
724 
725  real(RP), intent(inout) :: dens(KA)
726 
727  real(RP), intent(out) :: temp(KA)
728  real(RP), intent(out) :: pres(KA)
729 
730  real(RP) :: Rtot (KA)
731  real(RP) :: CPovCV(KA)
732  real(RP) :: CVtot
733  real(RP) :: CPtot
734 
735  real(RP) :: dens_s, dhyd, dgrd
736  integer :: ite
737  logical :: converged
738 
739  integer :: k
740  !---------------------------------------------------------------------------
741 
742  do k = ks, ke
743  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
744  + rvap * qv(k)
745  cvtot = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
746  + cv_vapor * qv(k) &
747  + cv_water * qc(k)
748  cptot = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
749  + cp_vapor * qv(k) &
750  + cp_water * qc(k)
751  cpovcv(k) = cptot / cvtot
752  enddo
753 
754  pres(ke) = p00 * ( dens(ke) * rtot(ke) * pott(ke) / p00 )**cpovcv(ke)
755 
756  do k = ke-1, ks, -1
757 
758  dens_s = 0.0_rp
759  dens(k) = dens(k+1) ! first guess
760 
761  converged = .false.
762  do ite = 1, itelim
763  if ( abs(dens(k)-dens_s) <= criteria ) then
764  converged = .true.
765  exit
766  endif
767 
768  dens_s = dens(k)
769  pres(k) = p00 * ( dens_s * rtot(k) * pott(k) / p00 )**cpovcv(k)
770 
771  dhyd = - ( pres(k+1) - pres(k) ) / dz(k) & ! dpdz
772  - grav * 0.5_rp * ( dens(k+1) + dens_s ) ! rho*g
773 
774  dgrd = + cpovcv(k) * pres(k) / dens_s / dz(k) &
775  - grav * 0.5_rp
776 
777  dens(k) = dens_s - dhyd/dgrd
778 
779  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
780  enddo
781 
782  if ( .NOT. converged ) then
783  log_error("ATMOS_HYDROSTATIC_buildrho_atmos_rev_1D",*) 'iteration not converged!', &
784  k,dens(k),ite,dens_s,dhyd,dgrd
785  call prc_abort
786  endif
787 
788  pres(k) = p00 * ( dens(k) * rtot(k) * pott(k) / p00 )**cpovcv(k)
789 
790  enddo
791 
792  do k = ks, ke
793  temp(k) = pres(k) / ( dens(k) * rtot(k) )
794  enddo
795 
796  dens( 1:ks-1) = dens(ks)
797  pres( 1:ks-1) = pres(ks)
798  temp( 1:ks-1) = temp(ks)
799 
800  return
801  end subroutine atmos_hydrostatic_buildrho_atmos_rev_1d
802 
803  !-----------------------------------------------------------------------------
805  subroutine atmos_hydrostatic_buildrho_atmos_2d( &
806  IA, IS, IE, JA, JS, JE, &
807  pott_L2, qv_L2, qc_L2, &
808  dens_L1, pott_L1, qv_L1, qc_L1, &
809  dz, k, &
810  dens_L2, temp_L2, pres_L2 )
811  implicit none
812  integer, intent(in) :: IA, IS, IE
813  integer, intent(in) :: JA, JS, JE
814 
815  real(RP), intent(in) :: pott_L2(IA,JA)
816  real(RP), intent(in) :: qv_L2 (IA,JA)
817  real(RP), intent(in) :: qc_L2 (IA,JA)
818  real(RP), intent(in) :: dens_L1(IA,JA)
819  real(RP), intent(in) :: pott_L1(IA,JA)
820  real(RP), intent(in) :: qv_L1 (IA,JA)
821  real(RP), intent(in) :: qc_L1 (IA,JA)
822  real(RP), intent(in) :: dz (IA,JA)
823  integer, intent(in) :: k
824 
825  real(RP), intent(out) :: dens_L2(IA,JA)
826  real(RP), intent(out) :: temp_L2(IA,JA)
827  real(RP), intent(out) :: pres_L2(IA,JA)
828 
829  integer :: i, j
830  !---------------------------------------------------------------------------
831 
832  !$omp parallel do OMP_SCHEDULE_ collapse(2)
833  do j = js, je
834  do i = is, ie
835  call atmos_hydrostatic_buildrho_atmos_0d( &
836  pott_l2(i,j), qv_l2(i,j), qc_l2(i,j), & ! [IN]
837  dens_l1(i,j), pott_l1(i,j), qv_l1(i,j), qc_l1(i,j), & ! [IN]
838  dz(i,j), k, & ! [IN]
839  dens_l2(i,j), temp_l2(i,j), pres_l2(i,j) ) ! [OUT]
840  enddo
841  enddo
842 
843  return
844  end subroutine atmos_hydrostatic_buildrho_atmos_2d
845 
846  !-----------------------------------------------------------------------------
848  subroutine atmos_hydrostatic_buildrho_atmos_rev_2d( &
849  IA, IS, IE, JA, JS, JE, &
850  pott_L1, qv_L1, qc_L1, &
851  dens_L2, pott_L2, qv_L2, qc_L2, &
852  dz, k, &
853  dens_L1, temp_L1, pres_L1 )
854  implicit none
855  integer, intent(in) :: IA, IS, IE
856  integer, intent(in) :: JA, JS, JE
857 
858  real(RP), intent(in) :: pott_L1(IA,JA)
859  real(RP), intent(in) :: qv_L1 (IA,JA)
860  real(RP), intent(in) :: qc_L1 (IA,JA)
861  real(RP), intent(in) :: dens_L2(IA,JA)
862  real(RP), intent(in) :: pott_L2(IA,JA)
863  real(RP), intent(in) :: qv_L2 (IA,JA)
864  real(RP), intent(in) :: qc_L2 (IA,JA)
865  real(RP), intent(in) :: dz (IA,JA)
866  integer, intent(in) :: k
867 
868  real(RP), intent(out) :: dens_L1(IA,JA)
869  real(RP), intent(out) :: temp_L1(IA,JA)
870  real(RP), intent(out) :: pres_L1(IA,JA)
871 
872  integer :: i, j
873  !---------------------------------------------------------------------------
874 
875  !$omp parallel do OMP_SCHEDULE_ collapse(2)
876  do j = js, je
877  do i = is, ie
878  call atmos_hydrostatic_buildrho_atmos_0d( &
879  pott_l1(i,j), qv_l1(i,j), qc_l1(i,j), & ! [IN]
880  dens_l2(i,j), pott_l2(i,j), qv_l2(i,j), qc_l2(i,j), & ! [IN]
881  -dz(i,j), k, & ! [IN]
882  dens_l1(i,j), temp_l1(i,j), pres_l1(i,j) ) ! [OUT]
883  enddo
884  enddo
885 
886  return
887  end subroutine atmos_hydrostatic_buildrho_atmos_rev_2d
888 
889  !-----------------------------------------------------------------------------
891  subroutine atmos_hydrostatic_buildrho_atmos_3d( &
892  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
893  pott, qv, qc, &
894  dz, &
895  dens, &
896  temp, pres, &
897  kref )
898  implicit none
899  integer, intent(in) :: KA, KS, KE
900  integer, intent(in) :: IA, IS, IE
901  integer, intent(in) :: JA, JS, JE
902 
903  real(RP), intent(in) :: pott(KA,IA,JA)
904  real(RP), intent(in) :: qv (KA,IA,JA)
905  real(RP), intent(in) :: qc (KA,IA,JA)
906  real(RP), intent(in) :: dz (KA,IA,JA)
907 
908  real(RP), intent(inout) :: dens(KA,IA,JA)
909 
910  real(RP), intent(out) :: temp(KA,IA,JA)
911  real(RP), intent(out) :: pres(KA,IA,JA)
912 
913  integer, intent(in), optional, target :: kref(IA,JA)
914 
915  integer, pointer :: kref_(:,:)
916 
917  integer :: i, j
918  !---------------------------------------------------------------------------
919 
920  if ( present(kref) ) then
921  kref_ => kref
922  else
923  allocate( kref_(ia,ja) )
924  !$omp parallel do OMP_SCHEDULE_ collapse(2)
925  do j = js, je
926  do i = is, ie
927  kref_(i,j) = ks
928  end do
929  end do
930  end if
931 
932 
933  !$omp parallel do OMP_SCHEDULE_ collapse(2)
934  do j = js, je
935  do i = is, ie
936  call atmos_hydrostatic_buildrho_atmos_1d( ka, kref_(i,j), ke, &
937  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
938  dz(:,i,j), & ! [IN]
939  dens(:,i,j), & ! [INOUT]
940  temp(:,i,j), pres(:,i,j) ) ! [OUT]
941  enddo
942  enddo
943 
944  if ( .not. present(kref) ) then
945  deallocate( kref_ )
946  end if
947 
948  return
949  end subroutine atmos_hydrostatic_buildrho_atmos_3d
950 
951  !-----------------------------------------------------------------------------
953  subroutine atmos_hydrostatic_buildrho_atmos_rev_3d( &
954  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
955  pott, qv, qc, &
956  dz, &
957  dens, &
958  temp, pres, &
959  kref )
960  implicit none
961 
962  integer, intent(in) :: KA, KS, KE
963  integer, intent(in) :: IA, IS, IE
964  integer, intent(in) :: JA, JS, JE
965  real(RP), intent(inout) :: dens(KA,IA,JA)
966  real(RP), intent(out) :: temp(KA,IA,JA)
967  real(RP), intent(out) :: pres(KA,IA,JA)
968  real(RP), intent(in) :: pott(KA,IA,JA)
969  real(RP), intent(in) :: qv (KA,IA,JA)
970  real(RP), intent(in) :: qc (KA,IA,JA)
971  real(RP), intent(in) :: dz (KA,IA,JA)
972  integer, intent(in), optional, target :: kref(IA,JA)
973 
974  integer, pointer :: kref_(:,:)
975  integer :: i, j
976  !---------------------------------------------------------------------------
977 
978  if ( present(kref) ) then
979  kref_ => kref
980 
981  !$omp parallel do OMP_SCHEDULE_ collapse(2)
982  do j = js, je
983  do i = is, ie
984  call atmos_hydrostatic_buildrho_atmos_rev_1d( ka, ks, kref_(i,j), &
985  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
986  dz(:,i,j), & ! [IN]
987  dens(:,i,j), & ! [INOUT]
988  temp(:,i,j), pres(:,i,j) ) ! [OUT]
989  enddo
990  enddo
991  else
992  !$omp parallel do OMP_SCHEDULE_ collapse(2)
993  do j = js, je
994  do i = is, ie
995  call atmos_hydrostatic_buildrho_atmos_rev_1d( ka, ks, ke, &
996  pott(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
997  dz(:,i,j), & ! [IN]
998  dens(:,i,j), & ! [INOUT]
999  temp(:,i,j), pres(:,i,j) ) ! [OUT]
1000  enddo
1001  enddo
1002  end if
1003 
1004  return
1005  end subroutine atmos_hydrostatic_buildrho_atmos_rev_3d
1006 
1007  !-----------------------------------------------------------------------------
1009 !OCL SERIAL
1010  subroutine atmos_hydrostatic_buildrho_bytemp_1d( &
1011  KA, KS, KE, &
1012  temp, qv, qc, &
1013  pres_sfc, temp_sfc, qv_sfc, qc_sfc, &
1014  cz, fz, &
1015  dens, pott, pres, pott_sfc )
1016  use scale_prc, only: &
1017  prc_abort
1018  use scale_atmos_hydrometeor, only: &
1019  cv_vapor, &
1020  cv_water, &
1021  cp_vapor, &
1022  cp_water
1023  implicit none
1024  integer, intent(in) :: KA, KS, KE
1025 
1026  real(RP), intent(in) :: temp(KA)
1027  real(RP), intent(in) :: qv (KA)
1028  real(RP), intent(in) :: qc (KA)
1029  real(RP), intent(in) :: pres_sfc
1030  real(RP), intent(in) :: temp_sfc
1031  real(RP), intent(in) :: qv_sfc
1032  real(RP), intent(in) :: qc_sfc
1033  real(RP), intent(in) :: cz( KA)
1034  real(RP), intent(in) :: fz(0:KA)
1035 
1036  real(RP), intent(out) :: dens(KA)
1037  real(RP), intent(out) :: pott(KA)
1038  real(RP), intent(out) :: pres(KA)
1039  real(RP), intent(out) :: pott_sfc
1040 
1041  real(RP) :: dens_sfc
1042 
1043  real(RP) :: Rtot_sfc
1044  real(RP) :: CVtot_sfc
1045  real(RP) :: CPtot_sfc
1046  real(RP) :: Rtot
1047  real(RP) :: CVtot
1048  real(RP) :: CPtot
1049 
1050  real(RP) :: RovCP_sfc
1051  real(RP) :: dens_s, dhyd, dgrd
1052  integer :: ite
1053  logical :: converged
1054 
1055  real(RP) :: dz(KA)
1056 
1057  integer :: k
1058  !---------------------------------------------------------------------------
1059 
1060  !--- from surface to lowermost atmosphere
1061 
1062  rtot_sfc = rdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1063  + rvap * qv_sfc
1064  cvtot_sfc = cvdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1065  + cv_vapor * qv_sfc &
1066  + cv_water * qc_sfc
1067  cptot_sfc = cpdry * ( 1.0_rp - qv_sfc - qc_sfc ) &
1068  + cp_vapor * qv_sfc &
1069  + cp_water * qc_sfc
1070 
1071  rtot = rdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1072  + rvap * qv(ks)
1073  cvtot = cvdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1074  + cv_vapor * qv(ks) &
1075  + cv_water * qc(ks)
1076  cptot = cpdry * ( 1.0_rp - qv(ks) - qc(ks) ) &
1077  + cp_vapor * qv(ks) &
1078  + cp_water * qc(ks)
1079 
1080  ! density at surface
1081  rovcp_sfc = rtot_sfc / cptot_sfc
1082  dens_sfc = pres_sfc / ( rtot_sfc * temp_sfc )
1083  pott_sfc = temp_sfc * ( p00/pres_sfc )**rovcp_sfc
1084 
1085  ! make density at lowermost cell center
1086  dens_s = 0.0_rp
1087  dens(ks) = dens_sfc ! first guess
1088 
1089  dz(ks-1) = cz(ks) - fz(ks-1)
1090  do k = ks, ke-1
1091  dz(k) = cz(k+1) - cz(k)
1092  end do
1093 
1094  converged = .false.
1095  do ite = 1, itelim
1096  if ( abs(dens(ks)-dens_s) <= criteria ) then
1097  converged = .true.
1098  exit
1099  endif
1100 
1101  dens_s = dens(ks)
1102 
1103  dhyd = + ( pres_sfc - dens_s * rtot * temp(ks) ) / dz(ks-1) & ! dp/dz
1104  - grav * 0.5_rp * ( dens_sfc + dens_s ) ! rho*g
1105 
1106  dgrd = - rtot * temp(ks) / dz(ks-1) &
1107  - grav * 0.5_rp
1108 
1109  dens(ks) = dens_s - dhyd/dgrd
1110 
1111  if ( dens(ks)*0.0_rp /= 0.0_rp ) exit
1112  enddo
1113 
1114  if ( .NOT. converged ) then
1115  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_1D",*) 'iteration not converged!', &
1116  dens(ks),ite,dens_s,dhyd,dgrd
1117  call prc_abort
1118  endif
1119 
1120  !--- from lowermost atmosphere to top of atmosphere
1121  call atmos_hydrostatic_buildrho_bytemp_atmos_1d( ka, ks, ke, &
1122  temp(:), qv(:), qc(:), & ! [IN]
1123  dz(:), & ! [IN]
1124  dens(:), & ! [INOUT]
1125  pott(:), pres(:) ) ! [OUT]
1126 
1127  return
1128  end subroutine atmos_hydrostatic_buildrho_bytemp_1d
1129 
1130  !-----------------------------------------------------------------------------
1132  subroutine atmos_hydrostatic_buildrho_bytemp_3d( &
1133  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1134  temp, qv, qc, &
1135  pres_sfc, temp_sfc, qv_sfc, qc_sfc, &
1136  cz, fz, &
1137  dens, pott, pres, &
1138  pott_sfc )
1139  implicit none
1140  integer, intent(in) :: KA, KS, KE
1141  integer, intent(in) :: IA, IS, IE
1142  integer, intent(in) :: JA, JS, JE
1143 
1144  real(RP), intent(in) :: temp(KA,IA,JA)
1145  real(RP), intent(in) :: qv (KA,IA,JA)
1146  real(RP), intent(in) :: qc (KA,IA,JA)
1147  real(RP), intent(in) :: pres_sfc(IA,JA)
1148  real(RP), intent(in) :: temp_sfc(IA,JA)
1149  real(RP), intent(in) :: qv_sfc (IA,JA)
1150  real(RP), intent(in) :: qc_sfc (IA,JA)
1151  real(RP), intent(in) :: cz(KA,IA,JA)
1152  real(RP), intent(in) :: fz(KA,IA,JA)
1153 
1154  real(RP), intent(out) :: dens(KA,IA,JA)
1155  real(RP), intent(out) :: pott(KA,IA,JA)
1156  real(RP), intent(out) :: pres(KA,IA,JA)
1157  real(RP), intent(out) :: pott_sfc(IA,JA)
1158 
1159  integer :: i, j
1160  !---------------------------------------------------------------------------
1161 
1162  !--- from surface to lowermost atmosphere
1163  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1164  do j = js, je
1165  do i = is, ie
1166  call atmos_hydrostatic_buildrho_bytemp_1d( ka, ks, ke, &
1167  temp(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
1168  pres_sfc(i,j), temp_sfc(i,j), qv_sfc(i,j), qc_sfc(i,j), & ! [IN]
1169  cz(:,i,j), fz(:,i,j), & ! [IN]
1170  dens(:,i,j), pott(:,i,j), pres(:,i,j), & ! [OUT]
1171  pott_sfc(i,j) ) ! [OUT]
1172  enddo
1173  enddo
1174 
1175  return
1176  end subroutine atmos_hydrostatic_buildrho_bytemp_3d
1177 
1178  !-----------------------------------------------------------------------------
1180 !OCL SERIAL
1181  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_1d( &
1182  KA, KS, KE, &
1183  temp, qv, qc, &
1184  dz, &
1185  dens, &
1186  pott, pres )
1187  use scale_prc, only: &
1188  prc_abort
1189  use scale_atmos_hydrometeor, only: &
1190  cv_vapor, &
1191  cv_water, &
1192  cp_vapor, &
1193  cp_water
1194  implicit none
1195  integer, intent(in) :: KA, KS, KE
1196 
1197  real(RP), intent(in) :: temp(KA)
1198  real(RP), intent(in) :: qv (KA)
1199  real(RP), intent(in) :: qc (KA)
1200  real(RP), intent(in) :: dz (KA)
1201 
1202  real(RP), intent(inout) :: dens(KA)
1203 
1204  real(RP), intent(out) :: pott(KA)
1205  real(RP), intent(out) :: pres(KA)
1206 
1207  real(RP) :: Rtot (KA)
1208  real(RP) :: CVtot (KA)
1209  real(RP) :: CPtot (KA)
1210 
1211  real(RP) :: RovCP
1212  real(RP) :: dens_s, dhyd, dgrd
1213  integer :: ite
1214  logical :: converged
1215 
1216  integer :: k
1217  !---------------------------------------------------------------------------
1218 
1219  do k = ks, ke
1220  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
1221  + rvap * qv(k)
1222  cvtot(k) = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
1223  + cv_vapor * qv(k) &
1224  + cv_water * qc(k)
1225  cptot(k) = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
1226  + cp_vapor * qv(k) &
1227  + cp_water * qc(k)
1228  enddo
1229 
1230  pres(ks) = dens(ks) * rtot(ks) * temp(ks)
1231 
1232  do k = ks+1, ke
1233 
1234  dens_s = 0.0_rp
1235  dens(k) = dens(k-1) ! first guess
1236 
1237  converged = .false.
1238  do ite = 1, itelim
1239  if ( abs(dens(k)-dens_s) <= criteria ) then
1240  converged = .true.
1241  exit
1242  endif
1243 
1244  dens_s = dens(k)
1245 
1246  dhyd = + ( pres(k-1) - dens_s * rtot(k) * temp(k) ) / dz(k-1) & ! dp/dz
1247  - grav * 0.5_rp * ( dens(k-1) + dens_s ) ! rho*g
1248 
1249  dgrd = - rtot(k) * temp(k) / dz(k-1) &
1250  - 0.5_rp * grav
1251 
1252  dens(k) = dens_s - dhyd/dgrd
1253 
1254  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
1255  enddo
1256 
1257  if ( .NOT. converged ) then
1258  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_atmos_1D",*) 'iteration not converged!', &
1259  k,dens(k),ite,dens_s,dhyd,dgrd
1260  call prc_abort
1261  endif
1262 
1263  pres(k) = dens(k) * rtot(k) * temp(k)
1264 
1265  enddo
1266 
1267  do k = ks, ke
1268  rovcp = rtot(k) / cptot(k)
1269  pott(k) = temp(k) * ( p00 / pres(k) )**rovcp
1270  enddo
1271 
1272  return
1273  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_1d
1274 
1275  !-----------------------------------------------------------------------------
1277 !OCL SERIAL
1278  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d( &
1279  KA, KS, KE, &
1280  temp, qv, qc, &
1281  dz, &
1282  dens, &
1283  pott, pres )
1284  use scale_prc, only: &
1285  prc_abort
1286  use scale_atmos_hydrometeor, only: &
1287  cv_vapor, &
1288  cv_water, &
1289  cp_vapor, &
1290  cp_water
1291  implicit none
1292  integer, intent(in) :: KA, KS, KE
1293 
1294  real(RP), intent(in) :: temp(KA)
1295  real(RP), intent(in) :: qv (KA)
1296  real(RP), intent(in) :: qc (KA)
1297  real(RP), intent(in) :: dz (KA)
1298 
1299  real(RP), intent(inout) :: dens(KA)
1300 
1301  real(RP), intent(out) :: pott(KA)
1302  real(RP), intent(out) :: pres(KA)
1303 
1304  real(RP) :: Rtot (KA)
1305  real(RP) :: CVtot (KA)
1306  real(RP) :: CPtot (KA)
1307 
1308  real(RP) :: RovCP
1309  real(RP) :: dens_s, dhyd, dgrd
1310  integer :: ite
1311  logical :: converged
1312 
1313  integer :: k
1314  !---------------------------------------------------------------------------
1315 
1316  do k = ks, ke
1317  rtot(k) = rdry * ( 1.0_rp - qv(k) - qc(k) ) &
1318  + rvap * qv(k)
1319  cvtot(k) = cvdry * ( 1.0_rp - qv(k) - qc(k) ) &
1320  + cv_vapor * qv(k) &
1321  + cv_water * qc(k)
1322  cptot(k) = cpdry * ( 1.0_rp - qv(k) - qc(k) ) &
1323  + cp_vapor * qv(k) &
1324  + cp_water * qc(k)
1325  enddo
1326 
1327  pres(ke) = dens(ke) * rtot(ke) * temp(ke)
1328 
1329  do k = ke-1, ke, -1
1330 
1331  dens_s = 0.0_rp
1332  dens(k) = dens(k+1) ! first guess
1333 
1334  converged = .false.
1335  do ite = 1, itelim
1336  if ( abs(dens(k)-dens_s) <= criteria ) then
1337  converged = .true.
1338  exit
1339  endif
1340 
1341  dens_s = dens(k)
1342 
1343  dhyd = - ( pres(k+1) - dens_s * rtot(k) * temp(k) ) / dz(k) & ! dp/dz
1344  - grav * 0.5_rp * ( dens(k+1) + dens_s ) ! rho*g
1345 
1346  dgrd = + rtot(k) * temp(k) / dz(k) &
1347  - 0.5_rp * grav
1348 
1349  dens(k) = dens_s - dhyd/dgrd
1350 
1351  if ( dens(k)*0.0_rp /= 0.0_rp ) exit
1352  enddo
1353 
1354  if ( .NOT. converged ) then
1355  log_error("ATMOS_HYDROSTATIC_buildrho_bytemp_atmos_rev_1D",*) 'iteration not converged!', &
1356  k,dens(k),ite,dens_s,dhyd,dgrd
1357  call prc_abort
1358  endif
1359 
1360  pres(k) = dens(k) * rtot(k) * temp(k)
1361 
1362  enddo
1363 
1364  do k = ks, ke
1365  rovcp = rtot(k) / cptot(k)
1366  pott(k) = temp(k) * ( p00 / pres(k) )**rovcp
1367  enddo
1368 
1369  return
1370  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_rev_1d
1371 
1372  !-----------------------------------------------------------------------------
1374  subroutine atmos_hydrostatic_buildrho_bytemp_atmos_3d( &
1375  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1376  temp, qv, qc, &
1377  dz, &
1378  dens, &
1379  pott, pres )
1380  implicit none
1381  integer, intent(in) :: KA, KS, KE
1382  integer, intent(in) :: IA, IS, IE
1383  integer, intent(in) :: JA, JS, JE
1384 
1385  real(RP), intent(in) :: temp(KA,IA,JA)
1386  real(RP), intent(in) :: qv (KA,IA,JA)
1387  real(RP), intent(in) :: qc (KA,IA,JA)
1388  real(RP), intent(in) :: dz (KA,IA,JA)
1389 
1390  real(RP), intent(inout) :: dens(KA,IA,JA)
1391  real(RP), intent(out) :: pott(KA,IA,JA)
1392  real(RP), intent(out) :: pres(KA,IA,JA)
1393 
1394  integer :: i, j
1395  !---------------------------------------------------------------------------
1396 
1397  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1398  do j = js, je
1399  do i = is, ie
1400  call atmos_hydrostatic_buildrho_bytemp_atmos_1d( ka, ks, ke, &
1401  temp(:,i,j), qv(:,i,j), qc(:,i,j), & ! [IN]
1402  dz(:,i,j), & ! [IN]
1403  dens(:,i,j), & ! [INOUT]
1404  pott(:,i,j), pres(:,i,j) ) ! [OUT]
1405  enddo
1406  enddo
1407 
1408  return
1409  end subroutine atmos_hydrostatic_buildrho_bytemp_atmos_3d
1410 
1411  !-----------------------------------------------------------------------------
1413 !OCL SERIAL
1414 !OCL NOSIMD
1415  subroutine atmos_hydrostatic_barometric_law_mslp_0d( &
1416  KA, KS, KE, &
1417  pres, temp, qv, &
1418  cz, &
1419  mslp )
1420  implicit none
1421  integer, intent(in) :: KA, KS, KE
1422 
1423  real(RP), intent(in) :: pres(KA)
1424  real(RP), intent(in) :: temp(KA)
1425  real(RP), intent(in) :: qv (KA)
1426  real(RP), intent(in) :: cz (KA)
1427 
1428  real(RP), intent(out) :: mslp
1429 
1430  ! work
1431  integer :: kref
1432  real(RP) :: vtemp
1433  !---------------------------------------------------------------------------
1434 
1435  kref = hydrostatic_barometric_law_mslp_kref + ks - 1
1436 
1437  ! virtual temperature
1438  vtemp = temp(kref) * (1.0_rp + epstvap * qv(kref) )
1439 
1440  ! barometric law assuming constant lapse rate
1441  mslp = pres(kref) * ( ( vtemp + laps * cz(kref) ) / vtemp ) ** ( grav / ( rdry * laps ) )
1442 
1443  return
1444  end subroutine atmos_hydrostatic_barometric_law_mslp_0d
1445 
1446  !-----------------------------------------------------------------------------
1448  subroutine atmos_hydrostatic_barometric_law_mslp_2d( &
1449  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1450  pres, temp, qv, &
1451  cz, &
1452  mslp )
1453  implicit none
1454  integer, intent(in) :: KA, KS, KE
1455  integer, intent(in) :: IA, IS, IE
1456  integer, intent(in) :: JA, JS, JE
1457 
1458  real(RP), intent(in) :: pres(KA,IA,JA)
1459  real(RP), intent(in) :: temp(KA,IA,JA)
1460  real(RP), intent(in) :: qv (KA,IA,JA)
1461  real(RP), intent(in) :: cz (KA,IA,JA)
1462 
1463  real(RP), intent(out) :: mslp(IA,JA)
1464 
1465  integer :: i, j
1466  !---------------------------------------------------------------------------
1467 
1468  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1469  do j = js, je
1470  do i = is, ie
1471  call atmos_hydrostatic_barometric_law_mslp_0d( ka, ks, ke, &
1472  pres(:,i,j), temp(:,i,j), qv(:,i,j), & ! [IN]
1473  cz(:,i,j), & ! [IN]
1474  mslp(i,j) ) ! [OUT]
1475  enddo
1476  enddo
1477 
1478  return
1479  end subroutine atmos_hydrostatic_barometric_law_mslp_2d
1480 
1481  !-----------------------------------------------------------------------------
1483 !OCL SERIAL
1484 !OCL NOSIMD
1485  subroutine atmos_hydrostatic_barometric_law_pres_0d( &
1486  mslp, temp, &
1487  dz, &
1488  pres )
1489  implicit none
1490  real(RP), intent(in) :: mslp
1491  real(RP), intent(in) :: temp
1492  real(RP), intent(in) :: dz
1493 
1494  real(RP), intent(out) :: pres
1495 
1496  ! work
1497  real(RP) :: TM
1498  !---------------------------------------------------------------------------
1499 
1500  tm = temp + laps * dz * 0.5_rp ! column-mean air temperature
1501 
1502  ! barometric law
1503  pres = mslp / exp( grav * dz / ( rdry * tm ) )
1504 
1505  return
1506  end subroutine atmos_hydrostatic_barometric_law_pres_0d
1507 
1508  !-----------------------------------------------------------------------------
1510  subroutine atmos_hydrostatic_barometric_law_pres_2d( &
1511  IA, IS, IE, JA, JS, JE, &
1512  mslp, temp, &
1513  dz, &
1514  pres )
1515  implicit none
1516  integer, intent(in) :: IA, IS, IE
1517  integer, intent(in) :: JA, JS, JE
1518 
1519  real(RP), intent(in) :: mslp(IA,JA)
1520  real(RP), intent(in) :: temp(IA,JA)
1521  real(RP), intent(in) :: dz (IA,JA)
1522 
1523  real(RP), intent(out) :: pres(IA,JA)
1524 
1525  integer :: i, j
1526  !---------------------------------------------------------------------------
1527 
1528  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1529  do j = js, je
1530  do i = is, ie
1531  call atmos_hydrostatic_barometric_law_pres_0d( mslp(i,j), temp(i,j), dz(i,j), & ! [IN]
1532  pres(i,j) ) ! [OUT]
1533  enddo
1534  enddo
1535 
1536  return
1537  end subroutine atmos_hydrostatic_barometric_law_pres_2d
1538 
1539 end module scale_atmos_hydrostatic
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_const::const_lapsdry
real(rp), public const_lapsdry
dry adiabatic lapse rate [K/m]
Definition: scale_const.F90:59
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:133
scale_precision
module PRECISION
Definition: scale_precision.F90:14
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:1284
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
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_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
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_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:130
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const
module CONSTANT
Definition: scale_const.F90:11
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:1380
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:57
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:1139
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:504
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_prof
module profiler
Definition: scale_prof.F90:11
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)
Build up density from upermost atmosphere (1D)
Definition: scale_atmos_hydrostatic.F90:708
scale_atmos_hydrostatic::atmos_hydrostatic_buildrho_atmos_1d
subroutine atmos_hydrostatic_buildrho_atmos_1d(KA, KS, KE, pott, qv, qc, dz, dens, temp, pres)
Build up density from lowermost atmosphere (1D)
Definition: scale_atmos_hydrostatic.F90:604
scale_atmos_hydrostatic
module atmosphere / hydrostatic barance
Definition: scale_atmos_hydrostatic.F90:12
scale_const::const_cvvap
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
Definition: scale_const.F90:65
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_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_atmos_hydrostatic::atmos_hydrostatic_setup
subroutine, public atmos_hydrostatic_setup
Setup.
Definition: scale_atmos_hydrostatic.F90:116
scale_const::const_laps
real(rp), public const_laps
lapse rate of ISA [K/m]
Definition: scale_const.F90:58
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
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_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
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::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:132
Generated by   doxygen 1.8.17