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_const, only: &
62  grav => const_grav
63  use scale_atmos_thermodyn, only: &
64  thermodyn_entr => atmos_thermodyn_entr
65  use scale_atmos_saturation, only: &
66  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq
67  use scale_history, only: &
68  hist_in
69  implicit none
70 
71  integer, intent(in) :: Kstr
72  real(RP), intent(in) :: DENS(ka,ia,ja)
73  real(RP), intent(in) :: TEMP(ka,ia,ja)
74  real(RP), intent(in) :: PRES(ka,ia,ja)
75  real(RP), intent(in) :: QTRC(ka,ia,ja,qa)
76  real(RP), intent(in) :: CZ ( ka,ia,ja)
77  real(RP), intent(in) :: FZ (0:ka,ia,ja)
78  real(RP), intent(out) :: CAPE(ia,ja)
79  real(RP), intent(out) :: CIN (ia,ja)
80  real(RP), intent(out) :: LCL (ia,ja)
81  real(RP), intent(out) :: LFC (ia,ja)
82  real(RP), intent(out) :: LNB (ia,ja)
83 
84  real(RP) :: DENS_p (ka,ia,ja)
85  real(RP) :: TEMP_p (ka,ia,ja)
86  real(RP) :: QTRC_p (ka,ia,ja,qa)
87  real(RP) :: ENTR_p (ka,ia,ja)
88  real(RP) :: BUOY_p (ka,ia,ja)
89  real(RP) :: BUOY_pf(ka,ia,ja)
90  real(RP) :: QSAT_p (ka,ia,ja)
91 
92  integer :: kLCL(ia,ja)
93  integer :: kLFC(ia,ja)
94  integer :: kLNB(ia,ja)
95 
96  integer :: k, i, j
97  !---------------------------------------------------------------------------
98 
99  ! entropy at start point
100  call thermodyn_entr( entr_p(kstr,:,:), & ! [OUT]
101  temp(kstr,:,:), & ! [IN]
102  pres(kstr,:,:), & ! [IN]
103  qtrc(kstr,:,:,:) ) ! [IN]
104 
105  ! lift parcel
106  call atmos_adiabat_liftparcel( kstr, & ! [IN]
107  temp(:,:,:), & ! [IN]
108  pres(:,:,:), & ! [IN]
109  qtrc(:,:,:,:), & ! [IN]
110  entr_p(kstr,:,:), & ! [IN]
111  dens_p(:,:,:), & ! [OUT]
112  temp_p(:,:,:), & ! [OUT]
113  qtrc_p(:,:,:,:) ) ! [OUT]
114 
115  ! entropy profile (lifted parcel)
116  call thermodyn_entr( entr_p(:,:,:), & ! [OUT]
117  temp_p(:,:,:), & ! [IN]
118  pres(:,:,:), & ! [IN]
119  qtrc_p(:,:,:,:) ) ! [IN]
120 
121  ! parcel buoyancy
122  do j = js, je
123  do i = is, ie
124  do k = ks, ke
125  buoy_p(k,i,j) = grav * ( dens(k,i,j) - dens_p(k,i,j) ) / dens_p(k,i,j)
126  enddo
127  buoy_p( 1:ks-1,i,j) = 0.0_rp
128  buoy_p(ke+1:ka ,i,j) = 0.0_rp
129 
130  do k = ks, ke-1
131  buoy_pf(k,i,j) = 0.5_rp * ( buoy_p(k+1,i,j) + buoy_p(k,i,j) )
132  enddo
133  buoy_p(ke,i,j) = buoy_p(ke-1,i,j)
134 
135  buoy_p( 1:ks-1,i,j) = 0.0_rp
136  buoy_p(ke+1:ka ,i,j) = 0.0_rp
137  enddo
138  enddo
139 
140  ! saturation point profile (lifted parcel)
141  call saturation_dens2qsat_liq( qsat_p(:,:,:), & ! [OUT]
142  temp_p(:,:,:), & ! [IN]
143  dens_p(:,:,:) ) ! [IN]
144 
145  call hist_in( dens_p(:,:,:), 'DENS_parcel', 'density profile in lifting parcel', 'kg/m3' )
146  call hist_in( temp_p(:,:,:), 'TEMP_parcel', 'temperature profile in lifting parcel', 'K' )
147  call hist_in( entr_p(:,:,:), 'ENTR_parcel', 'entropy profile in lifting parcel', 'J/K' )
148  call hist_in( buoy_p(:,:,:), 'BUOY_parcel', 'buoyancy profile in lifting parcel', 'm/s2' )
149  call hist_in( qsat_p(:,:,:), 'QSAT_parcel', 'saturation profile in lifting parcel', 'kg/kg' )
150 
151  ! detect layer number of LNB, LFC, LCL
152  klcl(:,:) = -1
153  do j = js, je
154  do i = is, ie
155  do k = kstr, ke
156  if ( qtrc_p(kstr,i,j,i_qv) >= qsat_p(k,i,j) ) then
157  klcl(i,j) = k
158  exit
159  endif
160  enddo
161  enddo
162  enddo
163 
164  klfc(:,:) = -1
165  do j = js, je
166  do i = is, ie
167  do k = klcl(i,j), ke
168  if ( buoy_p(k,i,j) >= 0.0_rp ) then
169  klfc(i,j) = k
170  exit
171  endif
172  enddo
173  enddo
174  enddo
175 
176  klnb(:,:) = -1
177  do j = js, je
178  do i = is, ie
179  do k = ke, kstr, -1
180  if ( buoy_p(k,i,j) >= 0.0_rp ) then
181  klnb(i,j) = k
182  exit
183  endif
184  enddo
185  enddo
186  enddo
187 
188  lcl(:,:) = 0.0_rp
189  do j = js, je
190  do i = is, ie
191  if ( klcl(i,j) >= kstr ) then
192  lcl(i,j) = cz(klcl(i,j),i,j)
193  endif
194  enddo
195  enddo
196 
197  lfc(:,:) = 0.0_rp
198  do j = js, je
199  do i = is, ie
200  if ( klfc(i,j) >= kstr ) then
201  lfc(i,j) = cz(klfc(i,j),i,j)
202  endif
203  enddo
204  enddo
205 
206  lnb(:,:) = 0.0_rp
207  do j = js, je
208  do i = is, ie
209  if ( klnb(i,j) >= kstr ) then
210  lnb(i,j) = cz(klnb(i,j),i,j)
211  endif
212  enddo
213  enddo
214 
215  cape(:,:) = 0.0_rp
216  do j = js, je
217  do i = is, ie
218  if ( klfc(i,j) >= kstr .AND. klnb(i,j) > kstr ) then
219  do k = klfc(i,j), klnb(i,j)
220  cape(i,j) = cape(i,j) + buoy_pf(k-1,i,j) * ( fz(k,i,j)-fz(k-1,i,j) )
221  enddo
222  endif
223  enddo
224  enddo
225 
226  cin(:,:) = 0.0_rp
227  do j = js, je
228  do i = is, ie
229  if ( klfc(i,j) >= kstr ) then
230  do k = kstr+1, klfc(i,j)
231  cin(i,j) = cin(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  return
238  end subroutine atmos_adiabat_cape
239 
240  !-----------------------------------------------------------------------------
244  subroutine atmos_adiabat_liftparcel( &
245  Kstr, &
246  TEMP, &
247  PRES, &
248  QTRC, &
249  ENTR_p, &
250  DENS_p, &
251  TEMP_p, &
252  QTRC_p )
253  use scale_const, only: &
254  eps => const_eps, &
255  rdry => const_rdry, &
256  cpdry => const_cpdry, &
257  rvap => const_rvap, &
258  cpvap => const_cpvap, &
259  lhv0 => const_lhv0, &
260  psat0 => const_psat0, &
261  pre00 => const_pre00, &
262  tem00 => const_tem00
263  use scale_atmos_thermodyn, only: &
264  thermodyn_entr => atmos_thermodyn_entr
265  use scale_atmos_saturation, only: &
266  saturation_pres2qsat_liq => atmos_saturation_pres2qsat_liq
267  implicit none
268 
269  integer, intent(in) :: Kstr
270  real(RP), intent(in) :: TEMP (ka,ia,ja)
271  real(RP), intent(in) :: PRES (ka,ia,ja)
272  real(RP), intent(in) :: QTRC (ka,ia,ja,qa)
273  real(RP), intent(in) :: ENTR_p(ia,ja)
274  real(RP), intent(out) :: DENS_p(ka,ia,ja)
275  real(RP), intent(out) :: TEMP_p(ka,ia,ja)
276  real(RP), intent(out) :: QTRC_p(ka,ia,ja,qa)
277 
278  real(RP) :: qsat_p(ka,ia,ja)
279  real(RP) :: qtot_p(ia,ja)
280 
281  real(RP) :: qdry_p, Rtot, CPtot
282  real(RP) :: pres_dry, pres_vap
283  real(RP) :: TEMP_unsat
284  real(RP) :: qsat_unsat
285 
286  real(RP) :: TEMP_ite
287  real(RP) :: qsat_ite
288  real(RP) :: QTRC_ite(qa)
289  real(RP) :: ENTR_ite
290  real(RP) :: TEMP_prev
291  real(RP) :: dENTR_dT
292 
293  real(RP), parameter :: criteria = 1.e-8_rp
294  integer, parameter :: itelim = 100
295  integer :: ite
296 
297  real(RP), parameter :: TEMMIN = 0.1_rp
298 
299  integer :: k, i, j, iqw
300  !---------------------------------------------------------------------------
301 
302  do j = js, je
303  do i = is, ie
304  do k = 1, kstr
305  temp_p(k,i,j) = temp(k,i,j)
306  do iqw = qqs, qqe
307  qtrc_p(k,i,j,iqw) = qtrc(k,i,j,iqw)
308  enddo
309  enddo
310  enddo
311  enddo
312 
313  ! vapor + cloud water at the start point
314  do j = js, je
315  do i = is, ie
316  qtot_p(i,j) = qtrc(kstr,i,j,i_qv) + qtrc(kstr,i,j,i_qc)
317  enddo
318  enddo
319 
320  do j = js, je
321  do i = is, ie
322  do k = kstr, ke
323  temp_p(k,i,j) = temp(k,i,j) ! first guess
324 
325  ! T1: unsaturated temperature, S = U1(PRES, TEMP_unsat, qtot_p)
326  qdry_p = 1.0_rp - qtot_p(i,j)
327  rtot = rdry * qdry_p + rvap * qtot_p(i,j)
328  cptot = cpdry * qdry_p + cpvap * qtot_p(i,j)
329 
330  ! dry air + vapor
331  pres_dry = max( pres(k,i,j) * qdry_p * rdry / rtot, eps )
332  pres_vap = max( pres(k,i,j) * qtot_p(i,j) * rvap / rtot, eps )
333 
334  temp_unsat = tem00 * exp( ( entr_p(i,j) + qdry_p * rdry * log( pres_dry / pre00 ) &
335  + qtot_p(i,j) * rvap * log( pres_vap / psat0 ) &
336  - qtot_p(i,j) * lhv0 / tem00 ) / cptot )
337 
338  call saturation_pres2qsat_liq( qsat_unsat, & ! [OUT]
339  temp_unsat, & ! [IN]
340  pres(k,i,j) ) ! [IN]
341 
342  ! T2: saturated temperature, S = U2(PRES, TEMP_ite, QTRC_ite)
343  if ( qtot_p(i,j) > qsat_unsat ) then
344 
345  temp_ite = tem00 * exp( ( entr_p(i,j) + rdry * log( pres(k,i,j) / pre00 ) ) / cpdry )
346 
347  do ite = 1, itelim
348 
349  call saturation_pres2qsat_liq( qsat_ite, & ! [OUT]
350  temp_ite, & ! [IN]
351  pres(k,i,j) ) ! [IN]
352 
353  qtrc_ite(:) = 0.0_rp ! Pseudo-adiabatic: no cloud water
354  qtrc_ite(i_qv) = min( qtot_p(i,j), qsat_ite )
355 
356  call thermodyn_entr( entr_ite, & ! [OUT]
357  temp_ite, & ! [IN]
358  pres(k,i,j), & ! [IN]
359  qtrc_ite(:) ) ! [IN]
360 
361  qdry_p = 1.0_rp - qtrc_ite(i_qv)
362  cptot = cpdry * qdry_p + cpvap * qtrc_ite(i_qv)
363 
364  dentr_dt = cptot / temp_ite &
365  - qtrc_ite(i_qv) * lhv0 / temp_ite**2 &
366  + qtrc_ite(i_qv) * lhv0**2 / temp_ite**3 / rvap
367  dentr_dt = max( dentr_dt, eps )
368 
369  temp_prev = temp_ite
370  temp_ite = temp_ite - ( entr_ite - entr_p(i,j) ) / dentr_dt
371  temp_ite = max( temp_ite, temmin )
372 
373  if( abs(temp_ite-temp_prev) < criteria ) exit
374 
375  enddo
376 
377  temp_p(k,i,j) = temp_ite
378 
379  endif
380 
381  ! parcel satulation point
382  call saturation_pres2qsat_liq( qsat_p(k,i,j), & ! [OUT]
383  temp_p(k,i,j), & ! [IN]
384  pres(k,i,j) ) ! [IN]
385 
386  ! update parcel vapor : remove condensed water
387  qtot_p(i,j) = min( qtot_p(i,j), qsat_p(k,i,j) )
388 
389  qtrc_p(k,i,j,i_qv) = qtot_p(i,j)
390  qtrc_p(k,i,j,i_qc) = 0.0_rp
391 
392  qdry_p = 1.0_rp - qtot_p(i,j)
393  rtot = rdry * qdry_p + rvap * qtot_p(i,j)
394 
395  dens_p(k,i,j) = pres(k,i,j) / ( rtot * temp_p(k,i,j) )
396  enddo
397  enddo
398  enddo
399 
400  return
401  end subroutine atmos_adiabat_liftparcel
402 
403 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:86
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_tracer::qqe
integer, public qqe
Definition: scale_tracer.F90:61
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:95
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
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:253
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:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_tracer::i_qv
integer, public i_qv
Definition: scale_tracer.F90:47
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:80
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:93
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
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
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_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_tracer::qqs
integer, public qqs
Definition: scale_tracer.F90:60
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_tracer::i_qc
integer, public i_qc
Definition: scale_tracer.F90:48
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13