SCALE-RM
scale_atmos_sub_adiabat.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 module scale_atmos_adiabat
11  !-----------------------------------------------------------------------------
12  !
13  !++ used modules
14  !
15  use scale_precision
16  use scale_stdio
17  use scale_prof
18  use scale_grid_index
19  use scale_tracer
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: atmos_adiabat_cape
28  public :: atmos_adiabat_liftparcel
29 
30  !-----------------------------------------------------------------------------
31  !
32  !++ Public parameters & variables
33  !
34  !-----------------------------------------------------------------------------
35  !
36  !++ Private procedure
37  !
38  !-----------------------------------------------------------------------------
39  !
40  !++ Private parameters & variables
41  !
42  !-----------------------------------------------------------------------------
43 contains
44  !-----------------------------------------------------------------------------
48  subroutine atmos_adiabat_cape( &
49  Kstr, &
50  DENS, &
51  TEMP, &
52  PRES, &
53  QTRC, &
54  CZ, &
55  FZ, &
56  CAPE, &
57  CIN, &
58  LCL, &
59  LFC, &
60  LNB )
61  use scale_process, only: &
62  prc_mpistop
63  use scale_const, only: &
64  grav => const_grav
65  use scale_atmos_hydrometeor, only: &
66  hydrometeor_entr => atmos_hydrometeor_entr, &
67  i_qv, &
68  i_qc
69  use scale_atmos_saturation, only: &
70  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq
71  use scale_history, only: &
72  hist_in
73  implicit none
74 
75  integer, intent(in) :: kstr
76  real(RP), intent(in) :: dens(ka,ia,ja)
77  real(RP), intent(in) :: temp(ka,ia,ja)
78  real(RP), intent(in) :: pres(ka,ia,ja)
79  real(RP), intent(in) :: qtrc(ka,ia,ja,qa)
80  real(RP), intent(in) :: cz ( ka,ia,ja)
81  real(RP), intent(in) :: fz (0:ka,ia,ja)
82  real(RP), intent(out) :: cape(ia,ja)
83  real(RP), intent(out) :: cin (ia,ja)
84  real(RP), intent(out) :: lcl (ia,ja)
85  real(RP), intent(out) :: lfc (ia,ja)
86  real(RP), intent(out) :: lnb (ia,ja)
87 
88  real(RP) :: dens_p (ka,ia,ja)
89  real(RP) :: temp_p (ka,ia,ja)
90  real(RP) :: qtrc_p (ka,ia,ja,qa)
91  real(RP) :: entr_p (ka,ia,ja)
92  real(RP) :: buoy_p (ka,ia,ja)
93  real(RP) :: buoy_pf(ka,ia,ja)
94  real(RP) :: qsat_p (ka,ia,ja)
95 
96  integer :: klcl(ia,ja)
97  integer :: klfc(ia,ja)
98  integer :: klnb(ia,ja)
99 
100  integer :: k, i, j
101  !---------------------------------------------------------------------------
102 
103  if ( i_qv < 0 .OR. i_qc < 0 ) then
104  write(*,*) 'xxx Qv & Qc do not exist. CAPE calculation is invalid. STOP'
105  call prc_mpistop
106  endif
107 
108  ! entropy at start point
109  call hydrometeor_entr( entr_p(kstr,:,:), & ! [OUT]
110  temp(kstr,:,:), & ! [IN]
111  pres(kstr,:,:), & ! [IN]
112  qtrc(kstr,:,:,:), & ! [IN]
113  tracer_r(:) ) ! [IN]
114 
115  ! lift parcel
116  call atmos_adiabat_liftparcel( kstr, & ! [IN]
117  temp(:,:,:), & ! [IN]
118  pres(:,:,:), & ! [IN]
119  qtrc(:,:,:,:), & ! [IN]
120  entr_p(kstr,:,:), & ! [IN]
121  dens_p(:,:,:), & ! [OUT]
122  temp_p(:,:,:), & ! [OUT]
123  qtrc_p(:,:,:,:) ) ! [OUT]
124 
125  ! entropy profile (lifted parcel)
126  call hydrometeor_entr( entr_p(:,:,:), & ! [OUT]
127  temp_p(:,:,:), & ! [IN]
128  pres(:,:,:), & ! [IN]
129  qtrc_p(:,:,:,:), & ! [IN]
130  tracer_r(:) ) ! [IN]
131 
132  ! parcel buoyancy
133  do j = js, je
134  do i = is, ie
135  do k = ks, ke
136  buoy_p(k,i,j) = grav * ( dens(k,i,j) - dens_p(k,i,j) ) / dens_p(k,i,j)
137  enddo
138  buoy_p( 1:ks-1,i,j) = 0.0_rp
139  buoy_p(ke+1:ka ,i,j) = 0.0_rp
140 
141  do k = ks, ke-1
142  buoy_pf(k,i,j) = 0.5_rp * ( buoy_p(k+1,i,j) + buoy_p(k,i,j) )
143  enddo
144  buoy_p(ke,i,j) = buoy_p(ke-1,i,j)
145 
146  buoy_p( 1:ks-1,i,j) = 0.0_rp
147  buoy_p(ke+1:ka ,i,j) = 0.0_rp
148  enddo
149  enddo
150 
151  ! saturation point profile (lifted parcel)
152  call saturation_dens2qsat_liq( qsat_p(:,:,:), & ! [OUT]
153  temp_p(:,:,:), & ! [IN]
154  dens_p(:,:,:) ) ! [IN]
155 
156  call hist_in( dens_p(:,:,:), 'DENS_parcel', 'density profile in lifting parcel', 'kg/m3' )
157  call hist_in( temp_p(:,:,:), 'TEMP_parcel', 'temperature profile in lifting parcel', 'K' )
158  call hist_in( entr_p(:,:,:), 'ENTR_parcel', 'entropy profile in lifting parcel', 'J/K' )
159  call hist_in( buoy_p(:,:,:), 'BUOY_parcel', 'buoyancy profile in lifting parcel', 'm/s2' )
160  call hist_in( qsat_p(:,:,:), 'QSAT_parcel', 'saturation profile in lifting parcel', 'kg/kg' )
161 
162  ! detect layer number of LNB, LFC, LCL
163  klcl(:,:) = -1
164  do j = js, je
165  do i = is, ie
166  do k = kstr, ke
167  if ( qtrc_p(kstr,i,j,i_qv) >= qsat_p(k,i,j) ) then
168  klcl(i,j) = k
169  exit
170  endif
171  enddo
172  enddo
173  enddo
174 
175  klfc(:,:) = -1
176  do j = js, je
177  do i = is, ie
178  do k = klcl(i,j), ke
179  if ( buoy_p(k,i,j) >= 0.0_rp ) then
180  klfc(i,j) = k
181  exit
182  endif
183  enddo
184  enddo
185  enddo
186 
187  klnb(:,:) = -1
188  do j = js, je
189  do i = is, ie
190  do k = ke, kstr, -1
191  if ( buoy_p(k,i,j) >= 0.0_rp ) then
192  klnb(i,j) = k
193  exit
194  endif
195  enddo
196  enddo
197  enddo
198 
199  lcl(:,:) = 0.0_rp
200  do j = js, je
201  do i = is, ie
202  if ( klcl(i,j) >= kstr ) then
203  lcl(i,j) = cz(klcl(i,j),i,j)
204  endif
205  enddo
206  enddo
207 
208  lfc(:,:) = 0.0_rp
209  do j = js, je
210  do i = is, ie
211  if ( klfc(i,j) >= kstr ) then
212  lfc(i,j) = cz(klfc(i,j),i,j)
213  endif
214  enddo
215  enddo
216 
217  lnb(:,:) = 0.0_rp
218  do j = js, je
219  do i = is, ie
220  if ( klnb(i,j) >= kstr ) then
221  lnb(i,j) = cz(klnb(i,j),i,j)
222  endif
223  enddo
224  enddo
225 
226  cape(:,:) = 0.0_rp
227  do j = js, je
228  do i = is, ie
229  if ( klfc(i,j) >= kstr .AND. klnb(i,j) > kstr ) then
230  do k = klfc(i,j), klnb(i,j)
231  cape(i,j) = cape(i,j) + buoy_pf(k-1,i,j) * ( fz(k,i,j)-fz(k-1,i,j) )
232  enddo
233  endif
234  enddo
235  enddo
236 
237  cin(:,:) = 0.0_rp
238  do j = js, je
239  do i = is, ie
240  if ( klfc(i,j) >= kstr ) then
241  do k = kstr+1, klfc(i,j)
242  cin(i,j) = cin(i,j) + buoy_pf(k-1,i,j) * ( fz(k,i,j)-fz(k-1,i,j) )
243  enddo
244  endif
245  enddo
246  enddo
247 
248  return
249  end subroutine atmos_adiabat_cape
250 
251  !-----------------------------------------------------------------------------
255  subroutine atmos_adiabat_liftparcel( &
256  Kstr, &
257  TEMP, &
258  PRES, &
259  QTRC, &
260  ENTR_p, &
261  DENS_p, &
262  TEMP_p, &
263  QTRC_p )
264  use scale_const, only: &
265  eps => const_eps, &
266  rdry => const_rdry, &
267  cpdry => const_cpdry, &
268  rvap => const_rvap, &
269  cpvap => const_cpvap, &
270  lhv0 => const_lhv0, &
271  psat0 => const_psat0, &
272  pre00 => const_pre00, &
273  tem00 => const_tem00
274  use scale_atmos_hydrometeor, only: &
275  hydrometeor_entr => atmos_hydrometeor_entr, &
276  i_qv, &
277  i_qc
278  use scale_atmos_saturation, only: &
279  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq
280  implicit none
281 
282  integer, intent(in) :: kstr
283  real(RP), intent(in) :: temp (ka,ia,ja)
284  real(RP), intent(in) :: pres (ka,ia,ja)
285  real(RP), intent(in) :: qtrc (ka,ia,ja,qa)
286  real(RP), intent(in) :: entr_p(ia,ja)
287  real(RP), intent(out) :: dens_p(ka,ia,ja)
288  real(RP), intent(out) :: temp_p(ka,ia,ja)
289  real(RP), intent(out) :: qtrc_p(ka,ia,ja,qa)
290 
291  real(RP) :: qsat_p(ka,ia,ja)
292  real(RP) :: qtot_p(ia,ja)
293 
294  real(RP) :: qdry_p, rtot, cptot
295  real(RP) :: pres_dry, pres_vap
296  real(RP) :: temp_unsat
297  real(RP) :: qsat_unsat
298 
299  real(RP) :: temp_ite
300  real(RP) :: qsat_ite
301  real(RP) :: qtrc_ite(qa)
302  real(RP) :: entr_ite
303  real(RP) :: temp_prev
304  real(RP) :: dentr_dt
305 
306  real(RP), parameter :: criteria = 1.e-8_rp
307  integer, parameter :: itelim = 100
308  integer :: ite
309 
310  real(RP), parameter :: temmin = 0.1_rp
311 
312  integer :: k, i, j, iqw
313  !---------------------------------------------------------------------------
314 
315  do j = js, je
316  do i = is, ie
317  do k = 1, kstr
318  temp_p(k,i,j) = temp(k,i,j)
319  do iqw = 1, qa
320  qtrc_p(k,i,j,iqw) = qtrc(k,i,j,iqw)
321  enddo
322  enddo
323  enddo
324  enddo
325 
326  ! vapor + cloud water at the start point
327  do j = js, je
328  do i = is, ie
329  qtot_p(i,j) = qtrc(kstr,i,j,i_qv) + qtrc(kstr,i,j,i_qc)
330  enddo
331  enddo
332 
333  do j = js, je
334  do i = is, ie
335  do k = kstr, ke
336  temp_p(k,i,j) = temp(k,i,j) ! first guess
337 
338  ! T1: unsaturated temperature, S = U1(PRES, TEMP_unsat, qtot_p)
339  qdry_p = 1.0_rp - qtot_p(i,j)
340  rtot = rdry * qdry_p + rvap * qtot_p(i,j)
341  cptot = cpdry * qdry_p + cpvap * qtot_p(i,j)
342 
343  ! dry air + vapor
344  pres_dry = max( pres(k,i,j) * qdry_p * rdry / rtot, eps )
345  pres_vap = max( pres(k,i,j) * qtot_p(i,j) * rvap / rtot, eps )
346 
347  temp_unsat = tem00 * exp( ( entr_p(i,j) + qdry_p * rdry * log( pres_dry / pre00 ) &
348  + qtot_p(i,j) * rvap * log( pres_vap / psat0 ) &
349  - qtot_p(i,j) * lhv0 / tem00 ) / cptot )
350 
351  call saturation_pres2qsat_liq( qsat_unsat, & ! [OUT]
352  temp_unsat, & ! [IN]
353  pres(k,i,j) ) ! [IN]
354 
355  ! T2: saturated temperature, S = U2(PRES, TEMP_ite, QTRC_ite)
356  if ( qtot_p(i,j) > qsat_unsat ) then
357 
358  temp_ite = tem00 * exp( ( entr_p(i,j) + rdry * log( pres(k,i,j) / pre00 ) ) / cpdry )
359 
360  do ite = 1, itelim
361 
362  call saturation_pres2qsat_liq( qsat_ite, & ! [OUT]
363  temp_ite, & ! [IN]
364  pres(k,i,j) ) ! [IN]
365 
366  qtrc_ite(:) = 0.0_rp ! Pseudo-adiabatic: no cloud water
367  qtrc_ite(i_qv) = min( qtot_p(i,j), qsat_ite )
368 
369  call hydrometeor_entr( entr_ite, & ! [OUT]
370  temp_ite, & ! [IN]
371  pres(k,i,j), & ! [IN]
372  qtrc_ite(:), & ! [IN]
373  tracer_r(:) ) ! [IN]
374 
375  qdry_p = 1.0_rp - qtrc_ite(i_qv)
376  cptot = cpdry * qdry_p + cpvap * qtrc_ite(i_qv)
377 
378  dentr_dt = cptot / temp_ite &
379  - qtrc_ite(i_qv) * lhv0 / temp_ite**2 &
380  + qtrc_ite(i_qv) * lhv0**2 / temp_ite**3 / rvap
381  dentr_dt = max( dentr_dt, eps )
382 
383  temp_prev = temp_ite
384  temp_ite = temp_ite - ( entr_ite - entr_p(i,j) ) / dentr_dt
385  temp_ite = max( temp_ite, temmin )
386 
387  if( abs(temp_ite-temp_prev) < criteria ) exit
388 
389  enddo
390 
391  temp_p(k,i,j) = temp_ite
392 
393  endif
394 
395  ! parcel satulation point
396  call saturation_pres2qsat_liq( qsat_p(k,i,j), & ! [OUT]
397  temp_p(k,i,j), & ! [IN]
398  pres(k,i,j) ) ! [IN]
399 
400  ! update parcel vapor : remove condensed water
401  qtot_p(i,j) = min( qtot_p(i,j), qsat_p(k,i,j) )
402 
403  qtrc_p(k,i,j,i_qv) = qtot_p(i,j)
404  qtrc_p(k,i,j,i_qc) = 0.0_rp
405 
406  qdry_p = 1.0_rp - qtot_p(i,j)
407  rtot = rdry * qdry_p + rvap * qtot_p(i,j)
408 
409  dens_p(k,i,j) = pres(k,i,j) / ( rtot * temp_p(k,i,j) )
410  enddo
411  enddo
412  enddo
413 
414  return
415  end subroutine atmos_adiabat_liftparcel
416 
417 end module scale_atmos_adiabat
scale_atmos_adiabat
module ATMOSPHERE / Adiabatic process
Definition: scale_atmos_sub_adiabat.F90:10
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_const::const_psat0
real(rp), parameter, public const_psat0
saturate pressure of water vapor at 0C [Pa]
Definition: scale_const.F90:83
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_atmos_hydrometeor::i_qc
integer, public i_qc
Definition: scale_atmos_sub_hydrometeor.F90:79
scale_tracer::tracer_r
real(rp), dimension(qa_max), public tracer_r
Definition: scale_tracer.F90:41
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:92
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:35
scale_atmos_adiabat::atmos_adiabat_liftparcel
subroutine, public atmos_adiabat_liftparcel(Kstr, TEMP, PRES, QTRC, ENTR_p, DENS_p, TEMP_p, QTRC_p)
Calc temperature profile with lifting parcel Method: Pseudo-adiabatic ascend from lowermost layer of ...
Definition: scale_atmos_sub_adiabat.F90:264
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_atmos_adiabat::atmos_adiabat_cape
subroutine, public atmos_adiabat_cape(Kstr, DENS, TEMP, PRES, QTRC, CZ, FZ, CAPE, CIN, LCL, LFC, LNB)
Calc CAPE and CIN Type of parcel method: Pseudo-adiabatic ascend from lowermost layer of the model Re...
Definition: scale_atmos_sub_adiabat.F90:61
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_grid_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_grid_index.F90:54
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:77
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:90
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_sub_hydrometeor.F90:68
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:66
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14