SCALE-RM
scale_atmos_adiabat.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
11 module scale_atmos_adiabat
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
19  !-----------------------------------------------------------------------------
20  implicit none
21  private
22  !-----------------------------------------------------------------------------
23  !
24  !++ Public procedure
25  !
26  public :: atmos_adiabat_setup
27  public :: atmos_adiabat_cape
28  public :: atmos_adiabat_liftparcel
29 
30  interface atmos_adiabat_cape
31  module procedure atmos_adiabat_cape_1d
32  module procedure atmos_adiabat_cape_3d
33  end interface atmos_adiabat_cape
34 
35  interface atmos_adiabat_liftparcel
36  module procedure atmos_adiabat_liftparcel_1d
37  module procedure atmos_adiabat_liftparcel_3d
38  end interface atmos_adiabat_liftparcel
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  !-----------------------------------------------------------------------------
53 contains
54  !-----------------------------------------------------------------------------
56  !-----------------------------------------------------------------------------
57  subroutine atmos_adiabat_setup
58  use scale_atmos_hydrometeor, only: &
59  atmos_hydrometeor_setup
60  use scale_atmos_saturation, only: &
61  atmos_saturation_setup
62 
63  call atmos_hydrometeor_setup
64  call atmos_saturation_setup
65 
66  return
67  end subroutine atmos_adiabat_setup
68 
69  !-----------------------------------------------------------------------------
73 !OCL SERIAL
74  subroutine atmos_adiabat_cape_1d( &
75  KA, KS, KE, &
76  Kstr, &
77  DENS, TEMP, PRES, &
78  QV, QC, Qdry, Rtot, CPtot, &
79  CZ, FZ, &
80  CAPE, CIN, LCL, LFC, LNB, &
81  DENS_p, TEMP_p, BUOY_p, QV_p, &
82  converged )
83  use scale_const, only: &
84  grav => const_grav
85  implicit none
86  integer, intent(in) :: KA, KS, KE
87 
88  integer, intent(in) :: Kstr
89  real(RP), intent(in) :: DENS (KA)
90  real(RP), intent(in) :: TEMP (KA)
91  real(RP), intent(in) :: PRES (KA)
92  real(RP), intent(in) :: QV (KA)
93  real(RP), intent(in) :: QC (KA)
94  real(RP), intent(in) :: Qdry (KA)
95  real(RP), intent(in) :: Rtot (KA)
96  real(RP), intent(in) :: CPtot(KA)
97  real(RP), intent(in) :: CZ ( KA)
98  real(RP), intent(in) :: FZ (0:KA)
99 
100  real(RP), intent(out) :: CAPE
101  real(RP), intent(out) :: CIN
102  real(RP), intent(out) :: LCL
103  real(RP), intent(out) :: LFC
104  real(RP), intent(out) :: LNB
105  real(RP), intent(out) :: DENS_p(KA)
106  real(RP), intent(out) :: TEMP_p(KA)
107  real(RP), intent(out) :: BUOY_p(KA)
108  real(RP), intent(out) :: QV_p (KA)
109  logical, intent(out) :: converged
110 
111  real(RP) :: BUOY_pf(KA)
112  integer :: kLCL, kLFC, kLNB
113 
114  integer :: k
115  !---------------------------------------------------------------------------
116 
117  ! lift parcel
118  call atmos_adiabat_liftparcel_1d( ka, ks, ke, &
119  kstr, & ! [IN]
120  dens(:), temp(:), pres(:), & ! [IN]
121  qv(:), qc(:), & ! [IN]
122  qdry(:), rtot(:), cptot(:), & ! [IN]
123  dens_p(:), temp_p(:), qv_p(:), & ! [OUT]
124  klcl, converged ) ! [OUT]
125 
126  if ( .not. converged ) return
127 
128  ! parcel buoyancy
129  do k = ks, ke
130  buoy_p(k) = grav * ( dens(k) - dens_p(k) ) / dens_p(k)
131  end do
132 
133  do k = ks, ke-1
134  buoy_pf(k) = 0.5_rp * ( buoy_p(k+1) + buoy_p(k) )
135  end do
136  buoy_p(ke) = buoy_p(ke-1)
137 
138  lcl = 0.0_rp
139  lfc = 0.0_rp
140  lnb = 0.0_rp
141  cape = 0.0_rp
142  cin = 0.0_rp
143 
144  klfc = -1
145  klnb = -1
146 
147  do k = klcl, ke
148  if ( buoy_p(k) >= 0.0_rp ) then
149  klfc = k
150  exit
151  endif
152  enddo
153 
154  do k = ke, kstr, -1
155  if ( buoy_p(k) >= 0.0_rp ) then
156  klnb = k
157  exit
158  endif
159  enddo
160 
161  if ( klcl >= kstr ) then
162  lcl = cz(klcl)
163  endif
164 
165  if ( klfc >= kstr ) then
166  lfc = cz(klfc)
167  endif
168 
169  if ( klnb >= kstr ) then
170  lnb = cz(klnb)
171  endif
172 
173  if ( klfc >= kstr .AND. klnb > kstr ) then
174  do k = klfc+1, klnb
175  cape = cape + buoy_pf(k-1) * ( fz(k)-fz(k-1) )
176  enddo
177  endif
178 
179  if ( klfc >= kstr ) then
180  do k = kstr+1, klfc
181  cin = cin + buoy_pf(k-1) * ( fz(k)-fz(k-1) )
182  enddo
183  endif
184 
185  return
186  end subroutine atmos_adiabat_cape_1d
187 
188  subroutine atmos_adiabat_cape_3d( &
189  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
190  Kstr, &
191  DENS, TEMP, PRES, &
192  QV, QC, Qdry, Rtot, CPtot, &
193  CZ, FZ, &
194  CAPE, CIN, LCL, LFC, LNB, &
195  DENS_p, TEMP_p, BUOY_p, QV_p )
196  use scale_prc, only: &
197  prc_abort
198  use scale_file_history, only: &
199  file_history_in
200  implicit none
201  integer, intent(in) :: KA, KS, KE
202  integer, intent(in) :: IA, IS, IE
203  integer, intent(in) :: JA, JS, JE
204 
205  integer, intent(in) :: Kstr
206  real(RP), intent(in) :: DENS (KA,IA,JA)
207  real(RP), intent(in) :: TEMP (KA,IA,JA)
208  real(RP), intent(in) :: PRES (KA,IA,JA)
209  real(RP), intent(in) :: QV (KA,IA,JA)
210  real(RP), intent(in) :: QC (KA,IA,JA)
211  real(RP), intent(in) :: Qdry (KA,IA,JA)
212  real(RP), intent(in) :: Rtot (KA,IA,JA)
213  real(RP), intent(in) :: CPtot(KA,IA,JA)
214  real(RP), intent(in) :: CZ ( KA,IA,JA)
215  real(RP), intent(in) :: FZ (0:KA,IA,JA)
216 
217  real(RP), intent(out) :: CAPE(IA,JA)
218  real(RP), intent(out) :: CIN (IA,JA)
219  real(RP), intent(out) :: LCL (IA,JA)
220  real(RP), intent(out) :: LFC (IA,JA)
221  real(RP), intent(out) :: LNB (IA,JA)
222 
223  real(RP), intent(out), optional, target :: DENS_p(KA,IA,JA)
224  real(RP), intent(out), optional, target :: TEMP_p(KA,IA,JA)
225  real(RP), intent(out), optional, target :: BUOY_p(KA,IA,JA)
226  real(RP), intent(out), optional, target :: QV_p (KA,IA,JA)
227 
228  logical :: converged, error
229 
230  real(RP), pointer :: P_DENS(:,:,:), P_TEMP(:,:,:), P_BUOY(:,:,:), P_QV(:,:,:)
231 
232  integer :: i, j
233  !---------------------------------------------------------------------------
234 
235  error = .false.
236 
237  if ( present(dens_p) ) then
238  p_dens => dens_p
239  else
240  allocate( p_dens(ka,ia,ja) )
241  end if
242  if ( present(temp_p) ) then
243  p_temp => temp_p
244  else
245  allocate( p_temp(ka,ia,ja) )
246  end if
247  if ( present(buoy_p) ) then
248  p_buoy => buoy_p
249  else
250  allocate( p_buoy(ka,ia,ja) )
251  end if
252  if ( present(qv_p) ) then
253  p_qv => qv_p
254  else
255  allocate( p_qv(ka,ia,ja) )
256  end if
257 
258 
259  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
260  !$omp private(converged)
261  do j = js, je
262  do i = is, ie
263 
264  call atmos_adiabat_cape_1d( ka, ks, ke, &
265  kstr, &
266  dens(:,i,j), temp(:,i,j), pres(:,i,j), & ! [IN]
267  qv(:,i,j), qc(:,i,j), qdry(:,i,j), & ! [IN]
268  rtot(:,i,j), cptot(:,i,j), & ! [IN]
269  cz(:,i,j), fz(:,i,j), & ! [IN]
270  cape(i,j), cin(i,j), & ! [OUT]
271  lcl(i,j), lfc(i,j), lnb(i,j), & ! [OUT]
272  p_dens(:,i,j), p_temp(:,i,j), & ! [OUT]
273  p_buoy(:,i,j), p_qv(:,i,j), & ! [OUT]
274  converged ) ! [OUT]
275 
276  if ( .not. converged ) then
277  log_error("ATMOS_ADIABAT_cape_3D",*) '[liftparcel] not converged! ', i, j
278  error = .true.
279  end if
280 
281  end do
282  end do
283 
284  if ( error ) call prc_abort
285 
286  call file_history_in( p_dens(:,:,:), 'DENS_parcel', 'density profile in lifting parcel', 'kg/m3' )
287  call file_history_in( p_temp(:,:,:), 'TEMP_parcel', 'temperature profile in lifting parcel', 'K' )
288  call file_history_in( p_buoy(:,:,:), 'BUOY_parcel', 'buoyancy profile in lifting parcel', 'm/s2' )
289  call file_history_in( p_qv(:,:,:), 'QV_parcel', 'humidity profile in lifting parcel', 'kg/kg' )
290 
291  if ( .not. present(dens_p) ) deallocate( p_dens )
292  if ( .not. present(temp_p) ) deallocate( p_temp )
293  if ( .not. present(buoy_p) ) deallocate( p_buoy )
294  if ( .not. present(qv_p ) ) deallocate( p_qv )
295 
296  return
297  end subroutine atmos_adiabat_cape_3d
298 
299  !-----------------------------------------------------------------------------
303 !OCL SERIAL
304  subroutine atmos_adiabat_liftparcel_1d( &
305  KA, KS, KE, &
306  Kstr, &
307  DENS, TEMP, PRES, QV, QC, &
308  QDRY, Rtot, CPtot, &
309  DENS_p1D, TEMP_p1D, QV_p1D, &
310  kLCL, converged )
311  use scale_atmos_hydrometeor, only: &
312  hydrometeor_entr => atmos_hydrometeor_entr, &
313  cp_water
314  use scale_atmos_saturation, only: &
315  atmos_saturation_moist_conversion_pres_liq
316  implicit none
317  integer, intent(in) :: KA, KS, KE
318 
319  integer, intent(in) :: Kstr
320  real(RP), intent(in) :: DENS (KA)
321  real(RP), intent(in) :: TEMP (KA)
322  real(RP), intent(in) :: PRES (KA)
323  real(RP), intent(in) :: QV (KA)
324  real(RP), intent(in) :: QC (KA)
325  real(RP), intent(in) :: QDRY (KA)
326  real(RP), intent(in) :: Rtot (KA)
327  real(RP), intent(in) :: CPtot(KA)
328 
329  real(RP), intent(out) :: DENS_p1D(KA)
330  real(RP), intent(out) :: TEMP_p1D(KA)
331  real(RP), intent(out) :: QV_p1D (KA)
332  integer, intent(out) :: kLCL
333  logical, intent(out) :: converged
334 
335  real(RP) :: ENTR_p
336  real(RP) :: QV_p, QC_p, Qdry_p
337  real(RP) :: Rtot_p, CPtot_p
338  real(RP) :: TEMP_p
339 
340  real(RP) :: fact
341 
342  integer :: k
343  !---------------------------------------------------------------------------
344 
345  do k = ks, kstr
346  dens_p1d(k) = dens(k)
347  temp_p1d(k) = temp(k)
348  qv_p1d(k) = qv(k)
349  enddo
350 
351  ! vapor + liquid water at the start point
352  qv_p = qv(kstr)
353  qc_p = qc(kstr)
354  qdry_p = qdry(kstr)
355  rtot_p = rtot(kstr)
356  cptot_p = cptot(kstr)
357 
358  call hydrometeor_entr( temp(kstr), pres(kstr), & ! [IN]
359  qv_p, qc_p, qdry_p, & ! [IN]
360  rtot_p, cptot_p, & ! [IN]
361  entr_p ) ! [OUT]
362 
363  klcl = -1
364  do k = kstr, ke
365 
366  call atmos_saturation_moist_conversion_pres_liq( &
367  pres(k), entr_p, qdry_p, & ! [IN]
368  qv_p, qc_p, rtot_p, cptot_p, & ! [INOUT]
369  temp_p, & ! [OUT]
370  converged ) ! [OUT]
371 
372  if ( .NOT. converged ) then
373  exit
374  endif
375 
376  if ( qc_p > 0.0_rp .and. klcl == -1 ) klcl = k
377 
378  ! remove condensed water
379  fact = 1.0_rp / ( 1.0_rp - qc_p )
380  cptot_p = ( cptot_p - cp_water * qc_p ) * fact
381  rtot_p = rtot_p * fact
382  qdry_p = qdry_p * fact
383  qv_p = qv_p * fact
384  ! Entr_p = Entr_p - QC_p * CP_WATER * log( TEMP_p / TEM00 )
385  qc_p = 0.0_rp
386 
387  dens_p1d(k) = pres(k) / ( rtot_p * temp_p )
388  temp_p1d(k) = temp_p
389  qv_p1d(k) = qv_p
390 
391  enddo
392 
393  return
394  end subroutine atmos_adiabat_liftparcel_1d
395 
396  subroutine atmos_adiabat_liftparcel_3d( &
397  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
398  Kstr, &
399  DENS, TEMP, PRES, QV, QC, &
400  QDRY, Rtot, CPtot, &
401  DENS_p3D, TEMP_p3D, QV_p3D, &
402  kLCL, converged )
403  use scale_prc, only: &
404  prc_abort
405  use scale_atmos_hydrometeor, only: &
406  hydrometeor_entr => atmos_hydrometeor_entr
407  use scale_atmos_saturation, only: &
408  atmos_saturation_moist_conversion_pres_liq
409  implicit none
410 
411  integer, intent(in) :: KA, KS, KE
412  integer, intent(in) :: IA, IS, IE
413  integer, intent(in) :: JA, JS, JE
414 
415  integer, intent(in) :: Kstr
416  real(RP), intent(in) :: DENS (KA,IA,JA)
417  real(RP), intent(in) :: TEMP (KA,IA,JA)
418  real(RP), intent(in) :: PRES (KA,IA,JA)
419  real(RP), intent(in) :: QV (KA,IA,JA)
420  real(RP), intent(in) :: QC (KA,IA,JA)
421  real(RP), intent(in) :: QDRY (KA,IA,JA)
422  real(RP), intent(in) :: Rtot (KA,IA,JA)
423  real(RP), intent(in) :: CPtot(KA,IA,JA)
424 
425  real(RP), intent(out) :: DENS_p3D(KA,IA,JA)
426  real(RP), intent(out) :: TEMP_p3D(KA,IA,JA)
427  real(RP), intent(out) :: QV_p3D (KA,IA,JA)
428  integer, intent(out) :: kLCL(IA,JA)
429  logical, intent(out) :: converged
430 
431  logical :: error
432 
433  integer :: i, j
434  !---------------------------------------------------------------------------
435 
436  error = .false.
437 
438  !$omp parallel do OMP_SCHEDULE_ collapse(2)
439  do j = js, je
440  do i = is, ie
441 
442  call atmos_adiabat_liftparcel_1d( ka, ks, ke, &
443  kstr, & ! [IN]
444  dens(:,i,j), temp(:,i,j), pres(:,i,j), & ! [IN]
445  qv(:,i,j), qc(:,i,j), & ! [IN]
446  qdry(:,i,j), rtot(:,i,j), cptot(:,i,j), & ! [IN]
447  dens_p3d(:,i,j), temp_p3d(:,i,j), qv_p3d(:,i,j), & ! [OUT]
448  klcl(i,j), converged ) ! [OUT]
449  if ( .not. converged ) then
450  log_error("ATMOS_ADIABAT_liftparcel_3D",*) '[liftparcel] not converged! ', i, j
451  error = .true.
452  end if
453 
454  enddo
455  enddo
456 
457  if ( error ) call prc_abort
458 
459  return
460  end subroutine atmos_adiabat_liftparcel_3d
461 
462 end module scale_atmos_adiabat
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
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_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_adiabat::atmos_adiabat_cape_1d
subroutine atmos_adiabat_cape_1d(KA, KS, KE, Kstr, DENS, TEMP, PRES, QV, QC, Qdry, Rtot, CPtot, CZ, FZ, CAPE, CIN, LCL, LFC, LNB, DENS_p, TEMP_p, BUOY_p, QV_p, converged)
Calc CAPE and CIN Type of parcel method: Pseudo-adiabatic ascend from lowermost layer of the model Re...
Definition: scale_atmos_adiabat.F90:83
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_adiabat
module atmosphere / adiabat
Definition: scale_atmos_adiabat.F90:11
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_adiabat::atmos_adiabat_liftparcel_3d
subroutine atmos_adiabat_liftparcel_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, Kstr, DENS, TEMP, PRES, QV, QC, QDRY, Rtot, CPtot, DENS_p3D, TEMP_p3D, QV_p3D, kLCL, converged)
Definition: scale_atmos_adiabat.F90:403
scale_atmos_adiabat::atmos_adiabat_liftparcel_1d
subroutine atmos_adiabat_liftparcel_1d(KA, KS, KE, Kstr, DENS, TEMP, PRES, QV, QC, QDRY, Rtot, CPtot, DENS_p1D, TEMP_p1D, QV_p1D, kLCL, converged)
Calc temperature profile with lifting parcel Method: Pseudo-adiabatic ascend from lowermost layer of ...
Definition: scale_atmos_adiabat.F90:311
scale_atmos_saturation::atmos_saturation_setup
subroutine, public atmos_saturation_setup
Setup.
Definition: scale_atmos_saturation.F90:215
scale_atmos_adiabat::atmos_adiabat_setup
subroutine, public atmos_adiabat_setup
Setup.
Definition: scale_atmos_adiabat.F90:58
scale_atmos_hydrometeor::atmos_hydrometeor_setup
subroutine, public atmos_hydrometeor_setup
Setup.
Definition: scale_atmos_hydrometeor.F90:154
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_adiabat::atmos_adiabat_cape_3d
subroutine atmos_adiabat_cape_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, Kstr, DENS, TEMP, PRES, QV, QC, Qdry, Rtot, CPtot, CZ, FZ, CAPE, CIN, LCL, LFC, LNB, DENS_p, TEMP_p, BUOY_p, QV_p)
Definition: scale_atmos_adiabat.F90:196
Generated by   doxygen 1.8.17