SCALE-RM
scale_atmos_phy_cp_kf.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
44 #include "scalelib.h"
45 module scale_atmos_phy_cp_kf
46  !------------------------------------------------------------------------------
47  !
48  !++ used modules
49  !
50  use scale_precision
51  use scale_io
52  use scale_prof
53  use scale_const, only: &
54  tem00 => const_tem00
55 
56  !------------------------------------------------------------------------------
57  implicit none
58  private
59 
60  !-----------------------------------------------------------------------------
61  !
62  !++ Public procedure
63  !
64  public :: atmos_phy_cp_kf_setup
65  public :: atmos_phy_cp_kf_tendency
66 
67  !-----------------------------------------------------------------------------
68  !
69  !++ Public pparameters & variabeles
70  !
71  !-----------------------------------------------------------------------------
72  !
73  !++ Private procedures
74  !
75  private :: cp_kf_param
76  private :: cp_kf_trigger
77  private :: cp_kf_updraft
78  private :: cp_kf_downdraft
79  private :: cp_kf_compensational
80  private :: cp_kf_calciexn
81  private :: cp_kf_precipitation_oc1973
82  private :: cp_kf_precipitation_kessler
83  private :: cp_kf_tpmix2
84  private :: cp_kf_dtfrznew
85  private :: cp_kf_prof5
86  private :: cp_kf_tpmix2dd
87  private :: cp_kf_envirtht
88  private :: cp_kf_lutab
89 
90  abstract interface
91  subroutine kf_precipitation( &
92  G, DZ, BOTERM, ENTERM, &
93  WTW, QLIQ, QICE, QNEWLQ, QNEWIC, QLQOUT, QICOUT )
94  use scale_precision
95  real(RP), INTENT(IN ) :: G, DZ,BOTERM,ENTERM
96  real(RP), INTENT(INOUT) :: WTW, QLIQ, QICE, QNEWLQ, QNEWIC, QLQOUT, QICOUT
97  end subroutine kf_precipitation
98  end interface
99 
100  !-----------------------------------------------------------------------------
101  !
102  !++ Private parameters & variables
103  !
104  ! KF subroutine look up tables VV
105  integer, private, parameter :: KFNT=250,kfnp=220
106  real(RP), private :: TTAB(KFNT,KFNP),QSTAB(KFNT,KFNP)
107  real(RP), private :: THE0K(KFNP)
108  real(RP), private :: ALU(200)
109  real(RP), private :: RDPR,RDTHK,PLUTOP,PBOT
110  real(RP), private :: GdCP
111  !
112 
115  real(RP), private, parameter :: ALIQ = 6.112e2_rp
116  real(RP), private, parameter :: BLIQ = 17.67_rp
117  real(RP), private, parameter :: CLIQ = bliq*tem00
118  real(RP), private, parameter :: DLIQ = 29.65_rp
119  real(RP), private, parameter :: XLV0 = 3.15e6_rp
120  real(RP), private, parameter :: XLV1 = 2370._rp
121  real(RP), private, parameter :: KF_EPS = 1.e-12_rp
124  real(RP), private :: RATE
126  integer , private :: TRIGGER_type
127  real(RP), private :: DELCAPE
128  real(RP), private :: DEEPLIFETIME
129  real(RP), private :: SHALLOWLIFETIME
130  real(RP), private :: DEPTH_USL
131  logical , private :: WARMRAIN
132  logical , private :: KF_LOG
133  real(RP), private :: kf_threshold
134  integer , private :: prec_type
135  logical, private :: DO_prec
136 
137  procedure(kf_precipitation), pointer,private :: CP_kf_precipitation => null()
138 
139  real(RP), private, allocatable :: lifetime (:,:)
140  integer , private, allocatable :: I_convflag(:,:)
141 
142  real(RP), private, allocatable :: deltaz (:,:,:)
143  real(RP), private, allocatable :: Z (:,:,:)
144  real(RP), private, allocatable :: deltax (:,:)
145 
146  ! history
147  integer, parameter :: hist_vmax = 14
148  integer, private :: hist_id(hist_vmax,0:1)
149  logical, private :: hist_flag
150  integer, parameter :: I_HIST_QV_FC = 1
151  integer, parameter :: I_HIST_QV_UE = 2
152  integer, parameter :: I_HIST_QV_DE = 3
153  integer, parameter :: I_HIST_QV_UD = 4
154  integer, parameter :: I_HIST_QV_DD = 5
155  integer, parameter :: I_HIST_QV_NF = 6
156  integer, parameter :: I_HIST_QC_FC = 7
157  integer, parameter :: I_HIST_QC_UD = 8
158  integer, parameter :: I_HIST_QI_FC = 9
159  integer, parameter :: I_HIST_QI_UD = 10
160  integer, parameter :: I_HIST_QR_FC = 11
161  integer, parameter :: I_HIST_QR_RF = 12
162  integer, parameter :: I_HIST_QS_FC = 13
163  integer, parameter :: I_HIST_QS_SF = 14
164  integer, parameter :: I_HIST_D = 0 ! deep convection
165  integer, parameter :: I_HIST_S = 1 ! shallow convection
166 
167  real(RP), private, allocatable :: hist_work(:,:,:,:,:)
168 
169  !------------------------------------------------------------------------------
170 contains
171  !------------------------------------------------------------------------------
175  subroutine atmos_phy_cp_kf_setup ( &
176  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
177  CZ, AREA, &
178  WARMRAIN_in )
179  use scale_prc, only: &
180  prc_abort
181  use scale_file_history, only: &
182  file_history_reg
183  implicit none
184  integer, intent(in) :: ka, ks, ke
185  integer, intent(in) :: ia, is, ie
186  integer, intent(in) :: ja, js, je
187 
188  real(rp), intent(in) :: cz(ka,ia,ja)
189  real(rp), intent(in) :: area(ia,ja)
190  logical, intent(in) :: warmrain_in
191 
192 
193  integer :: atmos_phy_cp_kf_trigger_type = 1
194  integer :: atmos_phy_cp_kf_prec_type = 1
195  real(rp) :: atmos_phy_cp_kf_rate = 0.03_rp
196  logical :: atmos_phy_cp_kf_do_prec = .true.
197  real(rp) :: atmos_phy_cp_kf_dlcape = 0.1_rp
198  real(rp) :: atmos_phy_cp_kf_dlifetime = 1800.0_rp
199  real(rp) :: atmos_phy_cp_kf_slifetime = 2400.0_rp
200  real(rp) :: atmos_phy_cp_kf_depth_usl = 300.0_rp
201  real(rp) :: atmos_phy_cp_kf_thres = 1.e-3_rp
202  logical :: atmos_phy_cp_kf_log = .false.
203 
204  namelist / param_atmos_phy_cp_kf / &
205  atmos_phy_cp_kf_rate, &
206  atmos_phy_cp_kf_trigger_type, &
207  atmos_phy_cp_kf_dlcape, &
208  atmos_phy_cp_kf_dlifetime, &
209  atmos_phy_cp_kf_slifetime, &
210  atmos_phy_cp_kf_depth_usl, &
211  atmos_phy_cp_kf_prec_type, &
212  atmos_phy_cp_kf_do_prec, &
213  atmos_phy_cp_kf_thres, &
214  atmos_phy_cp_kf_log
215 
216  integer :: k, i, j
217  integer :: n
218  integer :: ierr
219  !---------------------------------------------------------------------------
220 
221  log_newline
222  log_info("ATMOS_PHY_CP_kf_setup",*) 'Setup'
223  log_info("ATMOS_PHY_CP_kf_setup",*) 'Kain-Fritsch scheme'
224 
225  warmrain = warmrain_in
226 
227  !--- read namelist
228  rewind(io_fid_conf)
229  read(io_fid_conf,nml=param_atmos_phy_cp_kf,iostat=ierr)
230  if( ierr < 0 ) then !--- missing
231  log_info("ATMOS_PHY_CP_kf_setup",*) 'Not found namelist. Default used.'
232  elseif( ierr > 0 ) then !--- fatal error
233  log_error("ATMOS_PHY_CP_kf_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_CP_KF. Check!'
234  call prc_abort
235  endif
236  log_nml(param_atmos_phy_cp_kf)
237 
238  call cp_kf_lutab ! set up of KF look-up table
239 
240  call cp_kf_param( atmos_phy_cp_kf_prec_type, & ! [IN]
241  atmos_phy_cp_kf_rate, & ! [IN]
242  atmos_phy_cp_kf_do_prec, & ! [IN]
243  atmos_phy_cp_kf_dlcape, & ! [IN]
244  atmos_phy_cp_kf_dlifetime, & ! [IN]
245  atmos_phy_cp_kf_slifetime, & ! [IN]
246  atmos_phy_cp_kf_depth_usl, & ! [IN]
247  atmos_phy_cp_kf_thres, & ! [IN]
248  atmos_phy_cp_kf_log, & ! [IN]
249  atmos_phy_cp_kf_trigger_type ) ! [IN]
250 
251  ! output parameter lists
252  log_newline
253  log_info("ATMOS_PHY_CP_kf_setup",*) "Ogura-Cho condense material convert rate : ", atmos_phy_cp_kf_rate
254  log_info("ATMOS_PHY_CP_kf_setup",*) "Trigger function type, 1:KF2004 3:NO2007 : ", atmos_phy_cp_kf_trigger_type
255  log_info("ATMOS_PHY_CP_kf_setup",*) "CAPE decrease rate : ", atmos_phy_cp_kf_dlcape
256  log_info("ATMOS_PHY_CP_kf_setup",*) "Minimum lifetime scale of deep convection [sec] : ", atmos_phy_cp_kf_dlifetime
257  log_info("ATMOS_PHY_CP_kf_setup",*) "Lifetime of shallow convection [sec] : ", atmos_phy_cp_kf_slifetime
258  log_info("ATMOS_PHY_CP_kf_setup",*) "Updraft souce layer depth [hPa] : ", atmos_phy_cp_kf_depth_usl
259  log_info("ATMOS_PHY_CP_kf_setup",*) "Precipitation type 1:Ogura-Cho(1973) 2:Kessler : ", atmos_phy_cp_kf_prec_type
260  log_info("ATMOS_PHY_CP_kf_setup",*) "Kessler type precipitation's threshold : ", atmos_phy_cp_kf_thres
261  log_info("ATMOS_PHY_CP_kf_setup",*) "Enable sedimentation (precipitation) ? : ", atmos_phy_cp_kf_do_prec
262  log_info("ATMOS_PHY_CP_kf_setup",*) "Warm rain? : ", warmrain
263  log_info("ATMOS_PHY_CP_kf_setup",*) "Output log? : ", atmos_phy_cp_kf_log
264 
265  ! output variables
266  allocate( lifetime(ia,ja) )
267  allocate( i_convflag(ia,ja) )
268  lifetime(:,:) = 0.0_rp
269  i_convflag(:,:) = 2
270 
271  allocate( z(ka,ia,ja) )
272  z(:,:,:) = cz(:,:,:) ! because scale_atmos_phy_cp interface ,not use scale_grid
273 
274  allocate( deltaz(ka,ia,ja) )
275  ! deltaz is the interval of between model full levels(scalar point )
276  deltaz(:,:,:) = 0.0_rp
277  do j = js, je
278  do i = is, ie
279  do k = ks, ke
280  deltaz(k,i,j) = cz(k+1,i,j) - cz(k,i,j)
281  enddo
282  enddo
283  enddo
284  deltaz(ke,:,:) = 0.0_rp
285 
286  allocate( deltax(ia,ja) )
287  do j = js, je
288  do i = is, ie
289  deltax(i,j) = sqrt( area(i,j) )
290  end do
291  end do
292 
293 
294  ! history
295  hist_id(:,:) = -1
296  call file_history_reg( 'QV_t_KF_conv_DP', &
297  'QV tendency due to flux convergence (deep convection) in KF', &
298  'kg/m3/s', &
299  hist_id(i_hist_qv_fc,i_hist_d) )
300  call file_history_reg( 'QV_t_KF_conv_SH', &
301  'QV tendency due to flux convergence (shallow convection) in KF', &
302  'kg/m3/s', &
303  hist_id(i_hist_qv_fc,i_hist_s) )
304  call file_history_reg( 'QV_t_KF_upent_DP', &
305  'QV tendency due to entrainment of updraft (deep convection) in KF', &
306  'kg/m3/s', &
307  hist_id(i_hist_qv_ue,i_hist_d) )
308  call file_history_reg( 'QV_t_KF_upent_SH', &
309  'QV tendency due to entrainment of updraft (shallow convection) in KF', &
310  'kg/m3/s', &
311  hist_id(i_hist_qv_ue,i_hist_s) )
312  call file_history_reg( 'QV_t_KF_downent_DP', &
313  'QV tendency due to entrainment of downdraft (deep convection) in KF', &
314  'kg/m3/s', &
315  hist_id(i_hist_qv_de,i_hist_d) )
316  call file_history_reg( 'QV_t_KF_updet_DP', &
317  'QV tendency due to detrainment of updraft (deep convection) in KF', &
318  'kg/m3/s', &
319  hist_id(i_hist_qv_ud,i_hist_d) )
320  call file_history_reg( 'QV_t_KF_updet_SH', &
321  'QV tendency due to detrainment of updraft (shallow convection) in KF', &
322  'kg/m3/s', &
323  hist_id(i_hist_qv_ud,i_hist_s) )
324  call file_history_reg( 'QV_t_KF_downdet_DP', &
325  'QV tendency due to detrainment of downdraft (deep convection) in KF', &
326  'kg/m3/s', &
327  hist_id(i_hist_qv_dd,i_hist_d) )
328  call file_history_reg( 'QV_t_KF_negfix_DP', &
329  'QV tendency due to negative fixer (deep convection) in KF', &
330  'kg/m3/s', &
331  hist_id(i_hist_qv_nf,i_hist_d) )
332  call file_history_reg( 'QV_t_KF_negfix_SH', &
333  'QV tendency due to negative fixer (shallow convection) in KF', &
334  'kg/m3/s', &
335  hist_id(i_hist_qv_nf,i_hist_s) )
336 
337  call file_history_reg( 'QC_t_KF_conv_DP', &
338  'QC tendency due to flux convergence (deep convection) in KF', &
339  'kg/m3/s', &
340  hist_id(i_hist_qc_fc,i_hist_d) )
341  call file_history_reg( 'QC_t_KF_conv_SH', &
342  'QC tendency due to flux convergence (shallow convection) in KF', &
343  'kg/m3/s', &
344  hist_id(i_hist_qc_fc,i_hist_s) )
345  call file_history_reg( 'QC_t_KF_updet_DP', &
346  'QC tendency due to detrainment of updraft (deep convection) in KF', &
347  'kg/m3/s', &
348  hist_id(i_hist_qc_ud,i_hist_d) )
349  call file_history_reg( 'QC_t_KF_updet_SH', &
350  'QC tendency due to detrainment of updraft (shallow convection) in KF', &
351  'kg/m3/s', &
352  hist_id(i_hist_qc_ud,i_hist_s) )
353 
354  call file_history_reg( 'QI_t_KF_conv_DP', &
355  'QI tendency due to flux convergence (deep convection) in KF', &
356  'kg/m3/s', &
357  hist_id(i_hist_qi_fc,i_hist_d) )
358  call file_history_reg( 'QI_t_KF_conv_SH', &
359  'QI tendency due to flux convergence (shallow convection) in KF', &
360  'kg/m3/s', &
361  hist_id(i_hist_qi_fc,i_hist_s) )
362  call file_history_reg( 'QI_t_KF_updet_DP', &
363  'QI tendency due to detrainment of updraft (deep convection) in KF', &
364  'kg/m3/s', &
365  hist_id(i_hist_qi_ud,i_hist_d) )
366  call file_history_reg( 'QI_t_KF_updet_SH', &
367  'QI tendency due to detrainment of updraft (shallow convection) in KF', &
368  'kg/m3/s', &
369  hist_id(i_hist_qi_ud,i_hist_s) )
370 
371  call file_history_reg( 'QR_t_KF_conv_DP', &
372  'QR tendency due to flux convergence (deep convection) in KF', &
373  'kg/m3/s', &
374  hist_id(i_hist_qr_fc,i_hist_d) )
375  call file_history_reg( 'QR_t_KF_conv_SH', &
376  'QR tendency due to flux convergence (shallow convection) in KF', &
377  'kg/m3/s', &
378  hist_id(i_hist_qr_fc,i_hist_s) )
379  call file_history_reg( 'QR_t_KF_rainfall_DP', &
380  'QR tendency due to rain fall (deep convection) in KF', &
381  'kg/m3/s', &
382  hist_id(i_hist_qr_rf,i_hist_d) )
383  call file_history_reg( 'QR_t_KF_rainfall_SH', &
384  'QR tendency due to rain fall (shallow convection) in KF', &
385  'kg/m3/s', &
386  hist_id(i_hist_qr_rf,i_hist_s) )
387 
388  call file_history_reg( 'QS_t_KF_conv_DP', &
389  'QS tendency due to flux convergence (deep convection) in KF', &
390  'kg/m3/s', &
391  hist_id(i_hist_qs_fc,i_hist_d) )
392  call file_history_reg( 'QS_t_KF_conv_SH', &
393  'QS tendency due to flux convergence (shallow convection) in KF', &
394  'kg/m3/s', &
395  hist_id(i_hist_qs_fc,i_hist_s) )
396  call file_history_reg( 'QS_t_KF_snowfall_DP', &
397  'QS tendency due to snow fall (deep convection) in KF', &
398  'kg/m3/s', &
399  hist_id(i_hist_qs_sf,i_hist_d) )
400  call file_history_reg( 'QS_t_KF_snowfall_SH', &
401  'QS tendency due to snow fall (shallow convection) in KF', &
402  'kg/m3/s', &
403  hist_id(i_hist_qs_sf,i_hist_s) )
404 
405  do n = 1, hist_vmax
406  if ( hist_id(n,0) > 0 .or. hist_id(n,1) > 0 ) then
407  allocate( hist_work(ka,ia,ja,hist_vmax,0:1) )
408  hist_work(:,:,:,:,:) = 0.0_rp
409  exit
410  end if
411  end do
412 
413  return
414  end subroutine atmos_phy_cp_kf_setup
415 
416  !------------------------------------------------------------------------------
420  subroutine cp_kf_param( & ! set kf tuning parametres
421  prec_type_in, &
422  RATE_in, &
423  DO_PREC_in, &
424  DELCAPE_in , &
425  DEEPLIFETIME_in, &
426  SHALLOWLIFETIME_in, &
427  DEPTH_USL_in, &
428  KF_threshold_in, &
429  KF_LOG_in, &
430  TRIGGER_type_in )
431  use scale_prc, only: &
432  prc_abort
433  implicit none
434  integer, intent(in) :: prec_type_in
435  real(rp), intent(in) :: rate_in
436  logical, intent(in) :: do_prec_in
437  real(rp), intent(in) :: delcape_in
438  real(rp), intent(in) :: deeplifetime_in
439  real(rp), intent(in) :: shallowlifetime_in
440  real(rp), intent(in) :: depth_usl_in
441  real(rp), intent(in) :: kf_threshold_in
442  logical, intent(in) :: kf_log_in
443  integer, intent(in) :: trigger_type_in
444  !
445  rate = rate_in
446  delcape = delcape_in
447  deeplifetime = deeplifetime_in
448  shallowlifetime = shallowlifetime_in
449  depth_usl = depth_usl_in
450  prec_type = prec_type_in
451  do_prec = do_prec_in
452  kf_threshold = kf_threshold_in
453  kf_log = kf_log_in
454  trigger_type = trigger_type_in
455 
456  if ( trigger_type /= 1 .and. trigger_type /=3 ) then
457  log_error("CP_kf_param",*) "TRIGGER_type must be 1 or 3 but now : ",trigger_type
458  call prc_abort
459  end if
460 
461  select case ( prec_type )
462  case (1)
463  cp_kf_precipitation => cp_kf_precipitation_oc1973 ! Ogura and Cho (1973)
464  case (2)
465  cp_kf_precipitation => cp_kf_precipitation_kessler ! Kessler type
466  case default
467  log_error("CP_kf_param",*) 'KF namelist'
468  log_error_cont(*) 'prec_type must be 1 or 2 : ', prec_type
469  call prc_abort
470  end select
471 
472  return
473  end subroutine cp_kf_param
474 
475  !------------------------------------------------------------------------------
479  subroutine atmos_phy_cp_kf_tendency( &
480  KA, KS, KE, &
481  IA, IS, IE, &
482  JA, JS, JE, &
483  DENS, &
484  U, V, &
485  RHOT, &
486  TEMP, PRES, &
487  QDRY, QV_in, &
488  w0avg, &
489  FZ, &
490  KF_DTSECD, &
491  DENS_t, &
492  RHOT_t, &
493  RHOQV_t, &
494  RHOQ_t, &
495  nca, &
496  SFLX_rain, &
497  SFLX_snow, &
498  SFLX_engi, &
499  cloudtop, &
500  cloudbase, &
501  cldfrac_dp, &
502  cldfrac_sh )
503  use scale_file_history, only: &
504  file_history_query, &
505  file_history_put, &
506  file_history_in
507  use scale_const, only: &
508  grav => const_grav, &
509  pre00 => const_pre00, &
510  epsvap => const_epsvap
511  use scale_atmos_hydrometeor, only: &
512  cp_vapor, &
513  cp_water, &
514  cp_ice, &
515  cv_vapor, &
516  cv_water, &
517  cv_ice, &
518  n_hyd, &
519  i_hc, &
520  i_hr, &
521  i_hi, &
522  i_hs
523  use scale_atmos_saturation ,only :&
524  saturation_psat_liq => atmos_saturation_psat_liq
525  implicit none
526  integer, intent(in) :: ka, ks, ke
527  integer, intent(in) :: ia, is, ie
528  integer, intent(in) :: ja, js, je
529 
530  real(rp), intent(in) :: dens (ka,ia,ja)
531  real(rp), intent(in) :: u (ka,ia,ja)
532  real(rp), intent(in) :: v (ka,ia,ja)
533  real(rp), intent(in) :: rhot (ka,ia,ja)
534  real(rp), intent(in) :: temp (ka,ia,ja)
535  real(rp), intent(in) :: pres (ka,ia,ja)
536  real(rp), intent(in) :: qdry (ka,ia,ja)
537  real(rp), intent(in) :: qv_in(ka,ia,ja)
538  real(rp), intent(in) :: w0avg(ka,ia,ja)
539  real(rp), intent(in) :: fz (0:ka,ia,ja)
540  real(dp), intent(in) :: kf_dtsecd
541 
542  real(rp), intent(inout) :: dens_t (ka,ia,ja)
543  real(rp), intent(inout) :: rhot_t (ka,ia,ja)
544  real(rp), intent(inout) :: rhoqv_t(ka,ia,ja)
545  real(rp), intent(inout) :: rhoq_t (ka,ia,ja,n_hyd)
546  real(rp), intent(inout) :: nca (ia,ja)
547 
548  real(rp), intent(inout) :: sflx_rain(ia,ja)
549  real(rp), intent(inout) :: sflx_snow(ia,ja)
550  real(rp), intent(inout) :: sflx_engi(ia,ja)
551  real(rp), intent(inout) :: cloudtop (ia,ja)
552  real(rp), intent(inout) :: cloudbase (ia,ja)
553  real(rp), intent(inout) :: cldfrac_dp (ka,ia,ja)
554  real(rp), intent(inout) :: cldfrac_sh (ka,ia,ja)
555 
556  integer :: k, i, j, iq
557  integer :: nic
558  integer :: k_lcl
559  integer :: k_top
560  integer :: k_ml
561  integer :: k_lc,k_let,k_pbl
562  integer :: k_lfs
563 
564  real(rp) :: qv (ka)
565  real(rp) :: psat
566  real(rp) :: qsat (ka)
567  real(rp) :: rh (ka)
568  real(rp) :: deltap(ka)
569 
570  real(rp) :: cldfrac_kf(ka,2)
571  real(rp) :: dens_nw(ka)
572  real(rp) :: time_advec
573  real(rp) :: umf(ka)
574  real(rp) :: umflcl
575  real(rp) :: upent(ka)
576  real(rp) :: updet(ka)
577  ! updraft value
578  real(rp) :: temp_u(ka)
579  real(rp) :: tempv(ka)
580  real(rp) :: qv_u(ka)
581  real(rp) :: cape
582  ! water
583  real(rp) :: qc(ka)
584  real(rp) :: qi(ka)
585  real(rp) :: qvdet(ka)
586  real(rp) :: qcdet(ka)
587  real(rp) :: qidet(ka)
588  real(rp) :: flux_qr(ka)
589  real(rp) :: flux_qs(ka)
590  ! upward theta
591  real(rp) :: theta_eu(ka)
592  real(rp) :: theta_ee(ka)
593  real(rp) :: zmix
594  real(rp) :: presmix
595  real(rp) :: umfnewdold(ka)
596  real(rp) :: dpthmx
597  real(rp) :: ems(ka)
598  real(rp) :: emsd(ka)
599  ! output from downdraft
600  real(rp) :: wspd(3)
601  real(rp) :: dmf(ka)
602  real(rp) :: downent(ka)
603  real(rp) :: downdet(ka)
604  real(rp) :: theta_d(ka)
605  real(rp) :: qv_d(ka)
606  real(rp) :: rain_flux
607  real(rp) :: snow_flux
608  real(rp) :: prec_engi
609 
610  ! update variables
611  real(rp) :: theta_nw(ka)
612  real(rp) :: qv_nw(ka)
613  real(rp) :: qc_nw(ka)
614  real(rp) :: qi_nw(ka)
615  real(rp) :: qr_nw(ka)
616  real(rp) :: qs_nw(ka)
617 
618  real(rp) :: rhod(ka)
619 
620  real(rp) :: dqv, dqc, dqi, dqr, dqs
621 
622  real(rp) :: r_convflag(ia,ja)
623 
624  real(rp) :: kf_dtsec
625 
626  logical :: flag
627  integer :: n, m
628  ! ------
629 
630  log_progress(*) 'atmosphere / physics / cumulus / KF'
631  log_info("ATMOS_PHY_CP_kf_tendency",*) 'KF Convection Check '
632 
633  kf_dtsec = kf_dtsecd ! DP to RP
634 
635  call prof_rapstart('CP_kf', 3)
636 
637  do n = 1, hist_vmax
638  call file_history_query( hist_id(n,0), hist_flag )
639  if ( hist_flag ) exit
640  call file_history_query( hist_id(n,1), hist_flag )
641  if ( hist_flag ) exit
642  end do
643 
644  !$omp parallel do default(none) &
645  !$omp private(RHOD,tempv,qv_d,qc,qi,PSAT,QSAT,QV,rh, &
646  !$omp dens_nw,theta_nw,qv_nw,qc_nw,qi_nw,qr_nw,qs_nw, &
647  !$omp flux_qr,flux_qs,theta_eu,theta_ee,theta_d,umfnewdold,qvdet,qcdet,qidet, &
648  !$omp umf,umflcl,cape,presmix,upent,updet,temp_u,qv_u, &
649  !$omp dpthmx,zmix,rain_flux,snow_flux,prec_engi,time_advec, &
650  !$omp k_lcl,k_lc,k_lfs,k_pbl,k_top,k_let,k_ml, &
651  !$omp nic,deltap,cldfrac_KF,ems,emsd,wspd,dmf,downent,downdet, &
652  !$omp dQv,dQC,dQR,dQI,dQS) &
653  !$omp shared(DENS,RHOT,U,V,QDRY,TEMP,PRES,QV_in,FZ,Z,w0avg, &
654  !$omp DENS_t,RHOT_t,RHOQV_t,RHOQ_t, &
655  !$omp GRAV,PRE00,CP_VAPOR,CP_WATER,CP_ICE,CV_VAPOR,CV_WATER,CV_ICE,EPSvap, &
656  !$omp WARMRAIN,KF_DTSEC,SHALLOWLIFETIME, &
657  !$omp KA,KS,KE,IS,IE,JS,JE, &
658  !$omp nca,cloudtop,cloudbase,SFLX_rain,SFLX_snow,SFLX_engi, &
659  !$omp lifetime,deltaz,deltax,I_convflag, &
660  !$omp cldfrac_sh,cldfrac_dp,hist_flag,hist_work)
661  do j = js, je
662  do i = is, ie
663 
664  nca(i,j) = nca(i,j) - kf_dtsec
665 
666  ! check convection
667  if ( nca(i,j) .ge. 0.5_rp * kf_dtsec ) cycle
668 
669  do k = ks, ke
670  ! preparing a NON Hydriometeor condition to fit assumption in KF scheme
671  rhod(k) = dens(k,i,j) * qdry(k,i,j)
672  enddo
673 
674 
675  do k = ks, ke
676  ! temporary: WRF TYPE equations are used to maintain consistency with kf_main
677  !call SATURATION_psat_liq( PSAT, TEMP(k,i,j) )
678  !QSAT(k) = 0.622_RP * PSAT / ( PRES(k,i,j) - ( 1.0_RP-0.622_RP ) * PSAT )
679  psat = aliq*exp((bliq*temp(k,i,j)-cliq)/(temp(k,i,j)-dliq))
680  qsat(k) = epsvap * psat / ( pres(k,i,j) - psat )
681 
682  ! calculate water vaper and relative humidity
683  qv(k) = max( kf_eps, min( qsat(k), qv_in(k,i,j) / qdry(k,i,j) ) ) ! compare QSAT and QV, guess lower limit
684  rh(k) = qv(k) / qsat(k)
685  enddo
686 
687  ! calculate delta P by hydrostatic balance
688  ! deltap is the pressure interval between half levels(face levels) @ SCALE
689  do k = ks, ke
690  deltap(k) = rhod(k) * grav * ( fz(k,i,j) - fz(k-1,i,j) ) ! rho*g*dz
691  enddo
692 
693  call cp_kf_trigger ( &
694  ka, ks, ke, & ! [IN]
695  deltaz(:,i,j), z(:,i,j), & ! [IN]
696  qv(:), qsat(:), & ! [IN]
697  pres(:,i,j), & ! [IN]
698  deltap(:), deltax(i,j), & ! [IN]
699  temp(:,i,j), & ! [IN]
700  w0avg(:,i,j), & ! [IN]
701  i_convflag(i,j), & ! [OUT]
702  cloudtop(i,j), & ! [OUT]
703  temp_u(:), tempv(:), & ! [OUT]
704  qv_u(:), qc(:), qi(:), & ! [OUT]
705  qvdet(:), qcdet(:), qidet(:), & ! [OUT]
706  flux_qr(:), flux_qs(:), & ! [OUT]
707  theta_eu(:), theta_ee(:), & ! [OUT]
708  cape, & ! [OUT]
709  umf(:), umflcl, & ! [OUT]
710  upent(:), updet(:), & ! [OUT]
711  k_lcl, k_lc, k_pbl, & ! [OUT]
712  k_top, k_let, k_ml, & ! [OUT]
713  presmix, & ! [OUT]
714  dpthmx, & ! [OUT]
715  cloudbase(i,j), zmix, & ! [OUT]
716  umfnewdold(:) ) ! [OUT]
717 
718  if (i_convflag(i,j) /= 2) then ! convection allowed I_convflag=0 or 1
719 
720  ! calc ems(box weight[kg])
721  ems(k_top+1:ke) = 0._rp
722  emsd(k_top+1:ke) = 0._rp
723  do k = ks, k_top
724  ems(k) = deltap(k) * deltax(i,j)**2 / grav
725  emsd(k) = 1._rp/ems(k)
726  umfnewdold(k) = 1._rp/umfnewdold(k)
727  end do
728 
729  call cp_kf_downdraft ( &
730  ka, ks, ke, & ! [IN]
731  i_convflag(i,j), & ! [IN]
732  k_lcl, k_ml, k_top, k_pbl, k_let, k_lc, & ! [IN]
733  z(:,i,j), cloudbase(i,j), & ! [IN]
734  u(:,i,j), v(:,i,j), rh(:), qv(:), pres(:,i,j), & ! [IN]
735  deltap(:), deltax(i,j), & ! [IN]
736  ems(:), & ! [IN]
737  theta_ee(:), & ! [IN]
738  umf(:), temp_u(:), & ! [IN]
739  flux_qr(:), flux_qs(:), tempv(:), & ! [IN]
740  wspd(:), dmf(:), downent(:), downdet(:), & ! [OUT]
741  theta_d(:), qv_d(:), & ! [OUT]
742  rain_flux, snow_flux, prec_engi, & ! [OUT]
743  k_lfs ) ! [OUT]
744 
745  call cp_kf_compensational ( &
746  ka, ks, ke, & ! [IN]
747  k_top, k_lc, k_pbl, k_ml, k_lfs, & ! [IN]
748  deltaz(:,i,j), z(:,i,j), pres(:,i,j), deltap(:), deltax(i,j), & ! [IN]
749  temp(:,i,j), qv(:), ems(:), emsd(:), & ! [IN]
750  presmix, zmix, dpthmx, & ! [IN]
751  cape, & ! [IN]
752  temp_u(:), qvdet(:), umflcl, & ! [IN]
753  qc(:), qi(:), flux_qr(:), flux_qs(:), & ! [IN]
754  umfnewdold(:), & ! [IN]
755  wspd(:), & ! [IN]
756  qv_d(:), theta_d(:), & ! [IN]
757  prec_engi, & ! [IN]
758  kf_dtsec, & ! [IN]
759  rhod(:), i, j, & ! [IN]
760  i_convflag(i,j), k_lcl, & ! [INOUT]
761  umf(:), upent(:), updet(:), & ! [INOUT]
762  qcdet(:), qidet(:), dmf(:), downent(:), downdet(:), & ! [INOUT]
763  rain_flux, snow_flux, & ! [INOUT]
764  nic, & ! [OUT]
765  theta_nw(:), & ! [OUT]
766  qv_nw(:), qc_nw(:), qi_nw(:), qr_nw(:), qs_nw(:), & ! [OUT]
767  sflx_rain(i,j), sflx_snow(i,j), sflx_engi(i,j), & ! [OUT]
768  cldfrac_kf, lifetime(i,j), time_advec ) ! [OUT]
769 
770  end if
771 
772  if (i_convflag(i,j) == 2) then ! no convection
773 
774  sflx_rain(i,j) = 0.0_rp
775  sflx_snow(i,j) = 0.0_rp
776  sflx_engi(i,j) = 0.0_rp
777  cldfrac_kf(ks:ke,:) = 0.0_rp
778  lifetime(i,j) = 0.0_rp
779  cloudtop(i,j) = 0.0_rp
780  cloudbase(i,j) = 0.0_rp
781  nca(i,j) = -100.0_rp
782 
783  dens_t(ks:ke,i,j) = 0.0_rp
784  rhot_t(ks:ke,i,j) = 0.0_rp
785  rhoqv_t(ks:ke,i,j) = 0.0_rp
786  do iq = 1, n_hyd
787  rhoq_t(ks:ke,i,j,iq) = 0.0_rp
788  end do
789 
790  if ( hist_flag ) then
791  do m = 0, 1
792  do n = 1, hist_vmax
793  hist_work(:,i,j,n,m) = 0.0_rp
794  end do
795  end do
796  end if
797 
798  else
799 
800  ! compute tendencys
801 
802  ! check:
803  ! FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
804  ! IF THE ADVECTIVE TIME PERIOD (time_advec) IS LESS THAN SPECIFIED MINIMUM
805  ! lifetime, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC.
806  !
807  if (i_convflag(i,j) == 0) then ! deep
808  if (time_advec < lifetime(i,j)) nic=nint(time_advec/kf_dtsec)
809  nca(i,j) = nic * kf_dtsec ! convection feed back act this time span
810  elseif (i_convflag(i,j) == 1) then ! shallow
811  lifetime(i,j) = max(shallowlifetime, kf_dtsec)
812  nca(i,j) = kf_dtsec ! convection feed back act this time span
813  end if
814 
815  do k = ks, k_top
816  ! vapor
817  dqv = rhod(k) * ( qv_nw(k) - qv(k) )
818  rhoqv_t(k,i,j) = dqv / lifetime(i,j)
819  ! liquid water
820  dqc = qc_nw(k) * rhod(k)
821  dqr = qr_nw(k) * rhod(k)
822  rhoq_t(k,i,j,i_hc) = dqc / lifetime(i,j)
823  rhoq_t(k,i,j,i_hr) = dqr / lifetime(i,j)
824  ! ice water
825  if ( .not. warmrain ) then
826  dqi = qi_nw(k) * rhod(k)
827  dqs = qs_nw(k) * rhod(k)
828  else
829  dqi = 0.0_rp
830  dqs = 0.0_rp
831  end if
832  rhoq_t(k,i,j,i_hi) = dqi / lifetime(i,j)
833  rhoq_t(k,i,j,i_hs) = dqs / lifetime(i,j)
834 
835  dens_t(k,i,j) = rhoqv_t(k,i,j) &
836  + rhoq_t(k,i,j,i_hc) + rhoq_t(k,i,j,i_hr) &
837  + rhoq_t(k,i,j,i_hi) + rhoq_t(k,i,j,i_hs)
838  dens_nw(k) = dens(k,i,j) + dqv + dqc + dqr + dqi + dqs
839 
840  end do
841  do iq = i_hs+1, n_hyd
842  do k = ks, k_top
843  rhoq_t(k,i,j,iq) = 0.0_rp
844  end do
845  end do
846 
847  ! internal energy
848  do k = ks, k_top
849  rhot_t(k,i,j) = ( dens_nw(k) * theta_nw(k) - rhot(k,i,j) ) / lifetime(i,j)
850  end do
851 
852  do k=k_top+1, ke
853  dens_t(k,i,j) = 0.0_rp
854  rhot_t(k,i,j) = 0.0_rp
855  rhoqv_t(k,i,j) = 0.0_rp
856  do iq = 1, n_hyd
857  rhoq_t(k,i,j,iq) = 0.0_rp
858  end do
859  end do
860 
861  ! to keep conservation
862  ! if noconvection then nca is same value before call. nca only modifyed convectioned
863  end if
864 
865  cldfrac_sh(ks:ke,i,j) = cldfrac_kf(ks:ke,1)
866  cldfrac_dp(ks:ke,i,j) = cldfrac_kf(ks:ke,2)
867  end do
868  end do
869 
870  call prof_rapend('CP_kf', 3)
871 
872  ! history
873  call file_history_in( lifetime(:,:), 'KF_LIFETIME', 'lifetime of KF scheme', 's' )
874  r_convflag(:,:) = real(i_convflag(:,:),rp)
875  call file_history_in( r_convflag(:,:), 'KF_CONVFLAG', 'CONVECTION FLAG', '' )
876  do m = 0, 1
877  do n = 1, hist_vmax
878  call file_history_query( hist_id(n,m), flag )
879  if ( flag ) then
880  call file_history_put( hist_id(n,m), hist_work(:,:,:,n,m) )
881  end if
882  end do
883  end do
884 
885  return
886  end subroutine atmos_phy_cp_kf_tendency
887 
888  !------------------------------------------------------------------------------
892 !OCL SERIAL
893  subroutine cp_kf_trigger ( &
894  KA, KS, KE, &
895  dz_kf, z_kf, &
896  qv, qes, &
897  pres, &
898  deltap, deltax, &
899  temp, &
900  w0avg, &
901  I_convflag, &
902  cloudtop, &
903  temp_u, tempv, &
904  qv_u, qc, qi, &
905  qvdet, qcdet, qidet, &
906  flux_qr, flux_qs, &
907  theta_eu, theta_ee, &
908  cape, &
909  umf, umflcl, &
910  upent, updet, &
911  k_lcl, k_lc, k_pbl, &
912  k_top, k_let, k_ml, &
913  presmix, &
914  dpthmx, &
915  zlcl, zmix, &
916  umfnewdold )
917  use scale_const,only :&
918  tem00 => const_tem00, &
919  grav => const_grav, &
920  epsvap => const_epsvap, &
921  epstvap => const_epstvap
922  use scale_prc, only: &
923  prc_abort
924  implicit none
925  integer, intent(in) :: ka, ks, ke
926  real(rp), intent(in) :: dz_kf(ka),z_kf(ka)
927  real(rp), intent(in) :: qv(ka)
928  real(rp), intent(in) :: qes(ka)
929  real(rp), intent(in) :: pres(ka)
930  real(rp), intent(in) :: deltap(ka)
931  real(rp), intent(in) :: deltax
932  real(rp), intent(in) :: temp(ka)
933  real(rp), intent(in) :: w0avg(ka)
934 
935  real(rp), intent(out) :: umf(ka)
936  real(rp), intent(out) :: umflcl
937  real(rp), intent(out) :: upent(ka)
938  real(rp), intent(out) :: updet(ka)
939  real(rp), intent(out) :: temp_u(ka)
940  real(rp), intent(out) :: tempv(ka)
941  real(rp), intent(out) :: qv_u(ka)
942  real(rp), intent(out) :: cape
943  real(rp), intent(out) :: cloudtop
944  real(rp), intent(out) :: qc(ka)
945  real(rp), intent(out) :: qi(ka)
946  real(rp), intent(out) :: qvdet(ka)
947  real(rp), intent(out) :: qcdet(ka)
948  real(rp), intent(out) :: qidet(ka)
949  real(rp), intent(out) :: flux_qr(ka)
950  real(rp), intent(out) :: flux_qs(ka)
951  real(rp), intent(out) :: theta_eu(ka)
952  real(rp), intent(out) :: theta_ee(ka)
953  integer, intent(out) :: i_convflag
957  integer, intent(out) :: k_lcl
958  integer, intent(out) :: k_top
959  integer, intent(out) :: k_ml
960  real(rp), intent(out) :: zlcl
961  integer, intent(out) :: k_lc, k_let, k_pbl
962  real(rp), intent(out) :: zmix
963  real(rp), intent(out) :: presmix
964  real(rp), intent(out) :: umfnewdold(ka)
965  real(rp), intent(out) :: dpthmx
966  !---------------------------------------------------------------------------
967 
968  integer, parameter :: itr_max = 10000
969 
970  integer :: kk,nk
971  integer :: k_llfc
973  integer :: n_uslcheck
974  integer :: k_lclm1
975  integer :: k_start
976  integer :: k_check(ka)
977  integer :: n_check
978  integer :: n_layers
979  integer :: nchm
980  integer :: nuchm
981  integer :: nnn
982  integer :: itr
983 
984  real(rp) :: cloudhight
985  real(rp) :: qrout(ka)
986  real(rp) :: qsout(ka)
987  real(rp) :: pres300
988  real(rp) :: pres15
990  real(rp) :: theta_mix
991  real(rp) :: temp_mix
992  real(rp) :: qv_mix
993  real(rp) :: emix
994  real(rp) :: temp_dew
995  real(rp) :: templog
996  real(rp) :: deltaz
997  real(rp) :: temp_env
998  real(rp) :: tempv_env
999  real(rp) :: qvenv
1000  real(rp) :: w_cz
1001  real(rp) :: w_g
1002  real(rp) :: w_lcl
1003  real(rp) :: temp_lcl
1004  real(rp) :: tempv_lcl
1005  real(rp) :: pres_lcl
1006  real(rp) :: dtvv
1007  real(rp) :: dtrh
1008  real(rp) :: radius
1009  real(rp) :: dptt
1010  real(rp) :: dumfdp
1011  real(rp) :: d_min
1012  real(rp) :: umfnew,umfold
1013  real(rp) :: chmax
1014  real(rp) :: cldhgt(ka)
1015  real(rp) :: dpthmin
1016  real(rp) :: rh_lcl
1017  real(rp) :: u00
1018  real(rp) :: dqssdt
1019  real(rp) :: qs_lcl
1020  !real(RP) :: tempvq_u(KA)
1021  ! -----
1022  pres300 = pres(ks) - depth_usl*100._rp ! pressure @ surface - 300 mb. maybe 700mb or so default depth_usl is 300hPa
1023  do kk = ks, ke
1024  tempv(kk) = temp(kk) * ( 1.0_rp + epstvap * qv(kk) ) ! vertual temperature
1025  end do
1026 
1027  ! search above 300 hPa index to "k_llfc"
1028  do kk = ks, ke
1029  if (pres(kk) >= pres300) k_llfc = kk
1030  end do
1031  ! usl(updraft sourcer layer) has interval for 15mb
1032  n_check = ks ! first layer
1033  k_check(ks) = ks ! first layer
1034  pres15 = pres(ks) - 15.e2_rp
1035 ! k_check = KS
1036  ! calc 15 hpa interval Num of layer(n_check) and index of layer(k_check)
1037  do kk = ks+1, k_llfc
1038  if(pres(kk) < pres15 ) then
1039  n_check = n_check + 1
1040  k_check(n_check) = kk
1041  pres15 = pres15 - 15.e2_rp
1042  end if
1043  end do
1044 
1045  ! main routine
1046  dpthmin = 50.e2_rp ! set initialy USL layer depth (50 hPa)
1047  i_convflag = 2 ! no convection set initialy
1048  nuchm = ks-1 ! initialize (used shallow convection check) !like "0"
1049  n_uslcheck = ks-1 ! initialize index of usl check like "0"
1050  nchm = ks-1 ! initializelike "0"
1051  cldhgt(:) = 0._rp ! cloud hight is all 0
1052  do itr = 1, itr_max ! usl
1053  n_uslcheck = n_uslcheck + 1 ! nu in WRF
1054  ! calc shallow convection after all layer checked
1055  ! NOTE: This 'if block' is only used
1056  if (n_uslcheck > n_check) then ! if over potentially usl layer
1057  if (i_convflag == 1) then ! if schallow convection then
1058  chmax = 0._rp ! initialize Cloud Height Max
1059  nchm = ks-1 ! initialize index of Cloud Height Max "like 0"
1060  do kk = ks, n_check
1061  nnn = k_check(kk) ! index of checked layer (15 hpa interval layer)
1062  if (cldhgt(nnn) > chmax) then ! get max cloud height in shallow convection
1063  nchm = nnn
1064  nuchm = kk
1065  chmax = cldhgt(nnn)
1066  end if
1067  end do
1068  n_uslcheck = nuchm - 1
1069  cycle ! usl ! recalc updraft for shallow convection
1070  else ! not shallow convection then
1071  i_convflag = 2 ! no convection
1072  return
1073  end if ! end if schallow convection then
1074  end if ! end if over potentially usl layer
1075  k_lc = k_check(n_uslcheck)
1076  n_layers = 0 ! number of USL layer contain discretization layers
1077  dpthmx = 0._rp ! depth of USL (pressure [Pa])
1078  ! nk = k_lc - KS !< check k_lc -KS is not surface or bottom of ground
1079  nk = k_lc - 1
1080  ! calculate above k_lc layers depth of pressure
1081  ! usl layer depth is nessesary 50mb or more
1082  if ( nk + 1 < ks ) then
1083  if(kf_log) then
1084  log_info("CP_kf_trigger",*) 'would go off bottom: cu_kf',pres(ks),k_lc,nk+1,n_uslcheck !,nk i,j
1085  end if
1086  else
1087  do ! serach USL layer index. USL layer is nessesally more 50hPa
1088  nk = nk +1
1089  if ( nk > ke ) then
1090  if(kf_log) then
1091  log_info("CP_kf_trigger",*) 'would go off top: cu_kf'!,nk i,j
1092  end if
1093  exit
1094  end if
1095  dpthmx = dpthmx + deltap(nk) ! depth of pressure
1096  n_layers = n_layers + 1
1097  if (dpthmx > dpthmin) then
1098  exit
1099  end if
1100  end do
1101  end if
1102  if (dpthmx < dpthmin) then ! if dpthmx(USL layer depth) < dptmin(minimum USL layerdepth)
1103  i_convflag = 2
1104  return ! chenge at 3d version
1105  end if
1106  k_pbl = k_lc + n_layers -1
1107  ! initialize
1108  theta_mix = 0._rp; temp_mix = 0._rp; qv_mix = 0._rp; zmix = 0._rp;presmix = 0._rp
1109  ! clc mix variable (USL layer average valuable)
1110  do kk = k_lc, k_pbl
1111  temp_mix = temp_mix + deltap(kk)*temp(kk)
1112  qv_mix = qv_mix + deltap(kk)*qv(kk)
1113  zmix = zmix + deltap(kk)*z_kf(kk)
1114  presmix = presmix + deltap(kk)*pres(kk)
1115  end do
1116  ! mix tmp calculate
1117  temp_mix = temp_mix/dpthmx
1118  qv_mix = qv_mix/dpthmx
1119  presmix = presmix/dpthmx
1120  zmix = zmix/dpthmx
1121  emix = qv_mix * presmix / ( epsvap + qv_mix ) ! water vapor pressure
1122  ! calculate dewpoint and LCL temperature not use look up table
1123  ! LCL is Lifted condensation level
1124  ! this dewpoint formulation is Bolton's formuration (see Emanuel 1994 4.6.2)
1125  templog = log(emix/aliq)
1126  temp_dew = (cliq - dliq*templog)/(bliq - templog) ! dew point temperature
1127  ! temperature @ lcl setting need dewpoit
1128  ! this algolizm proposed by Davies-Jones (1983)
1129  temp_lcl = temp_dew - (0.212_rp + 1.571e-3_rp*(temp_dew - tem00) &
1130  - 4.36e-4_rp*(temp_mix - tem00))*(temp_mix - temp_dew) ! LCL temperature
1131  temp_lcl = min(temp_lcl,temp_mix)
1132  tempv_lcl = temp_lcl * ( 1.0_rp + epstvap * qv_mix ) ! LCL vertual temperature
1133  zlcl = zmix + (temp_lcl - temp_mix)/gdcp ! height of LCL
1134  ! index of z@lcl
1135  ! z(k_lclm1) < zlcl < z(k_lcl)
1136  do kk = k_lc, ke
1137  k_lcl = kk
1138  if( zlcl <= z_kf(kk) ) exit
1139  end do
1140  k_lclm1 = k_lcl - 1
1141  if (zlcl > z_kf(ke)) then
1142  i_convflag = 2 ! no convection
1143  return
1144  end if
1145  !! interpolate environment temperature and vapor at LCL
1146  deltaz = ( zlcl - z_kf(k_lclm1) )/( z_kf(k_lcl)- z_kf(k_lclm1 ) )
1147  temp_env = temp(k_lclm1) + ( temp(k_lcl) - temp(k_lclm1) )*deltaz
1148  qvenv = qv(k_lclm1) + ( qv(k_lcl) - qv(k_lclm1) )*deltaz
1149  tempv_env = temp_env * ( 1.0_rp + epstvap * qvenv )
1150  ! w_cz set review Kain(2004) eq.2
1151  ! w_g set
1152  ! dtlcl setting
1153  ! wg is
1154  if (zlcl < 2.e3_rp) then ! Kain 2004 eq.2
1155 ! w_cz = 0.02_RP*zlcl/2.e3_RP!
1156  w_cz = 1.e-5_rp*zlcl !
1157  else
1158  w_cz = 0.02_rp
1159  end if
1160  !! wg is iapproximate running mean grid resolved vertical velocity at the lcl(m/s)
1161  !! need w0avg
1162  !!wg = wg-c(z)
1163  w_g = (w0avg(k_lclm1) + (w0avg(k_lcl) - w0avg(k_lclm1))*deltaz )*deltax/25.e3_rp - w_cz ! need w0avg
1164  if ( w_g < 1.e-4_rp) then ! too small wg
1165  dtvv = 0._rp
1166  else
1167  dtvv = 4.64_rp*w_g**0.33_rp ! Kain (2004) eq.1 **1/3
1168  end if
1169  !! triggers will be make
1170  dtrh = 0._rp
1171  !! in WRF has 3 type of trigger function
1172  if ( trigger_type == 2 ) then ! new function none
1173  else if ( trigger_type == 3 ) then ! NO2007 trigger function
1174  !...for ETA model, give parcel an extra temperature perturbation based
1175  !...the threshold RH for condensation (U00)...
1176  ! as described in Narita and Ohmori (2007), 12th Mesoscale Conf.
1177  !
1178  !...for now, just assume U00=0.75...
1179  !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
1180  u00 = 0.75_rp
1181  if(u00 < 1._rp) then
1182  qs_lcl=qes(k_lclm1)+(qes(k_lcl)-qes(k_lclm1))*deltaz
1183  rh_lcl = qvenv/qs_lcl
1184  dqssdt = qv_mix*(cliq-bliq*dliq)/((temp_lcl-dliq)*(temp_lcl-dliq))
1185  if(rh_lcl >= 0.75_rp .and. rh_lcl <=0.95_rp)then
1186  dtrh = 0.25_rp*(rh_lcl-0.75_rp)*qv_mix/dqssdt
1187  elseif(rh_lcl > 0.95_rp) then
1188  dtrh = (1._rp/rh_lcl-1._rp)*qv_mix/dqssdt
1189  else
1190  dtrh = 0._rp
1191  end if
1192  end if
1193  end if
1194 
1195  ! check...
1196  if (temp_lcl + dtvv + dtrh < temp_env) then ! kf triggerfucn dtrh is used @ NHM trigger func
1197  ! parcel is not bouyant
1198  ! cycle and check one more up layer(15 hPa )
1199  cycle
1200  else ! parcel is bouyant ,determin updraft
1201  ! theta_lclm1 is need
1202 
1203  ! calc updraft theta_E
1204  call cp_kf_envirtht( presmix, temp_mix, qv_mix, & ! [IN]
1205  theta_eu(k_lclm1) ) ! [OUT]
1206  !
1207  if (dtvv + dtrh > 1.e-4_rp ) then
1208  w_lcl = 1._rp + 0.5_rp*sqrt(2._rp*grav*(dtvv + dtrh)*500._rp/tempv_env)! Kain(2004) eq. 3??
1209  w_lcl = min(w_lcl,3._rp)
1210  else
1211  w_lcl = 1._rp
1212  end if
1213  k_let = k_lcl ! add k_let
1214  ! interpolate pressure
1215  pres_lcl = pres(k_lclm1) + (pres(k_lcl) - pres(k_lclm1))*deltaz
1216  ! denslcl = pres_lcl/(R*tempv_lcl)
1217  ! temp_lcl is already calculated
1218  if (w_g < 0._rp ) then
1219  radius = 1000._rp
1220  elseif (w_g > 0.1_rp) then
1221  radius = 2000._rp
1222  else
1223  radius = 1000._rp + 1000._rp*w_g*10._rp ! wg/0.1 -> w_g*10
1224  end if
1225 
1226  call cp_kf_updraft ( &
1227  ka, ks, ke, & ! [IN]
1228  k_lcl, k_pbl, dz_kf(:), z_kf(:), & ! [IN]
1229  w_lcl, temp_lcl, tempv_lcl, pres_lcl, & ! [IN]
1230  qv_mix, qv(:), temp(:), tempv_env, & ! [IN]
1231  zlcl, pres(:), deltap(:), & ! [IN]
1232  deltax, radius, dpthmx, & ! [IN]
1233  k_let, theta_eu(:), & ! [INOUT]
1234  k_top, & ! [OUT]
1235  umf(:), umflcl, & ! [OUT]
1236  upent(:), updet(:), & ! [OUT]
1237  umfnewdold(:), umfnew, umfold, & ! [OUT]
1238  temp_u(:), theta_ee(:), & ! [OUT]
1239  cloudhight, cloudtop, & ! [OUT]
1240  qv_u(:), qc(:), qi(:), qrout(:), qsout(:), & ! [OUT]
1241  qvdet(:), qcdet(:), qidet(:), & ! [OUT]
1242  cape, & ! [OUT]
1243  flux_qr(:), flux_qs(:) ) ! [OUT]
1244 
1245  ! temporaly cloud height for shallow convection
1246  cldhgt(k_lc) = cloudhight
1247  ! minimum cloud height for deep convection
1248  ! Kain (2004) eq.7
1249  if(temp_lcl > 293._rp) then
1250  d_min = 4.e3_rp
1251  elseif (temp_lcl < 293._rp .and. temp_lcl > 273._rp) then
1252  d_min = 2.e3_rp + 100._rp*(temp_lcl - tem00)
1253  else
1254  d_min = 2.e3_rp
1255  end if
1256 
1257  ! convection type check
1258  if ( k_top <= k_lcl .or. &
1259  k_top <= k_pbl .or. &
1260  k_let+1 <= k_pbl ) then ! no convection allowd
1261  cloudhight = 0._rp
1262  cldhgt(k_lc) = cloudhight ! 0
1263  do kk = k_lclm1,k_top
1264  umf(kk) = 0._rp
1265  upent(kk) = 0._rp
1266  updet(kk) = 0._rp
1267  qcdet(kk) = 0._rp
1268  qidet(kk) = 0._rp
1269  flux_qr(kk) = 0._rp
1270  flux_qs(kk) = 0._rp
1271  end do
1272  elseif (cloudhight > d_min .and. cape > 1._rp) then ! deep convection
1273  i_convflag = 0 ! deep convection
1274  exit ! exit usl loop
1275  else ! shallow convection
1280  i_convflag = 1
1281  if(n_uslcheck == nuchm) then !
1282  exit ! exit usl loop
1283  else
1284  do kk = k_lclm1,k_top
1285  umf(kk) = 0._rp
1286  upent(kk) = 0._rp
1287  updet(kk) = 0._rp
1288  qcdet(kk) = 0._rp
1289  qidet(kk) = 0._rp
1290  flux_qr(kk) = 0._rp
1291  flux_qs(kk) = 0._rp
1292  end do
1293  end if ! if(n_uslcheck == NUCHM) then
1294  end if ! convection type
1295  end if ! triggeer
1296  end do ! usl
1297  if ( itr .ge. itr_max ) then
1298  log_error("CP_kf_trigger",*) 'iteration max count was reached in the USL loop in the KF scheme'
1299  call prc_abort
1300  end if
1301 
1302 
1303  if (i_convflag == 1) then ! shallow convection
1304  k_start = max(k_pbl,k_lcl)
1305  k_let = k_start
1306  end if
1307  !
1308  ! IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1309  ! THIS LEVEL.
1310  if (k_let == k_top) then
1311  updet(k_top) = umf(k_top) + updet(k_top) - upent(k_top)
1312  qcdet(k_top) = qc(k_top)*updet(k_top)*umfnew/umfold
1313  qidet(k_top) = qi(k_top)*updet(k_top)*umfnew/umfold
1314  upent(k_top) = 0._rp
1315  umf(k_top) = 0._rp
1316  else
1317  !! begin total detrainment at the level above the k_let
1318  dptt = 0._rp
1319  do kk = k_let+1,k_top
1320  dptt = dptt + deltap(kk)
1321  end do
1322  dumfdp = umf(k_let)/dptt
1323  !!
1324 
1331  do kk = k_let+1,k_top
1332  if(kk == k_top) then
1333  updet(kk) = umf(kk-1)
1334  upent(kk) = 0._rp
1335  qcdet(kk) = updet(kk)*qc(kk)*umfnewdold(kk)
1336  qidet(kk) = updet(kk)*qi(kk)*umfnewdold(kk)
1337  else
1338  umf(kk) = umf(kk-1) - deltap(kk)*dumfdp
1339  upent(kk) = umf(kk)*(1._rp - 1._rp/umfnewdold(kk))
1340  updet(kk) = umf(kk-1) - umf(kk) + upent(kk)
1341  qcdet(kk) = updet(kk)*qc(kk)*umfnewdold(kk)
1342  qidet(kk) = updet(kk)*qi(kk)*umfnewdold(kk)
1343  end if
1344  if (kk >= k_let+2) then
1345  flux_qr(kk) = umf(kk-1)*qrout(kk)
1346  flux_qs(kk) = umf(kk-1)*qsout(kk)
1347  end if
1348  !
1349  end do
1350 
1351  end if ! let== ktop
1352 
1353 
1355  do kk = ks,k_top
1356  if(temp(kk) > tem00) k_ml = kk !! melt layer
1357  end do
1358  !
1359  do kk = ks,k_lclm1
1360  !!
1361  if(kk >= k_lc) then
1362  !
1363  if (kk == k_lc) then ! if bototm of USL(pbl)
1364  upent(kk) = umflcl*deltap(kk)/dpthmx
1365  umf(kk) = umflcl*deltap(kk)/dpthmx
1366  elseif (kk <= k_pbl) then ! if in USL(pbl)
1367  upent(kk) = umflcl*deltap(kk)/dpthmx
1368  umf(kk) = umf(kk-1) + upent(kk) !! assume no detrain
1369  else
1370  upent(kk) = 0._rp
1371  umf(kk) = umflcl
1372  end if
1373  !
1374  temp_u(kk) = temp_mix + (z_kf(kk) - zmix)*gdcp ! layner assumption of temp_u
1375  qv_u(kk) = qv_mix
1376  ! wu(kk) = w_lcl no use @wrf
1377  else ! below USL layer no updraft
1378  temp_u(kk) = 0._rp
1379  qv_u(kk) = 0._rp
1380  umf(kk) = 0._rp
1381  upent(kk) = 0._rp
1382  end if
1383  !
1384  updet(kk) = 0._rp
1385  qvdet(kk) = 0._rp
1386  qc(kk) = 0._rp
1387  qi(kk) = 0._rp
1388  qrout(kk) = 0._rp
1389  qsout(kk) = 0._rp
1390  flux_qr(kk) = 0._rp
1391  flux_qs(kk) = 0._rp
1392  qcdet(kk) = 0._rp
1393  qidet(kk) = 0._rp
1394  !!calc theta_E environment
1395  call cp_kf_envirtht( pres(kk), temp(kk), qv(kk), & ! [IN]
1396  theta_ee(kk) ) ! [OUT]
1397  !!
1398  end do
1399 
1400  ! define variables above cloud top
1401  do kk = k_top+1,ke
1402  umf(kk) = 0._rp
1403  upent(kk) = 0._rp
1404  updet(kk) = 0._rp
1405  qvdet(kk) = 0._rp
1406  qc(kk) = 0._rp
1407  qi(kk) = 0._rp
1408  qrout(kk) = 0._rp
1409  qsout(kk) = 0._rp
1410  flux_qr(kk) = 0._rp
1411  flux_qs(kk) = 0._rp
1412  qcdet(kk) = 0._rp
1413  qidet(kk) = 0._rp
1414  end do
1415  do kk = k_top+2,ke
1416  temp_u(kk) = 0._rp
1417  qv_u(kk) = 0._rp
1418  end do
1419 
1420  return
1421  end subroutine cp_kf_trigger
1422 
1423  !------------------------------------------------------------------------------
1427 !OCL SERIAL
1428  subroutine cp_kf_updraft ( &
1429  KA, KS, KE, &
1430  k_lcl, k_pbl, dz_kf, z_kf, &
1431  w_lcl, temp_lcl, tempv_lcl, pres_lcl, &
1432  qv_mix, qv, temp, tempv_env, &
1433  zlcl, pres, deltap, &
1434  deltax, radius, dpthmx, &
1435  k_let, theta_eu, &
1436  k_top, &
1437  umf, umflcl, &
1438  upent, updet, &
1439  umfnewdold, umfnew,umfold, &
1440  temp_u, theta_ee, &
1441  cloudhight, cloudtop, &
1442  qv_u, qc, qi, qrout, qsout, &
1443  qvdet, qcdet, qidet, &
1444  cape, &
1445  flux_qr, flux_qs )
1446  use scale_const,only :&
1447  grav => const_grav, &
1448  rdry => const_rdry, &
1449  epstvap => const_epstvap
1450  implicit none
1451  integer, intent(in) :: ka, ks, ke
1452  integer, intent(in) :: k_lcl
1453  integer, intent(in) :: k_pbl
1454  real(rp), intent(in) :: dz_kf(ka), z_kf(ka)
1455  real(rp), intent(in) :: w_lcl
1456  real(rp), intent(in) :: temp_lcl
1457  real(rp), intent(in) :: tempv_lcl
1458  real(rp), intent(in) :: pres_lcl
1459  real(rp), intent(in) :: qv_mix
1460  real(rp), intent(in) :: qv(ka)
1461  real(rp), intent(in) :: temp(ka)
1462  real(rp), intent(in) :: tempv_env
1463  real(rp), intent(in) :: zlcl
1464  real(rp), intent(in) :: pres(ka)
1465  real(rp), intent(in) :: deltap(ka)
1466  real(rp), intent(in) :: deltax
1467  real(rp), intent(in) :: radius
1468  real(rp), intent(in) :: dpthmx
1469 
1470  integer, intent(inout) :: k_let
1471  real(rp), intent(inout) :: theta_eu(ka)
1472 
1473  integer, intent(out) :: k_top
1474  real(rp), intent(out) :: umf(ka)
1475  real(rp), intent(out) :: umflcl
1476  real(rp), intent(out) :: upent(ka)
1477  real(rp), intent(out) :: updet(ka)
1478  real(rp), intent(out) :: umfnewdold(ka)
1479  real(rp), intent(out) :: umfnew,umfold
1480  real(rp), intent(out) :: temp_u(ka)
1481  real(rp), intent(out) :: theta_ee(ka)
1482  real(rp), intent(out) :: cloudhight
1483  real(rp), intent(out) :: cloudtop
1484  real(rp), intent(out) :: qv_u(ka)
1485  real(rp), intent(out) :: qc(ka)
1486  real(rp), intent(out) :: qi(ka)
1487  real(rp), intent(out) :: qrout(ka)
1488  real(rp), intent(out) :: qsout(ka)
1489  real(rp), intent(out) :: qvdet(ka)
1490  real(rp), intent(out) :: qcdet(ka)
1491  real(rp), intent(out) :: qidet(ka)
1492  real(rp), intent(out) :: cape
1493  real(rp), intent(out) :: flux_qr(ka)
1494  real(rp), intent(out) :: flux_qs(ka)
1495 
1496  integer :: kk,kkp1
1497  integer :: k_lclm1
1498  real(rp) :: tempv_u(ka)
1499  real(rp) :: tempvq_u(ka)
1500  real(rp) :: denslcl
1501  real(rp) :: ee1,ud1, ee2,ud2
1502  real(rp) :: f_eq(ka)
1503  real(rp) :: f_mix1,f_mix2
1504  real(rp) :: rei,dilbe
1505  real(rp) :: qcnew,qinew
1506  real(rp) :: qfrz
1507  real(rp) :: f_frozen1
1508  real(rp) :: temptmp
1509  real(rp) :: temptmp_ice
1510  real(rp) :: tempv(ka)
1511  real(rp) :: wtw
1512  real(rp) :: boeff
1513  real(rp) :: boterm
1514  real(rp) :: dztmp
1515  real(rp) :: entterm
1516  real(rp) :: theta_tmp
1517  real(rp) :: wu(ka)
1518  real(rp) :: qvtmp, qctmp, qitmp
1519  real(rp) :: temp_u95, temp_u10
1520  real(rp) :: qold
1521  real(rp),parameter :: temp_frzt = 268.16_rp
1522  real(rp),parameter :: temp_frzb = 248.16_rp
1523  ! -----
1524  ! initialize
1525  ! only qv exist other water variables is none
1526  umf(:) = 0._rp
1527  f_eq(:) = 0._rp
1528  upent(:) = 0._rp
1529  updet(:) = 0._rp
1530  temp_u(:) = 0._rp
1531  qv_u(:) = 0._rp
1532  qc(:) = 0._rp
1533  qi(:) = 0._rp
1534  qrout(:) = 0._rp
1535  qsout(:) = 0._rp
1536  qvdet(:) = 0._rp
1537  qcdet(:) = 0._rp
1538  qidet(:) = 0._rp
1539  flux_qr(:) = 0._rp
1540  flux_qs(:) = 0._rp
1541  !
1542  ee1 = 1._rp
1543  ud1 = 0._rp
1544  rei = 0._rp
1545  dilbe = 0._rp
1546  cape = 0._rp
1547  do kk = ks, ke
1548  tempv(kk) = temp(kk) *( 1.0_rp + epstvap * qv(kk) ) ! vertual temperature
1549  end do
1550  ! initial updraft mass flux
1551  umfnewdold(:) = 1._rp
1552  k_lclm1 = k_lcl - 1
1553  wtw = w_lcl*w_lcl
1554  denslcl = pres_lcl/(rdry*tempv_lcl)
1555  umf(k_lclm1) = denslcl*1.e-2_rp*deltax**2 ! (A0)
1556  umflcl = umf(k_lclm1)
1557  umfold = umflcl
1558  umfnew = umfold
1559  ! initial value setup
1560  temp_u(k_lclm1) = temp_lcl
1561  tempv_u(k_lclm1) = tempv_lcl
1562  qv_u(k_lclm1) = qv_mix
1563 
1569  temptmp_ice = temp_frzt
1570 
1573  do kk = k_lclm1,ke-1 ! up_main original(wrf cood K is k_lclm1)
1574  kkp1 = kk + 1
1575  ! temporaly use below layer valuables
1576  f_frozen1 = 0._rp ! frozen rate initialize this variable 0 (water)to 1(ice)
1577  theta_eu(kkp1) = theta_eu(kk) ! up parcel theta_E
1578  qv_u(kkp1) = qv_u(kk) ! up parcel vapor
1579  qc(kkp1) = qc(kk) ! up parcel water
1580  qi(kkp1) = qi(kk) ! up parcel ice
1581 
1585  call cp_kf_tpmix2(pres(kkp1), theta_eu(kkp1), & ! [IN]
1586  qv_u(kkp1), qc(kkp1), qi(kkp1), & ! [INOUT]
1587  temp_u(kkp1), qcnew, qinew ) ! [OUT]
1593  if (temp_u(kkp1) <= temp_frzt ) then ! if frozen temperature then ice is made
1594  if (temp_u(kkp1) > temp_frzb) then
1595  if (temptmp_ice > temp_frzt) temptmp_ice = temp_frzt
1596  !!# liner assumption determin frozen ratio
1597  f_frozen1 = (temptmp_ice - temp_u(kkp1))/(temptmp_ice - temp_frzb)
1598  else
1599  f_frozen1 = 1._rp ! all water is frozen
1600  endif
1601  f_frozen1 = min(1._rp, max(0._rp, f_frozen1)) ! [add] 2017/05/19
1602  temptmp_ice = temp_u(kkp1)
1603  ! calc how much ice is a layer ?
1604  qfrz = (qc(kkp1) + qcnew)*f_frozen1 ! all ice
1605  qinew = qinew + qcnew*f_frozen1 ! new create ice
1606  qcnew = qcnew - qcnew*f_frozen1 ! new create liquit
1607  qi(kkp1) = qi(kkp1) + qc(kkp1)*f_frozen1 ! ice old + convert liquit to ice
1608  qc(kkp1) = qc(kkp1) - qc(kkp1)*f_frozen1 ! liquit old - convet liquit to ice
1609  ! calculate effect of freezing
1610  ! and determin new create frozen
1611  call cp_kf_dtfrznew( pres(kkp1), qfrz, & ! [IN]
1612  temp_u(kkp1), theta_eu(kkp1), qv_u(kkp1), qi(kkp1) ) ! [OUT]
1613  end if
1614  tempv_u(kkp1) = temp_u(kkp1) * ( 1.0_rp + epstvap * qv_u(kkp1) ) ! updraft vertual temperature
1615  ! calc bouyancy term for verticl velocity
1616  if (kk == k_lclm1) then !! lcl layer exist between kk and kk+1 layer then use interporate value
1617  boeff = (tempv_lcl + tempv_u(kkp1))/(tempv_env + tempv(kkp1)) - 1._rp
1618  boterm = 2._rp*(z_kf(kkp1) - zlcl)*grav*boeff/1.5_rp
1619  dztmp = z_kf(kkp1) - zlcl
1620  else
1621  boeff = (tempv_u(kk) + tempv_u(kkp1))/(tempv(kk) + tempv(kkp1)) - 1._rp
1622  boterm = 2._rp*(dz_kf(kk) )*grav*boeff/1.5_rp
1623  dztmp = dz_kf(kk)
1624  end if
1625  ! entrainment term
1626  entterm = 2._rp*rei*wtw/umfold
1627  ! calc precp(qr) , snow(qr) and vertical virocity
1628  call cp_kf_precipitation( &
1629  grav, dztmp, boterm, entterm, & ! [IN]
1630  wtw, qc(kkp1), qi(kkp1), qcnew, qinew, qrout(kkp1), qsout(kkp1) ) ! [INOUT]
1631 
1636  if (wtw < 1.e-3_rp ) then ! if no vertical velocity
1637  exit ! end calc updraft
1638  else
1639  wu(kkp1) = sqrt(wtw)
1640  end if
1641  !!# calc tehta_e in environment to entrain into updraft
1642  call cp_kf_envirtht( pres(kkp1), temp(kkp1), qv(kkp1), & ! [IN]
1643  theta_ee(kkp1) ) ! [OUT]
1644  ! rei is the rate of environment inflow
1645  rei = umflcl*deltap(kkp1)*0.03_rp/radius !!# Kain 1990 eq.1 ;Kain 2004 eq.5
1646 
1647  ! calc cape
1648  tempvq_u(kkp1) = temp_u(kkp1) * ( 1.0_rp + epstvap * qv_u(kkp1) - qc(kkp1) - qi(kkp1) )
1649  if (kk == k_lclm1) then!! lcl layer exist between kk and kk+1 then use interporate value
1650  dilbe = ((tempv_lcl + tempvq_u(kkp1))/(tempv_env + tempv(kkp1)) - 1._rp)*dztmp
1651  else
1652  dilbe = ((tempvq_u(kk) + tempvq_u(kkp1))/(tempv(kk) + tempv(kkp1)) - 1._rp)*dztmp
1653  end if
1654  if(dilbe > 0._rp) cape = cape + dilbe*grav
1655 
1659  ! calc entrainment/detrainment
1660  if(tempvq_u(kkp1) <= tempv(kkp1)) then ! if entrain and detrain
1661  ! original KF90 no entrainment allow
1662  ee2 = 0.5_rp ! Kain (2004) eq.4
1663  ud2 = 1._rp
1664  f_eq(kkp1) = 0._rp
1665  else
1666  k_let = kkp1
1667  ! determine teh critical mixed fraction of updraft and environmental air ...
1668  ! if few mix moisture air
1669  f_mix1 = 0.95_rp
1670  f_mix2 = 1._rp - f_mix1
1671  theta_tmp = f_mix1*theta_ee(kkp1) + f_mix2*theta_eu(kkp1)
1672  qvtmp = f_mix1*qv(kkp1) + f_mix2*qv_u(kkp1)
1673  qctmp = f_mix2*qc(kkp1)
1674  qitmp = f_mix2*qi(kkp1)
1675  ! need only temptmp because calc bouyancy
1676  call cp_kf_tpmix2( pres(kkp1), theta_tmp, & ! [IN]
1677  qvtmp, qctmp, qitmp, & ! [INOUT]
1678  temptmp, qcnew, qinew ) ! [OUT]
1679  ! qinew and qcnew is damy valuavle(not use )
1680  temp_u95 = temptmp * ( 1.0_rp + epstvap * qvtmp - qctmp - qitmp )
1681  ! TU95 in old coad
1682  if ( temp_u95 > tempv(kkp1)) then ! few mix but bouyant then ! if95
1683  ee2 = 1._rp ! rate of entrain is 1 -> all entrain
1684  ud2 = 0._rp
1685  f_eq(kkp1) = 1._rp
1686  else
1687  f_mix1 = 0.10_rp
1688  f_mix2 = 1._rp - f_mix1
1689  theta_tmp = f_mix1*theta_ee(kkp1) + f_mix2*theta_eu(kkp1)
1690  qvtmp = f_mix1*qv(kkp1) + f_mix2*qv_u(kkp1)
1691  qctmp = f_mix2*qc(kkp1)
1692  qitmp = f_mix2*qi(kkp1)
1693  ! need only temptmp because calc bouyancy
1694  call cp_kf_tpmix2( pres(kkp1), theta_tmp, & ! [IN]
1695  qvtmp, qctmp, qitmp, & ! [INOUT]
1696  temptmp, qcnew, qinew ) ! [OUT]
1697  ! qinew and qcnew is damy valuavle(not use )
1698  temp_u10 = temptmp * (1.0 + epstvap * qvtmp - qctmp - qitmp )
1699  if (abs(temp_u10 - tempvq_u(kkp1)) < 1.e-3_rp ) then !if10%
1700  ee2 = 1._rp ! all entrain
1701  ud2 = 0._rp
1702  f_eq(kkp1) = 1._rp
1703  else
1704  f_eq(kkp1) = (tempv(kkp1) - tempvq_u(kkp1))*f_mix1 &
1705  & /(temp_u10 - tempvq_u(kkp1))
1706  f_eq(kkp1) = max(0._rp,f_eq(kkp1) )
1707  f_eq(kkp1) = min(1._rp,f_eq(kkp1) )
1708  if (f_eq(kkp1) == 1._rp) then ! f_eq
1709  ee2 = 1._rp
1710  ud2 = 0._rp
1711  elseif (f_eq(kkp1) == 0._rp) then
1712  ee2 = 0._rp
1713  ud2 = 1._rp
1714  else
1715 
1717  call cp_kf_prof5( f_eq(kkp1), & ! [IN]
1718  ee2, ud2 ) ! [INOUT]
1719  end if ! end of f_iq
1720  end if ! end of if10%
1721  end if ! end of if95%
1722  end if! end of entrain/detrain
1723  ee2 = max(ee2,0.5_rp)
1724  ud2 = 1.5_rp*ud2
1725  !
1726  upent(kkp1) = 0.5_rp*rei*(ee1 + ee2)
1727  updet(kkp1) = 0.5_rp*rei*(ud1 + ud2)
1728 
1730  if (umf(kk) - updet(kkp1) < 10._rp) then ! why 10.???
1734  if (dilbe > 0._rp) then
1735  cape = cape - dilbe*grav
1736  end if
1737  k_let = kk ! then k_let = k_top
1738  exit
1739  else
1740  ee1 = ee2
1741  ud1 = ud2
1742  umfold = umf(kk) - updet(kkp1)
1743  umfnew = umfold + upent(kkp1)
1744  umf(kkp1) = umfnew
1745  umfnewdold(kkp1) = umfnew/umfold
1746  qcdet(kkp1) = qc(kkp1)*updet(kkp1)
1747  qidet(kkp1) = qi(kkp1)*updet(kkp1)
1748  qvdet(kkp1) = qv_u(kkp1)
1749  qold = qv_u(kkp1) + qc(kkp1) + qi(kkp1) + qrout(kkp1) + qsout(kkp1)
1750  ! below layer updraft qv and entrain q /new updraft massflux
1751  qv_u(kkp1) = ( umfold*qv_u(kkp1) + upent(kkp1)*qv(kkp1) ) / umfnew
1752  qc(kkp1) = qc(kkp1)*umfold/umfnew
1753  qi(kkp1) = qi(kkp1)*umfold/umfnew
1754  theta_eu(kkp1) = ( umfold * ( 1.0_rp + qold ) * theta_eu(kkp1) &
1755  + upent(kkp1) * ( 1.0_rp + qv(kkp1) ) * theta_ee(kkp1) &
1756  ) / ( umfnew * ( 1.0_rp + qv_u(kkp1) + qc(kkp1) + qi(kkp1) ) )
1757  ! flux_qr is ratio of generation of liquid fallout(RAIN)
1758  ! flux_qi is ratio of generation of ice fallout(SNOW)
1759  flux_qr(kkp1) = qrout(kkp1)*umf(kk)
1760  flux_qs(kkp1) = qsout(kkp1)*umf(kk)
1761 
1762  ! if below cloud base then
1763  ! updraft entrainment is umf@LCL*dp/depth
1764  if(kkp1 <= k_pbl ) upent(kkp1) = upent(kkp1) + umflcl*deltap(kkp1)/dpthmx
1765  end if
1766  end do ! this loop is exit at w=0.
1767  k_top = kk ! vertical coodinate index @ w=0
1768 
1769  cloudhight = z_kf(k_top) - zlcl
1770  cloudtop = z_kf(k_top)
1771 
1772  return
1773  end subroutine cp_kf_updraft
1774 
1775  !------------------------------------------------------------------------------
1779 !OCL SERIAL
1780  subroutine cp_kf_downdraft ( &
1781  KA, KS, KE, &
1782  I_convflag, &
1783  k_lcl, k_ml, k_top, k_pbl, k_let, k_lc, &
1784  z_kf, zlcl, &
1785  u, v, rh, qv, pres, &
1786  deltap, deltax, &
1787  ems, &
1788  theta_ee, &
1789  umf, temp_u, &
1790  flux_qr, flux_qs, tempv, &
1791  wspd, dmf, downent, downdet, &
1792  theta_d, qv_d, &
1793  rain_flux, snow_flux, &
1794  prec_engi, &
1795  k_lfs )
1796  use scale_const,only :&
1797  pre00 => const_pre00, &
1798  rdry => const_rdry, &
1799  rvap => const_rvap, &
1800  cpdry => const_cpdry, &
1801  cpvap => const_cpvap, &
1802  emelt => const_emelt , &
1803  epsvap => const_epsvap, &
1804  epstvap => const_epstvap
1805  use scale_atmos_hydrometeor, only: &
1806  lhf, &
1807  cv_water, &
1808  cv_ice
1809  use scale_atmos_saturation ,only :&
1810  atmos_saturation_psat_liq
1811  use scale_prc, only: &
1812  prc_abort
1813  implicit none
1814  integer, intent(in) :: ka, ks, ke
1815  integer, intent(in) :: i_convflag
1816  integer, intent(in) :: k_lcl
1817  integer, intent(in) :: k_top
1818  integer, intent(in) :: k_pbl
1819  integer, intent(in) :: k_let
1820  integer, intent(in) :: k_lc
1821  integer, intent(in) :: k_ml
1822  real(rp), intent(in) :: z_kf(ka)
1823  real(rp), intent(in) :: u(ka)
1824  real(rp), intent(in) :: v(ka)
1825  real(rp), intent(in) :: rh(ka)
1826  real(rp), intent(in) :: qv(ka)
1827  real(rp), intent(in) :: pres(ka)
1828  real(rp), intent(in) :: deltap(ka)
1829  real(rp), intent(in) :: deltax
1830  real(rp), intent(in) :: ems(ka)
1831  real(rp), intent(in) :: zlcl
1832  real(rp), intent(in) :: umf(ka)
1833  real(rp), intent(in) :: temp_u(ka)
1834  real(rp), intent(in) :: flux_qr(ka)
1835  real(rp), intent(in) :: flux_qs(ka)
1836  real(rp), intent(in) :: tempv(ka)
1837  real(rp), intent(in) :: theta_ee(ka)
1838 
1839  real(rp), intent(out) :: wspd(3)
1840  real(rp), intent(out) :: dmf(ka)
1841  real(rp), intent(out) :: downent(ka)
1842  real(rp), intent(out) :: downdet(ka)
1843  real(rp), intent(out) :: theta_d(ka)
1844  real(rp), intent(out) :: qv_d(ka)
1845  real(rp), intent(out) :: rain_flux
1846  real(rp), intent(out) :: snow_flux
1847  real(rp), intent(out) :: prec_engi
1848  integer, intent(out) :: k_lfs
1849 
1850  integer :: kk, kkp1
1851  integer :: k_z5
1852  integer :: k_dstart
1853  integer :: k_lfstmp
1854  integer :: k_ldt
1855  integer :: k_ldb
1856  real(rp) :: shsign
1857  real(rp) :: vws
1858  real(rp) :: pef
1859  real(rp) :: pef_wind
1860  real(rp) :: pef_cloud
1861  real(rp) :: cbh
1862  real(rp) :: rcbh
1863  real(rp) :: dens_d
1864  real(rp) :: rhbar
1865  real(rp) :: rhh
1866  real(rp) :: dssdt
1867  real(rp) :: dtmp
1868  real(rp) :: qsrh
1869  real(rp) :: rl
1870  real(rp) :: t1rh
1871  real(rp) :: dpthtmp
1872  real(rp) :: dpthdet
1873  real(rp) :: temp_d(ka)
1874  real(rp) :: tempv_d(ka)
1875  real(rp) :: theta_ed(ka)
1876  real(rp) :: qvsd(ka)
1877  real(rp) :: qvs_tmp
1878  real(rp) :: iexn(ka)
1879  real(rp) :: f_dmf
1880  real(rp) :: dq
1881  real(rp) :: es
1882  real(rp) :: ddinc
1883  real(rp) :: total_prec
1884  real(rp) :: total_rain
1885  real(rp) :: total_snow
1886  logical :: flag_rain
1887  real(rp) :: tder
1888  real(rp) :: qvold
1889  ! -----
1890 
1891  wspd(:) = 0._rp
1892  dmf(:) = 0._rp
1893  downent(:) = 0._rp
1894  downdet(:) = 0._rp
1895  theta_d(:) = 0._rp
1896  temp_d(:) = 0._rp
1897  qv_d(:) = 0._rp
1898  rain_flux = 0._rp
1899  snow_flux = 0._rp
1900  total_rain = 0._rp
1901  total_snow = 0._rp
1902  prec_engi = 0._rp
1903  k_lfs = 0
1904 
1905  ! if no convection
1906  if (i_convflag == 2) return ! if 3d, may be change
1907 
1908  ! determine
1909  do kk = ks, ke
1910  if (pres(kk) >= pres(ks)*0.5_rp) k_z5 = kk
1911  end do
1912  ! calc wind speed at LCL and cloud top and
1913  wspd(1) = sqrt(u(k_lcl)*u(k_lcl) + v(k_lcl)*v(k_lcl))
1914  wspd(2) = sqrt(u(k_z5)*u(k_z5) + v(k_z5)*v(k_z5))
1915  wspd(3) = sqrt(u(k_top)*u(k_top) + v(k_top)*v(k_top))
1920  if(wspd(3) > wspd(1)) then
1921  shsign = 1._rp
1922  else
1923  shsign = -1._rp
1924  end if
1925  vws = (u(k_top) - u(k_lcl))*(u(k_top) - u(k_lcl)) &
1926  + (v(k_top) - v(k_lcl))*(v(k_top) - v(k_lcl))
1927 
1928  vws = 1.e3_rp*shsign*sqrt(vws)/(z_kf(k_top) - z_kf(k_lcl))
1929 
1930  ! wind shear type
1931  pef_wind = 1.591_rp + vws*(-0.639_rp + vws*(9.53e-2_rp - vws*4.96e-3_rp) )
1932  ! 0.2 < pef_wind< 0.9
1933  pef_wind = max(pef_wind,0.2_rp)
1934  pef_wind = min(pef_wind,0.9_rp)
1935 
1936 
1938  cbh = (zlcl - z_kf(ks))*3.281e-3_rp ! convert m-> feet that's Amaerican
1939  if( cbh < 3._rp) then
1940  rcbh = 0.02_rp
1941  else
1942  rcbh = 0.96729352_rp + cbh*(-0.70034167_rp + cbh*(0.162179896_rp + cbh*(- &
1943  1.2569798e-2_rp + cbh*(4.2772e-4_rp - cbh*5.44e-6_rp))))
1944  end if
1945  if(cbh > 25.0_rp) rcbh = 2.4_rp
1946  pef_cloud = 1._rp/(1._rp + rcbh)
1947  pef_cloud = min(pef_cloud,0.9_rp)
1948  ! pef is mean of wind shear type and cloud base heigt type
1949  pef = 0.5_rp*(pef_wind + pef_cloud)
1950 
1951 
1952  total_rain = 0.0_rp
1953  total_snow = 0.0_rp
1954  do kk = k_lcl, k_top
1955  total_rain = total_rain + flux_qr(kk)
1956  total_snow = total_snow + flux_qs(kk)
1957  end do
1958 
1959  tder = 0._rp ! initialize evapolation valuavle
1960  if(i_convflag == 1) then ! shallow convection
1961  k_lfs = ks
1962  else ! deep convection
1963  !! start downdraft about 150 hPa above cloud base
1964  k_dstart = k_pbl + 1 ! index of usl layer top +1
1965  k_lfstmp = k_let - 1
1966  do kk = k_dstart+1, ke
1967  dpthtmp = pres(k_dstart) - pres(kk)
1968  if(dpthtmp > 150.e2_rp) then ! if dpth > 150hpa then
1969  k_lfstmp = kk
1970  exit
1971  end if
1972  end do
1973  k_lfstmp = min(k_lfstmp, k_let - 1) ! k_let is equivalent temperature layer index then
1974  k_lfs = k_lfstmp
1975 
1976 
1978  ! dondraft layer
1979  ! downdraft is saturated
1980  if(pres(k_dstart) - pres(k_lfs) > 50.e2_rp) then ! LFS > 50mb(minimum of downdraft source layer)
1981  theta_ed(k_lfs) = theta_ee(k_lfs)
1982  qv_d(k_lfs) = qv(k_lfs)
1983  ! find wet-bulb temperature and qv
1984  call cp_kf_tpmix2dd( pres(k_lfs), theta_ed(k_lfs), & ! [IN]
1985  temp_d(k_lfs), qvs_tmp ) ! [OUT]
1986  call cp_kf_calciexn( pres(k_lfs), qvs_tmp, & ! [IN]
1987  iexn(k_lfs) ) ! [OUT]
1988  theta_d(k_lfs) = temp_d(k_lfs) * iexn(k_lfs)
1989  ! take a first guess at hte initial downdraft mass flux
1990  tempv_d(k_lfs) = temp_d(k_lfs) * ( 1.0_rp + epstvap * qvs_tmp )
1991  dens_d = pres(k_lfs)/(rdry*tempv_d(k_lfs))
1992  dmf(k_lfs) = -(1._rp - pef)*1.e-2_rp*deltax**2*dens_d ! AU0 = 1.e-2*dx**2
1993  downent(k_lfs) = dmf(k_lfs)
1994  downdet(k_lfs) = 0._rp
1995  rhbar = rh(k_lfs)*deltap(k_lfs)
1996  dpthtmp = deltap(k_lfs)
1997  ! calc downdraft entrainment and (detrainment=0.)
1998  ! and dmf
1999  ! downdraft theta and q
2000  do kk = k_lfs-1,k_dstart,-1
2001  kkp1 = kk + 1
2002  downent(kk) = downent(k_lfs)*ems(kk)/ems(k_lfs)
2003  downdet(kk) = 0._rp
2004  dmf(kk) = dmf(kkp1) + downent(kk)
2005  qvold = qv_d(kk)
2006  qv_d(kk) = ( qv_d(kkp1)*dmf(kkp1) + qv(kk)*downent(kk) ) / dmf(kk)
2007  theta_ed(kk) = ( theta_ed(kkp1) * dmf(kkp1) * ( 1.0_rp + qvold ) &
2008  + theta_ee(kk) * downent(kk) * ( 1.0_rp + qv(kk) ) &
2009  ) / ( dmf(kk) * ( 1.0_rp + qv_d(kk) ) )
2010  dpthtmp = dpthtmp + deltap(kk)
2011  rhbar = rhbar + rh(kk)*deltap(kk) ! rh average@ downdraft layer
2012  end do
2013  rhbar = rhbar/dpthtmp
2014  f_dmf = 2._rp*(1._rp - rhbar) ! Kain 2004 eq.11
2015 
2016  call cp_kf_tpmix2dd( pres(k_dstart), theta_ed(k_dstart), & ! [IN]
2017  temp_d(k_dstart), qvs_tmp ) ! [OUT]
2018  if ( do_prec ) then
2019  ! dmf = -umf(k_lcl) (OK?)
2020  dq = 0.0_rp
2021  do kk = k_lcl, k_top
2022  dq = dq + ( flux_qr(kk) * cv_water + flux_qs(kk) * cv_ice ) * ( temp_d(k_dstart) - temp_u(kk) )
2023  end do
2024  temp_d(k_dstart) = temp_d(k_dstart) - dq / ( umf(k_lcl) * cpdry )
2025 
2026  ! calc melting effect
2027  if ( k_lc < k_ml ) then ! if below melt level layer then
2028  flag_rain = .true.
2029  temp_d(k_dstart) = temp_d(k_dstart) - emelt * total_snow / ( umf(k_lcl) * cpdry )
2030  total_rain = total_rain + total_snow
2031  total_snow = 0.0_rp
2032  else
2033  flag_rain = .false.
2034  end if
2035  end if
2036 
2037  ! use check theis subroutine is this
2038  ! temporary: WRF TYPE equations are used to maintain consistency
2039  ! call ATMOS_SATURATION_psat_liq(es,temp_d(k_dstart)) !saturation vapar pressure
2040  es = aliq*exp((bliq*temp_d(k_dstart)-cliq)/(temp_d(k_dstart)-dliq))
2041 
2042  qvs_tmp = epsvap * es / ( pres(k_dstart) - es )
2043  ! Bolton 1980 pseudoequivalent potential temperature
2044  theta_ed(k_dstart) = temp_d(k_dstart)*(pre00/pres(k_dstart))**( ( rdry + rvap * qvs_tmp ) / ( cpdry + cpvap * qvs_tmp ) ) &
2045  * exp((3374.6525_rp/temp_d(k_dstart)-2.5403_rp)*qvs_tmp*(1._rp + 0.81_rp*qvs_tmp))
2046  k_ldt = min(k_lfs-1, k_dstart-1)
2047  dpthdet = 0._rp
2048  do kk = k_ldt,ks,-1 ! start calc downdraft detrain layer index
2049  !!
2050  dpthdet = dpthdet + deltap(kk)
2051  theta_ed(kk) = theta_ed(k_dstart)
2052  qv_d(kk) = qv_d(k_dstart)
2053  call cp_kf_tpmix2dd( pres(kk), theta_ed(kk), & ! [IN]
2054  temp_d(kk), qvs_tmp ) ! [OUT]
2055  if ( do_prec ) then
2056 
2057  qvsd(kk) = qvs_tmp
2058  ! specify RH decrease of 20%/km indowndraft
2059  rhh = 1._rp - 2.e-4_rp*(z_kf(k_dstart) -z_kf(kk) ) ! 0.2/1000.
2060  !
2061 
2063  if( rhh < 1._rp ) then
2064  !
2065  dssdt = (cliq - bliq*dliq)/((temp_d(kk) - dliq)*(temp_d(kk) - dliq))
2066  rl = xlv0 - xlv1*temp_d(kk)
2067  dtmp = rl * qvs_tmp * ( 1.0_rp - rhh ) / ( cpdry + rl * rhh * qvs_tmp * dssdt )
2068  t1rh = temp_d(kk) + dtmp
2069 
2070  !temporary: WRF TYPE equations are used to maintain consistency
2071  !call ATMOS_SATURATION_psat_liq(es,T1rh) !saturation vapar pressure
2072  es = aliq*exp((bliq*t1rh-cliq)/(t1rh-dliq))
2073 
2074  es = rhh*es
2075  qsrh = epsvap * es / ( pres(kk) - es )
2076  if(qsrh < qv_d(kk) ) then
2077  qsrh = qv_d(kk)
2078  t1rh = temp_d(kk) + (qvs_tmp - qsrh) * rl / cpdry
2079  end if
2080 
2081  temp_d(kk) = t1rh
2082  qvs_tmp = qsrh
2083  qvsd(kk) = qvs_tmp
2084  !
2085  end if
2086  else ! do_prec == .false.
2087  qvsd(kk) = qv_d(kk)
2088  end if
2089 
2090  tempv_d(kk) = temp_d(kk) * ( 1.0_rp + epstvap * qvsd(kk) )
2091  if(tempv_d(kk) > tempv(kk) .or. kk == ks) then
2092  k_ldb = kk
2093  exit
2094  end if
2095 
2096  end do ! end calc downdraft detrain layer index
2097 
2098  if( (pres(k_ldb)-pres(k_lfs)) > 50.e2_rp ) then ! minimum downdraft depth !
2099  do kk = k_ldt,k_ldb,-1
2100  kkp1 = kk + 1
2101  downdet(kk) = -dmf(k_dstart)*deltap(kk)/dpthdet
2102  downent(kk) = 0._rp
2103  dmf(kk) = dmf(kkp1) + downdet(kk)
2104  tder = tder + (qvsd(kk) - qv_d(kk))*downdet(kk)
2105  qv_d(kk) = qvsd(kk)
2106  call cp_kf_calciexn( pres(kk), qv_d(kk), & ! [IN]
2107  iexn(kk) ) ! [OUT]
2108  theta_d(kk) = temp_d(kk)*iexn(kk)
2109  end do
2110  end if
2111  end if ! LFS>50mb
2112  end if ! down devap
2113  !!
2117  if (tder < 1._rp) then ! dmf modify
2118  rain_flux = total_rain
2119  snow_flux = total_snow
2120  tder = 0._rp
2121  k_ldb = k_lfs
2122  do kk = ks, k_top
2123  dmf(kk) = 0._rp
2124  downdet(kk) = 0._rp
2125  downent(kk) = 0._rp
2126  theta_d(kk) = 0._rp
2127  temp_d(kk) = 0._rp
2128  qv_d(kk) = 0._rp
2129  end do
2130  else
2131  ddinc = -f_dmf*umf(k_lcl)/dmf(k_dstart) ! downdraft coef Kain 2004 eq.11
2132  if ( flag_rain ) then ! rain
2133  total_prec = total_rain
2134  else
2135  total_prec = total_snow
2136  end if
2137  if( tder*ddinc > total_prec ) then
2138  ddinc = total_prec / tder ! rate of prcp/evap
2139  end if
2140  tder = tder*ddinc
2141  do kk = k_ldb,k_lfs
2142  dmf(kk) = dmf(kk)*ddinc
2143  downent(kk) = downent(kk)*ddinc
2144  downdet(kk) = downdet(kk)*ddinc
2145  end do
2146  ! precipitation - evapolate water
2147  if ( flag_rain ) then
2148  rain_flux = total_rain - tder
2149  snow_flux = total_snow
2150  else
2151  rain_flux = total_rain
2152  snow_flux = total_snow - tder
2153  end if
2154  if ( total_prec < kf_eps ) then ! to avoid FPE w/ Kessler Precip
2155  pef = 0.0_rp
2156  else
2157  pef = ( rain_flux + snow_flux ) / total_prec
2158  endif
2159  prec_engi = rain_flux * cv_water * temp_d(k_ldb) &
2160  + snow_flux * ( cv_ice * temp_d(k_ldb) - lhf )
2161 
2162  !
2166  ! DO NK=LC,LFS
2167  ! UMF(NK)=UMF(NK)*UPDINC
2168  ! UDR(NK)=UDR(NK)*UPDINC
2169  ! UER(NK)=UER(NK)*UPDINC
2170  ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
2171  ! PPTICE(NK)=PPTICE(NK)*UPDINC
2172  ! DETLQ(NK)=DETLQ(NK)*UPDINC
2173  ! DETIC(NK)=DETIC(NK)*UPDINC
2174  ! ENDDO
2175 
2176  if (k_ldb > ks) then
2177  do kk = ks,k_ldb-1
2178  dmf(kk) = 0._rp
2179  downdet(kk) = 0._rp
2180  downent(kk) = 0._rp
2181  theta_d(kk) = 0._rp
2182  temp_d(kk) = 0._rp
2183  qv_d(kk) = 0._rp
2184  end do
2185  end if
2186  ! no dmf is above LFS layer
2187  do kk = k_lfs+1,ke
2188  dmf(kk) = 0._rp
2189  downdet(kk) = 0._rp
2190  downent(kk) = 0._rp
2191  theta_d(kk) = 0._rp
2192  temp_d(kk) = 0._rp
2193  qv_d(kk) = 0._rp
2194  end do
2195  ! no temperature and qv avave downdraft detrainment startlayer (k_ldt)
2196  do kk = k_ldt+1,k_lfs-1
2197  temp_d(kk) = 0._rp
2198  qv_d(kk) = 0._rp
2199  theta_d(kk) = 0._rp
2200  end do
2201  end if ! dmf modify
2202 
2203  return
2204  end subroutine cp_kf_downdraft
2205 
2206  !------------------------------------------------------------------------------
2210 !OCL SERIAL
2211  subroutine cp_kf_compensational (&
2212  KA, KS, KE, &
2213  k_top, k_lc, k_pbl, k_ml, k_lfs, &
2214  dz_kf, z_kf, pres, deltap, deltax, &
2215  temp_bf, qv, &
2216  ems, emsd, &
2217  presmix, zmix, dpthmx, &
2218  cape, &
2219  temp_u, qvdet, umflcl, &
2220  qc, qi, flux_qr, flux_qs, &
2221  umfnewdold, &
2222  wspd, &
2223  qv_d, theta_d, &
2224  prec_engi, &
2225  KF_DTSEC, &
2226  RHOD, i, j, &
2227  I_convflag, k_lcl_bf, &
2228  umf, upent, updet, &
2229  qcdet, qidet, dmf, downent, downdet, &
2230  rain_flux, snow_flux, &
2231  nic, &
2232  theta_nw, &
2233  qv_g, qc_nw, qi_nw, qr_nw, qs_nw, &
2234  sflx_rain, sflx_snow, sflx_engi, &
2235  cldfrac_KF, timecp, time_advec )
2236  use scale_const,only :&
2237  pre00 => const_pre00, &
2238  grav => const_grav, &
2239  rdry => const_rdry, &
2240  rvap => const_rvap, &
2241  cpdry => const_cpdry, &
2242  cpvap => const_cpvap, &
2243  emelt => const_emelt, &
2244  epsvap => const_epsvap, &
2245  epstvap => const_epstvap
2246  use scale_atmos_saturation, only: &
2247  atmos_saturation_psat_liq
2248  use scale_prc, only: &
2249  prc_abort
2250  implicit none
2251  integer, intent(in) :: ka, ks, ke
2252  integer, intent(in) :: k_top
2253  integer, intent(in) :: k_lc
2254  integer, intent(in) :: k_pbl
2255  integer, intent(in) :: k_ml
2256  integer, intent(in) :: k_lfs
2257  real(rp), intent(in) :: dz_kf(ka),z_kf(ka)
2258  real(rp), intent(in) :: pres(ka)
2259  real(rp), intent(in) :: deltap(ka)
2260  real(rp), intent(in) :: deltax
2261  real(rp), intent(in) :: temp_bf(ka)
2262  real(rp), intent(in) :: qv(ka)
2263  real(rp), intent(in) :: ems(ka)
2264  real(rp), intent(in) :: emsd(ka)
2265  real(rp), intent(in) :: presmix
2266  real(rp), intent(in) :: zmix
2267  real(rp), intent(in) :: dpthmx
2268  real(rp), intent(in) :: cape
2269  real(rp), intent(in) :: temp_u(ka)
2270  real(rp), intent(in) :: qvdet(ka)
2271  real(rp), intent(in) :: umflcl
2272  real(rp), intent(in) :: qc(ka)
2273  real(rp), intent(in) :: qi(ka)
2274  real(rp), intent(in) :: flux_qr(ka)
2275  real(rp), intent(in) :: flux_qs(ka)
2276  real(rp), intent(in) :: umfnewdold(ka)
2277  real(rp), intent(in) :: wspd(3)
2278  real(rp), intent(in) :: qv_d(ka)
2279  real(rp), intent(in) :: theta_d(ka)
2280  real(rp), intent(in) :: prec_engi
2281  real(rp), intent(in) :: kf_dtsec
2282  real(rp), intent(in) :: rhod(ka)
2283  integer, intent(in) :: i, j
2284 
2285  integer, intent(inout) :: i_convflag
2286  integer, intent(inout) :: k_lcl_bf
2287  real(rp), intent(inout) :: umf(ka)
2288  real(rp), intent(inout) :: upent(ka)
2289  real(rp), intent(inout) :: updet(ka)
2290  real(rp), intent(inout) :: qcdet(ka)
2291  real(rp), intent(inout) :: qidet(ka)
2292  real(rp), intent(inout) :: dmf(ka)
2293  real(rp), intent(inout) :: downent(ka)
2294  real(rp), intent(inout) :: downdet(ka)
2295  real(rp), intent(inout) :: rain_flux
2296  real(rp), intent(inout) :: snow_flux
2297 
2298  integer, intent(out) :: nic
2301  real(rp), intent(out) :: theta_nw(ka)
2302  real(rp), intent(out) :: qv_g(ka)
2303  real(rp), intent(out) :: qc_nw(ka)
2304  real(rp), intent(out) :: qi_nw(ka)
2305  real(rp), intent(out) :: qr_nw(ka)
2306  real(rp), intent(out) :: qs_nw(ka)
2307  real(rp), intent(out) :: sflx_rain
2308  real(rp), intent(out) :: sflx_snow
2309  real(rp), intent(out) :: sflx_engi
2310  real(rp), intent(out) :: cldfrac_kf(ka,2)
2311  real(rp), intent(out) :: timecp
2312  real(rp), intent(out) :: time_advec
2313 
2314  integer :: ncount
2315  integer :: ntimecount
2316  integer :: nstep
2317  integer :: noiter
2318  integer :: kk,kkp1
2319  integer :: k_lmax
2320  integer :: k_lcl, k_lclm1
2321 
2322  real(rp) :: umf2(ka), dmf2(ka)
2323  real(rp) :: upent2(ka), updet2(ka)
2324  real(rp) :: qcdet2(ka), qidet2(ka)
2325  real(rp) :: downent2(ka), downdet2(ka)
2326  real(rp) :: rain_flux2
2327  real(rp) :: snow_flux2
2328  real(rp) :: tkemax
2329  real(rp) :: z_lcl
2330  real(rp) :: theta(ka)
2331  real(rp) :: theta_u(ka)
2332  real(rp) :: theta_eu(ka)
2333  real(rp) :: theta_eg(ka)
2334  real(rp) :: iexn(ka)
2335  real(rp) :: qv_env
2336  real(rp) :: qv_mix
2337  real(rp) :: qv_gu(ka)
2338  real(rp) :: temp_g(ka)
2339  real(rp) :: temp_env
2340  real(rp) :: tempv_env
2341  real(rp) :: temp_lcl
2342  real(rp) :: tempv_lcl
2343  real(rp) :: temp_mix
2344  real(rp) :: temp_gu(ka)
2345  real(rp) :: tempv_g(ka)
2346  real(rp) :: tempvq_u(ka)
2347  real(rp) :: es
2348  real(rp) :: qvss
2349  real(rp) :: dq, tdpt, dssdt, emix, rl, tlog
2350  real(rp) :: vconv
2351  real(rp) :: dzz
2352  real(rp) :: deltaz
2353  real(rp) :: dilbe
2354  real(rp) :: theta_g(ka)
2355  real(rp) :: qv_nw(ka)
2356  real(rp) :: dpth
2357  real(rp) :: cape_g
2358  real(rp) :: dcape
2359  real(rp) :: fxm(ka)
2360  real(rp) :: f_cape ,f_capeold
2361  real(rp) :: stab
2362  real(rp) :: dtt_tmp
2363  real(rp) :: dtt
2364  real(rp) :: deltat
2365  real(rp) :: tma,tmb,tmm
2366  real(rp) :: evac
2367  real(rp) :: ainc,ainctmp, aincmx,aincold
2368  real(rp) :: aincfin
2369  real(rp) :: omg(ka)
2370  real(rp) :: topomg
2371  real(rp) :: fbfrc
2372  real(rp) :: dfda
2373  real(rp) :: domg_dp(ka)
2374  real(rp) :: absomgtc,absomg
2375  real(rp) :: f_dp
2376 
2377  real(rp) :: theta_fx(ka)
2378  real(rp) :: qv_fx(ka)
2379  real(rp) :: qc_fx(ka)
2380  real(rp) :: qi_fx(ka)
2381  real(rp) :: qr_fx(ka)
2382  real(rp) :: qs_fx(ka)
2383  real(rp) :: rainfb(ka), snowfb(ka)
2384 
2385  real(rp) :: err
2386  real(rp) :: qinit
2387  real(rp) :: qvfnl
2388  real(rp) :: qhydr
2389  real(rp) :: qpfnl
2390  real(rp) :: qfinl
2391  real(rp) :: cpm
2392  real(rp) :: umf_tmp
2393  real(rp) :: xcldfrac
2394 
2395  real(rp) :: qvold
2396 
2397  integer :: m
2398  ! -----
2399 
2400 
2401  sflx_rain = 0.0_rp
2402  sflx_snow = 0.0_rp
2403  sflx_engi = 0.0_rp
2404 
2405  if(i_convflag == 2) return ! noconvection
2406  k_lcl = k_lcl_bf
2407  if ( do_prec .and. i_convflag == 0 ) then ! deep convection
2408  fbfrc = 0.0_rp
2409  else
2410  fbfrc = 1.0_rp
2411  end if
2412  ! time scale of adjustment
2413  vconv = 0.5_rp*(wspd(1) + wspd(2)) ! k_lcl + k_z5
2414  timecp = deltax/vconv
2415  time_advec = timecp ! advection time sclale (30 minutes < timecp < 60 minutes)
2416  timecp = max(deeplifetime, timecp)
2417  timecp = min(3600._rp, timecp)
2418  if(i_convflag == 1) timecp = shallowlifetime ! shallow convection timescale is 40 minutes
2419  nic = nint(timecp/kf_dtsec)
2420  timecp = nic * kf_dtsec ! determin timecp not change below
2421  !
2422  ! maximam of ainc calculate
2423  aincmx = 1000._rp
2424  k_lmax = max(k_lcl,k_lfs)
2425  do kk = k_lc,k_lmax
2426  if((upent(kk)-downent(kk)) > 1.e-3_rp) then
2427  ainctmp = ems(kk)/( (upent(kk) -downent(kk))*timecp )
2428  aincmx = min(aincmx,ainctmp)
2429  end if
2430  end do
2431  ainc = 1.0_rp
2432  aincfin = 1.0_rp
2433  if(aincmx < ainc ) ainc = aincmx
2434  ! for interpolation save original variable
2435  rain_flux2 = rain_flux
2436  snow_flux2 = snow_flux
2437  do kk = ks,k_top
2438  umf2(kk) = umf(kk)
2439  upent2(kk) = upent(kk)
2440  updet2(kk) = updet(kk)
2441  qcdet2(kk) = qcdet(kk)
2442  qidet2(kk) = qidet(kk)
2443  dmf2(kk) = dmf(kk)
2444  downent2(kk) = downent(kk)
2445  downdet2(kk) = downdet(kk)
2446  end do
2447  f_cape = 1._rp
2448  stab = 1.05_rp - delcape ! default 0.95
2449  noiter = 0 ! noiter=1 then stop iteration
2450  if (i_convflag == 1) then ! shallow convection
2451  ! refarence Kain 2004
2452  tkemax = 5._rp
2453  evac = 0.50_rp*tkemax*1.e-1_rp
2454  ainc = evac*dpthmx*deltax**2/( umflcl*grav*timecp)
2455  rain_flux = rain_flux2*ainc
2456  snow_flux = snow_flux2*ainc
2457  do kk = ks,k_top
2458  umf(kk) = umf2(kk)*ainc
2459  upent(kk) = upent2(kk)*ainc
2460  updet(kk) = updet2(kk)*ainc
2461  qcdet(kk) = qcdet2(kk)*ainc
2462  qidet(kk) = qidet2(kk)*ainc
2463  dmf(kk) = dmf2(kk)*ainc
2464  downent(kk) = downent2(kk)*ainc
2465  downdet(kk) = downdet2(kk)*ainc
2466  end do
2467  aincfin = ainc
2468  end if
2469  ! theta set up by Emanuel Atmospheric convection, 1994 111p
2470  ! original KF theta is calced apploximatly.
2471  do kk = ks,k_top
2472  call cp_kf_calciexn( pres(kk), qv(kk), & ! [IN]
2473  iexn(kk) ) ! [OUT]
2474  theta(kk) = temp_bf(kk) * iexn(kk)
2475  call cp_kf_calciexn( pres(kk), qvdet(kk), & ! [IN]
2476  iexn(kk) ) ! [OUT]
2477  theta_u(kk) = temp_u(kk) * iexn(kk)
2478  end do
2479  temp_g(k_top+1:ke) = temp_bf(k_top+1:ke)
2480  qv_g(k_top+1:ke) = qv(k_top+1:ke)
2481  omg(ks-1) = 0.0_rp
2482 
2483  ! main loop of compensation
2484 #ifdef QUICKDEBUG
2485 ! (ocl for bug of the K compilar)
2486 !OCL NOEVAL
2487 #endif
2488  do ncount = 1,10 ! iteration
2489  dtt = timecp
2490  do kk = ks, k_top-1
2491  domg_dp(kk) = - ( upent(kk) - downent(kk) - updet(kk) - downdet(kk) ) * emsd(kk)
2492  omg(kk) = omg(kk-1) - deltap(kk) * domg_dp(kk) ! downward positive
2493  absomg = abs(omg(kk))
2494  absomgtc = absomg*timecp
2495  f_dp = 0.75_rp*deltap(kk)
2496  if(absomgtc > f_dp)THEN
2497  dtt_tmp = f_dp/absomg
2498  dtt=min(dtt,dtt_tmp) ! it is use determin nstep
2499  end if
2500  end do
2501  !! theta_nw and qv_nw has valus only in 1 to k_top
2502  do kk = ks, k_top
2503  theta_nw(kk) = theta(kk)
2504  qv_nw(kk) = qv(kk)
2505  end do
2506  do kk = ks, k_top-1
2507  fxm(kk) = - omg(kk) * deltax**2 / grav ! mass flux (upward positive)
2508  end do
2509  nstep = nint(timecp/dtt + 1) ! how many time step forwad
2510  deltat = timecp/real(nstep,rp) ! deltat*nstep = timecp
2511  if ( hist_flag ) then
2512  do m = 0, 1
2513  do kk = ks, k_top
2514  hist_work(kk,i,j,i_hist_qv_fc,m) = 0.0_rp
2515  hist_work(kk,i,j,i_hist_qv_ue,m) = 0.0_rp
2516  hist_work(kk,i,j,i_hist_qv_ud,m) = 0.0_rp
2517  end do
2518  end do
2519  ! only for deep convection
2520  do kk = ks, k_top
2521  hist_work(kk,i,j,i_hist_qv_de,0) = 0.0_rp
2522  hist_work(kk,i,j,i_hist_qv_dd,0) = 0.0_rp
2523  end do
2524  end if
2525  theta_fx(ks-1) = 0.0_rp
2526  qv_fx(ks-1) = 0.0_rp
2527  theta_fx(k_top) = 0.0_rp
2528  qv_fx(k_top) = 0.0_rp
2529  do ntimecount = 1, nstep
2533  do kk = ks, k_top-1 ! calc flux variable
2534  if( omg(kk) <= 0.0_rp ) then ! upward
2535  theta_fx(kk) = fxm(kk) * theta_nw(kk) * ( 1.0_rp + qv_nw(kk) )
2536  qv_fx(kk) = fxm(kk) * qv_nw(kk)
2537  else ! downward
2538  theta_fx(kk) = fxm(kk) * theta_nw(kk+1) * ( 1.0_rp + qv_nw(kk+1) )
2539  qv_fx(kk) = fxm(kk) * qv_nw(kk+1)
2540  end if
2541  end do
2542 
2544  do kk = ks, k_top
2545  qvold = qv_nw(kk)
2546  qv_nw(kk) = qv_nw(kk) &
2547  + ( - ( qv_fx(kk) - qv_fx(kk-1) ) &
2548  + updet(kk)*qvdet(kk) + downdet(kk)*qv_d(kk) &
2549  - ( upent(kk) - downent(kk) )*qv(kk) )*deltat*emsd(kk)
2550  theta_nw(kk) = ( theta_nw(kk) * ( 1.0_rp + qvold + qc(kk) + qi(kk) ) &
2551  + ( - ( theta_fx(kk) - theta_fx(kk-1) ) &
2552  + updet(kk) * theta_u(kk) * ( 1.0_rp + qvdet(kk) ) &
2553  + downdet(kk) * theta_d(kk) * ( 1.0_rp + qv_d(kk) ) &
2554  - ( upent(kk) - downent(kk) ) * theta(kk) * ( 1.0_rp + qv(kk) ) &
2555  ) * deltat * emsd(kk) &
2556  ) / ( 1.0_rp + qv_nw(kk) + qc(kk) + qi(kk) )
2557  end do
2558  if ( hist_flag ) then
2559  do kk = ks, k_top
2560  hist_work(kk,i,j,i_hist_qv_fc,i_convflag) = hist_work(kk,i,j,i_hist_qv_fc,i_convflag) &
2561  - ( qv_fx(kk) - qv_fx(kk-1) ) * emsd(kk) * rhod(kk) / nstep
2562  hist_work(kk,i,j,i_hist_qv_ud,i_convflag) = hist_work(kk,i,j,i_hist_qv_ud,i_convflag) &
2563  + updet(kk) * qvdet(kk) * emsd(kk) * rhod(kk) / nstep
2564  hist_work(kk,i,j,i_hist_qv_ue,i_convflag) = hist_work(kk,i,j,i_hist_qv_ue,i_convflag) &
2565  - upent(kk) * qv(kk) * emsd(kk) * rhod(kk) / nstep
2566  end do
2567  if ( i_convflag == 0 ) then ! only for deep convection
2568  do kk = ks, k_top
2569  hist_work(kk,i,j,i_hist_qv_dd,0) = hist_work(kk,i,j,i_hist_qv_dd,0) &
2570  + downdet(kk) * qv_d(kk) * emsd(kk) * rhod(kk) / nstep
2571  hist_work(kk,i,j,i_hist_qv_de,0) = hist_work(kk,i,j,i_hist_qv_de,0) &
2572  + downent(kk) * qv(kk) * emsd(kk) * rhod(kk) / nstep
2573  end do
2574  end if
2575  end if
2576  end do ! ntimecount
2577 
2578  do kk = ks, k_top
2579  theta_g(kk) = theta_nw(kk)
2580  qv_g(kk) = qv_nw(kk)
2581  end do
2582 
2584  if ( qv_g(k_top) < kf_eps ) then
2585  qv_g(k_top-1) = qv_g(k_top-1) + ( qv_g(k_top) - kf_eps ) * ems(k_top) * emsd(k_top-1)
2586  qv_g(k_top) = kf_eps
2587  end if
2588  do kk = k_top-1, ks+1, -1
2589  if( qv_g(kk) < kf_eps ) then ! negative moisture
2590  tma = qv_g(kk+1) * ems(kk+1)
2591  tmb = max( qv_g(kk-1), kf_eps ) * ems(kk-1)
2592  tmm = ( qv_g(kk) - kf_eps ) * ems(kk)
2593  tma = tma * ( 1.0_rp + tmm / ( tma + tmb ) )
2594  tma = max( tma, kf_eps * ems(kk+1) )
2595  qv_g(kk-1) = qv_g(kk-1) + ( qv_g(kk+1) * ems(kk+1) - tma + tmm ) * emsd(kk-1)
2596  qv_g(kk+1) = tma * emsd(kk+1)
2597  qv_g(kk) = kf_eps
2598  end if
2599  end do
2600  if( qv_g(ks) < kf_eps ) then
2601  log_error("CP_kf_compensational",*) "error qv<0 @ Kain-Fritsch cumulus parameterization"
2602  call prc_abort
2603  end if
2604  if ( hist_flag ) then
2605  do kk = ks, k_top
2606  hist_work(kk,i,j,i_hist_qv_nf,i_convflag) = ( qv_g(kk) - qv_nw(kk) ) * rhod(kk) / timecp
2607  end do
2608  end if
2609  ! calculate top layer omega and conpare determ omg
2610  topomg = (updet(k_top) - upent(k_top))*deltap(k_top)*emsd(k_top)
2611  if( abs(topomg - omg(k_top-1)) > 1.e-3_rp) then ! not same omega velocity error
2612  log_error("CP_kf_compensational",*) "KF omega is not consistent",ncount
2613  log_error_cont(*) "omega error",abs(topomg - omg(k_top-1)),k_top,topomg,omg(k_top-1)
2614  call prc_abort
2615  end if
2616  ! convert theta to T
2617  do kk = ks,k_top
2618  call cp_kf_calciexn( pres(kk), qv_g(kk), & ! [IN]
2619  iexn(kk) ) ! [OUT]
2620  temp_g(kk) = theta_g(kk) / iexn(kk)
2621  tempv_g(kk) = temp_g(kk) * ( 1.0_rp + epstvap * qv_g(kk) )
2622  end do
2623 
2624 
2625  if(i_convflag == 1) then
2626  exit ! if shallow convection , calc Cape is not used
2627  end if
2628  temp_mix = 0._rp
2629  qv_mix = 0._rp
2630  do kk = k_lc, k_pbl
2631  temp_mix = temp_mix + deltap(kk)*temp_g(kk)
2632  qv_mix = qv_mix + deltap(kk)*qv_g(kk)
2633  ! presmix = presmix + deltap(kk)
2634  end do
2635  temp_mix = temp_mix/dpthmx
2636  qv_mix = qv_mix/dpthmx
2637 
2641  es = aliq*exp((bliq*temp_mix-cliq)/(temp_mix-dliq))
2642 
2643  qvss = epsvap * es / (presmix -es) ! saturate watervapor
2644 
2645  ! Remove supersaturation for diagnostic purposes, if necessary..
2646  if (qv_mix > qvss) then ! saturate then
2647  rl = xlv0 -xlv1*temp_mix
2648  cpm = cpdry + ( cpvap - cpdry ) * qv_mix
2649  dssdt = qvss*(cliq-bliq*dliq)/( (temp_mix-dliq)**2)
2650  dq = (qv_mix -qvss)/(1._rp + rl*dssdt/cpm)
2651  temp_mix = temp_mix + rl / cpdry * dq
2652  qv_mix = qv_mix - dq
2653  temp_lcl = temp_mix
2654  else ! same as detern trigger
2655  qv_mix = max(qv_mix,0._rp)
2656  emix = qv_mix * presmix / ( epsvap + qv_mix )
2657  tlog = log(emix/aliq)
2658  ! dew point temperature Bolton(1980)
2659  tdpt = (cliq - dliq*tlog)/(bliq - tlog)
2660  temp_lcl = tdpt - (0.212_rp + 1.571e-3_rp*(tdpt - tem00) - 4.36e-4_rp*(temp_mix - tem00))*(temp_mix - tdpt)
2661  temp_lcl = min(temp_lcl,temp_mix)
2662  end if
2663  tempv_lcl = temp_lcl * ( 1.0_rp + epstvap * qv_mix )
2664  z_lcl = zmix + (temp_lcl - temp_mix)/gdcp
2665  do kk = k_lc, ke
2666  k_lcl = kk
2667  if( z_lcl <= z_kf(kk) ) exit
2668  end do
2669  ! estimate environmental tempeature and mixing ratio
2670  ! interpolate environment temperature and vapor at LCL
2671  k_lclm1 = k_lcl - 1
2672  deltaz = ( z_lcl - z_kf(k_lclm1) )/( z_kf(k_lcl)- z_kf(k_lclm1 ) )
2673  temp_env = temp_g(k_lclm1) + ( temp_g(k_lcl) - temp_g(k_lclm1) )*deltaz
2674  qv_env = qv_g(k_lclm1) + ( qv_g(k_lcl) - qv_g(k_lclm1) )*deltaz
2675  tempv_env = temp_env * ( 1.0_rp + epstvap * qv_env )
2676  !! pres_lcl=pres(k_lcl-1)+(pres(k_lcl-1)-pres(k_lcl-1))*deltaz
2677  theta_eu(k_lclm1)=temp_mix*(pre00/presmix)**( ( rdry + rvap * qv_mix ) / ( cpdry + cpvap * qv_mix ) ) &
2678  * exp((3374.6525_rp/temp_lcl-2.5403_rp)*qv_mix*(1._rp + 0.81_rp*qv_mix))
2679 
2680  ! COMPUTE ADJUSTED ABE(ABEG).(CAPE)
2681  cape_g = 0._rp ! cape "_g" add because
2682  do kk=k_lclm1,k_top-1 ! LTOPM1
2683  kkp1=kk+1
2684  theta_eu(kkp1) = theta_eu(kk)
2685  ! get temp_gu and qv_gu
2686  call cp_kf_tpmix2dd( pres(kkp1), theta_eu(kkp1), & ! [IN]
2687  temp_gu(kkp1), qv_gu(kkp1) ) ! [OUT]
2688  tempvq_u(kkp1) = temp_gu(kkp1) * ( 1.0_rp + epstvap * qv_gu(kkp1) - qc(kkp1)- qi(kkp1) )
2689  if(kk == k_lclm1) then ! interporate
2690  dzz = z_kf(k_lcl) - z_lcl
2691  dilbe = ((tempv_lcl + tempvq_u(kkp1))/(tempv_env + tempv_g(kkp1)) - 1._rp)*dzz
2692  else
2693  dzz = dz_kf(kk)
2694  dilbe = ((tempvq_u(kk) + tempvq_u(kkp1))/(tempv_g(kk) + tempv_g(kkp1)) - 1._rp)*dzz
2695  end if
2696  if(dilbe > 0._rp) cape_g = cape_g + dilbe*grav
2697 
2698  ! DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
2699  call cp_kf_envirtht( pres(kkp1), temp_g(kkp1), qv_g(kkp1), & ! [IN]
2700  theta_eg(kkp1) ) ! [OUT]
2701  !! theta_eg(environment theta_E )
2702  !! theta_eu(kkp1) = theta_eu(kkp1)*(1._RP/umfnewdold(kkp1)) + theta_eg(kkp1)*(1._RP - (1._RP/umfnewdold(kkp1)))
2703  theta_eu(kkp1) = theta_eu(kkp1)*(umfnewdold(kkp1)) + theta_eg(kkp1)*(1._rp - (umfnewdold(kkp1)))
2704  end do
2705  if (noiter == 1) exit ! noiteration
2706  dcape = max(cape - cape_g,cape*0.1_rp) ! delta cape
2707  f_cape = cape_g/cape ! ratio of cape new/old
2708  if(f_cape > 1._rp .and. i_convflag == 0) then ! if deep convection and cape is inclease this loop
2709  i_convflag = 2
2710  return
2711  end if
2712  if(ncount /= 1) then
2713  if(abs(ainc - aincold) < 1.e-4_rp) then ! IN cycle not change facter then exit iter next loop
2714  noiter = 1 ! exit this loop in nex step
2715  ainc = aincold
2716  cycle ! iter
2717  end if
2718  dfda = (f_cape - f_capeold)/(ainc - aincold)
2719  if (dfda > 0._rp ) then
2720  noiter = 1 ! exit this loop @next loop step
2721  ainc = aincold
2722  cycle ! iter
2723  end if
2724  end if
2725  aincold = ainc
2726  f_capeold = f_cape
2727  ! loop exit
2728  ! aincmx is upper limit of massflux factor
2729  ! if massflux factor 'ainc' is near "aincmax" then exit
2730  ! but need CAPE is less than 90% of original
2731  if (ainc/aincmx > 0.999_rp .and. f_cape > 1.05_rp-stab ) then
2732  exit
2733  end if
2734  ! loop exit 1. or 2.
2735  ! 1. NEW cape is less than 10% of oliginal cape
2736  ! 2. ncount = 10
2737  if( (f_cape <= 1.05_rp-stab .and. f_cape >= 0.95_rp-stab) .or. ncount == 10) then
2738  exit
2739  else ! no exit
2740  if(ncount > 10) exit ! sayfty ?? ncount musn't over 10...
2741  if(f_cape == 0._rp) then ! f_cape is 0 -> new cape is 0 : too reducted
2742  ainc = ainc*0.5_rp
2743  else
2744  if(dcape < 1.e-4) then ! too small cape then exit at next loop(iter loop) step
2745  noiter = 1
2746  ainc = aincold
2747  cycle
2748  else ! calculate new factor ainc
2749  ainc = ainc*stab*cape/max(dcape,kf_eps)
2750  end if
2751  end if
2752  ainc = min(aincmx,ainc) ! ainc must be less than aincmx
2753  ! if ainc becomes very small, effects of convection will be minimal so just ignore it
2754  if (ainc < 0.05_rp) then ! noconvection
2755  i_convflag = 2
2756  return
2757  end if
2758  !!update valuables use factar ainc
2759  rain_flux = rain_flux2*ainc
2760  snow_flux = snow_flux2*ainc
2761  do kk = ks,k_top
2762  umf(kk) = umf2(kk)*ainc
2763  upent(kk) = upent2(kk)*ainc
2764  updet(kk) = updet2(kk)*ainc
2765  qcdet(kk) = qcdet2(kk)*ainc
2766  qidet(kk) = qidet2(kk)*ainc
2767  dmf(kk) = dmf2(kk)*ainc
2768  downent(kk) = downent2(kk)*ainc
2769  downdet(kk) = downdet2(kk)*ainc
2770  end do
2771  aincfin = ainc
2772  ! go back up for another iteration
2773  end if
2774 
2775  end do ! iter(ncount)
2776 
2777  ! get the cloud fraction
2778  cldfrac_kf(:,:) = 0._rp
2779  if (i_convflag == 1) then
2780  do kk = k_lcl-1, k_top
2781  umf_tmp = umf(kk)/(deltax**2)
2782  xcldfrac = 0.07_rp*log(1._rp+(500._rp*umf_tmp))
2783  xcldfrac = max(1.e-2_rp,xcldfrac)
2784  cldfrac_kf(kk,1) = min(2.e-2_rp,xcldfrac) ! shallow
2785  end do
2786  else
2787  do kk = k_lcl-1, k_top
2788  umf_tmp = umf(kk)/(deltax**2)
2789  xcldfrac = 0.14_rp*log(1._rp+(500._rp*umf_tmp))
2790  xcldfrac = max(1.e-2_rp,xcldfrac)
2791  cldfrac_kf(kk,2) = min(6.e-1_rp,xcldfrac) ! deep
2792  end do
2793  end if
2794 
2796  do kk = ks,k_top
2797  rainfb(kk) = flux_qr(kk) * fbfrc * aincfin
2798  end do
2799  do kk = ks,k_top
2800  snowfb(kk) = flux_qs(kk) * fbfrc * aincfin
2801  end do
2802 
2803  ! no qc qi qr qs inputted in KF scheme
2804  qc_nw(ks:ke) = 0._rp
2805  qi_nw(ks:ke) = 0._rp
2806  qr_nw(ks:ke) = 0._rp
2807  qs_nw(ks:ke) = 0._rp
2808 
2809  if ( hist_flag ) then
2810  do m = 0, 1
2811  do kk = ks, k_top
2812  hist_work(kk,i,j,i_hist_qc_fc,m) = 0.0_rp
2813  hist_work(kk,i,j,i_hist_qc_ud,m) = 0.0_rp
2814  hist_work(kk,i,j,i_hist_qi_fc,m) = 0.0_rp
2815  hist_work(kk,i,j,i_hist_qi_ud,m) = 0.0_rp
2816  hist_work(kk,i,j,i_hist_qr_fc,m) = 0.0_rp
2817  hist_work(kk,i,j,i_hist_qr_rf,m) = 0.0_rp
2818  hist_work(kk,i,j,i_hist_qs_fc,m) = 0.0_rp
2819  hist_work(kk,i,j,i_hist_qs_sf,m) = 0.0_rp
2820  end do
2821  end do
2822  end if
2823  qc_fx(ks-1) = 0.0_rp
2824  qi_fx(ks-1) = 0.0_rp
2825  qr_fx(ks-1) = 0.0_rp
2826  qs_fx(ks-1) = 0.0_rp
2827  qc_fx(k_top) = 0.0_rp
2828  qi_fx(k_top) = 0.0_rp
2829  qr_fx(k_top) = 0.0_rp
2830  qs_fx(k_top) = 0.0_rp
2831  do ntimecount = 1, nstep ! same as T, QV
2832  do kk = ks, k_top-1 ! calc flux variable
2833  if( omg(kk) <= 0.0_rp ) then ! upward
2834  qc_fx(kk) = fxm(kk) * qc_nw(kk)
2835  qi_fx(kk) = fxm(kk) * qi_nw(kk)
2836  qr_fx(kk) = fxm(kk) * qr_nw(kk)
2837  qs_fx(kk) = fxm(kk) * qs_nw(kk)
2838  else ! downward
2839  qc_fx(kk) = fxm(kk) * qc_nw(kk+1)
2840  qi_fx(kk) = fxm(kk) * qi_nw(kk+1)
2841  qr_fx(kk) = fxm(kk) * qr_nw(kk+1)
2842  qs_fx(kk) = fxm(kk) * qs_nw(kk+1)
2843  end if
2844  end do
2845  do kk = ks, k_top
2846  qc_nw(kk) = qc_nw(kk) + ( - ( qc_fx(kk) - qc_fx(kk-1) ) + qcdet(kk) )*deltat*emsd(kk)
2847  qi_nw(kk) = qi_nw(kk) + ( - ( qi_fx(kk) - qi_fx(kk-1) ) + qidet(kk) )*deltat*emsd(kk)
2848  qr_nw(kk) = qr_nw(kk) + ( - ( qr_fx(kk) - qr_fx(kk-1) ) + rainfb(kk) )*deltat*emsd(kk)
2849  qs_nw(kk) = qs_nw(kk) + ( - ( qs_fx(kk) - qs_fx(kk-1) ) + snowfb(kk) )*deltat*emsd(kk)
2850  end do
2851  if ( hist_flag ) then
2852  do kk = ks, k_top
2853  hist_work(kk,i,j,i_hist_qc_fc,i_convflag) = hist_work(kk,i,j,i_hist_qc_fc,i_convflag) &
2854  - ( qc_fx(kk) - qc_fx(kk-1) ) * emsd(kk) * rhod(kk) / nstep
2855  hist_work(kk,i,j,i_hist_qc_ud,i_convflag) = hist_work(kk,i,j,i_hist_qc_ud,i_convflag) &
2856  + qcdet(kk) * emsd(kk) * rhod(kk) / nstep
2857  hist_work(kk,i,j,i_hist_qi_fc,i_convflag) = hist_work(kk,i,j,i_hist_qi_fc,i_convflag) &
2858  - ( qi_fx(kk) - qi_fx(kk-1) ) * emsd(kk) * rhod(kk) / nstep
2859  hist_work(kk,i,j,i_hist_qi_ud,i_convflag) = hist_work(kk,i,j,i_hist_qi_ud,i_convflag) &
2860  + qidet(kk) * emsd(kk) * rhod(kk) / nstep
2861  hist_work(kk,i,j,i_hist_qr_fc,i_convflag) = hist_work(kk,i,j,i_hist_qr_fc,i_convflag) &
2862  - ( qr_fx(kk) - qr_fx(kk-1) ) * emsd(kk) * rhod(kk) / nstep
2863  hist_work(kk,i,j,i_hist_qr_rf,i_convflag) = hist_work(kk,i,j,i_hist_qr_rf,i_convflag) &
2864  + rainfb(kk) * emsd(kk) * rhod(kk) / nstep
2865  hist_work(kk,i,j,i_hist_qs_fc,i_convflag) = hist_work(kk,i,j,i_hist_qs_fc,i_convflag) &
2866  - ( qs_fx(kk) - qs_fx(kk-1) ) * emsd(kk) * rhod(kk) / nstep
2867  hist_work(kk,i,j,i_hist_qs_sf,i_convflag) = hist_work(kk,i,j,i_hist_qs_sf,i_convflag) &
2868  + snowfb(kk) * emsd(kk) * rhod(kk) / nstep
2869  end do
2870  end if
2871 
2872  ! in original qlg= qlpa but it is not nessesary
2873  end do
2874 
2875  ! cumulus parameterization rain is determine
2876  ! if shallow convection then fbfrc = 1. -> noprcpitation
2877  sflx_rain = rain_flux*(1._rp - fbfrc)/(deltax**2)
2878  sflx_snow = snow_flux*(1._rp - fbfrc)/(deltax**2)
2879  sflx_engi = prec_engi*(1._rp - fbfrc)/(deltax**2)
2880 
2881  ! evaluate moisuture budget
2882  qinit = 0._rp ! initial qv
2883  qvfnl = 0._rp ! final qv
2884  qhydr = 0._rp ! final hydrometeor
2885  dpth = 0._rp !
2886  do kk = ks,k_top
2887  dpth = dpth + deltap(kk)
2888  qinit = qinit + qv(kk)*ems(kk)
2889  qvfnl = qvfnl + qv_g(kk)*ems(kk) ! qv
2890  qhydr = qhydr + (qc_nw(kk) + qi_nw(kk) + qr_nw(kk) + qs_nw(kk))*ems(kk)
2891  end do
2892  qpfnl = (rain_flux+snow_flux)*timecp*(1._rp - fbfrc)
2893  qfinl = qvfnl + qhydr + qpfnl
2894  err = ( qfinl - qinit )*100._rp/qinit
2895  if ( abs(err) > 0.05_rp ) then
2896  ! write error message
2897  ! moisture budget error
2898  log_error("CP_kf_compensational",*) "@KF,MOISTURE i,j=",i,j
2899  log_error_cont(*) "--------------------------------------"
2900  log_error_cont('("vert accum rho*qhyd : ",ES20.12)') qhydr
2901  log_error_cont('("vert accum rho*qv : ",ES20.12)') qvfnl-qinit
2902  log_error_cont('("precipitation rate : ",ES20.12)') qpfnl
2903  log_error_cont('("conserv qhyd + qv : ",ES20.12)') qhydr + qpfnl
2904  log_error_cont('("conserv total : ",ES20.12)') qfinl-qinit
2905  log_error_cont(*) "--------------------------------------"
2906  call prc_abort
2907  end if
2908 
2915  if (warmrain) then
2916  !!
2917  !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2918  !!
2919  do kk = ks,ke
2920  theta_nw(kk) = theta_nw(kk) - emelt * ( qi_nw(kk) + qs_nw(kk) ) * rhod(kk) / ( cpdry + ( cpvap - cpdry ) * qv_g(kk) ) * iexn(kk)
2921  qc_nw(kk) = qc_nw(kk) + qi_nw(kk)
2922  qr_nw(kk) = qr_nw(kk) + qs_nw(kk)
2923  qi_nw(kk) = 0.0_rp
2924  qs_nw(kk) = 0.0_rp
2925  end do
2926 
2927  if ( hist_flag ) then
2928  do kk = ks, ke
2929  hist_work(kk,i,j,i_hist_qc_fc,i_convflag) = hist_work(kk,i,j,i_hist_qc_fc,i_convflag) + hist_work(kk,i,j,i_hist_qi_fc,i_convflag)
2930  hist_work(kk,i,j,i_hist_qc_ud,i_convflag) = hist_work(kk,i,j,i_hist_qc_ud,i_convflag) + hist_work(kk,i,j,i_hist_qi_ud,i_convflag)
2931  hist_work(kk,i,j,i_hist_qr_fc,i_convflag) = hist_work(kk,i,j,i_hist_qr_fc,i_convflag) + hist_work(kk,i,j,i_hist_qs_fc,i_convflag)
2932  hist_work(kk,i,j,i_hist_qr_rf,i_convflag) = hist_work(kk,i,j,i_hist_qr_rf,i_convflag) + hist_work(kk,i,j,i_hist_qs_sf,i_convflag)
2933  end do
2934  end if
2935 
2936  return
2937 !!$ elseif( ???? ) then
2938 !!$ ! IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
2939 !!$ ! BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2940 !!$ do kk = KS,KE
2941 !!$ cpm = cp*(1._RP + 0.887*qv_g(kk))
2942 !!$ if(kk < k_ml) then
2943 !!$ temp_g(kk) = temp_g(kk) - (qi_nw(kk) + qs_nw(kk))*EMELT/CPM
2944 !!$ elseif(kk > k_ml) then! kk == k_ml no melt
2945 !!$ temp_g(kk) = temp_g(kk) + (qi_nw(kk) + qs_nw(kk))*EMELT/CPM
2946 !!$ end if
2947 !!$ qc_nw(kk) = qc_nw(kk) + qi_nw(kk)
2948 !!$ qr_nw(kk) = qr_nw(kk) + qs_nw(kk)
2949 !!$ qi_nw(kk) = 0._RP
2950 !!$ qs_nw(kk) = 0._RP
2951 !!$ end do
2952 !!$ return
2953  ! IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
2954  ! OF HYDROMETEORS DIRECTLY.
2955  else
2956 !!$ if( I_QI < 1 ) then
2957 !!$ do kk = KS, KE
2958 !!$ qs_nw(kk) = qs_nw(kk) + qi_nw(kk)
2959 !!$ end do
2960 !!$ end if
2961  return
2962  end if
2963 
2964  return
2965  end subroutine cp_kf_compensational
2966 
2967  !------------------------------------------------------------------------------
2972  subroutine cp_kf_calciexn( pres, qv, iexn )
2973  use scale_const,only :&
2974  pre00 => const_pre00, &
2975  rdry => const_rdry, &
2976  rvap => const_rvap, &
2977  cpdry => const_cpdry, &
2978  cpvap => const_cpvap
2979  implicit none
2980  real(rp),intent(in) :: pres
2981  real(rp),intent(in) :: qv
2982  real(rp),intent(out) :: iexn
2983 
2984  iexn = (pre00/pres)**( ( rdry + rvap * qv ) / ( cpdry + cpvap * qv ) )
2985 
2986  return
2987  end subroutine cp_kf_calciexn
2988 
2989  !------------------------------------------------------------------------------
2998  subroutine cp_kf_precipitation_oc1973( &
2999  G, DZ, BOTERM, ENTERM, &
3000  WTW, QLIQ, QICE, QNEWLQ, QNEWIC, QLQOUT, QICOUT )
3001  implicit none
3002  real(rp), intent(in) :: g
3003  real(rp), intent(in) :: dz
3004  real(rp), intent(in) :: boterm
3005  real(rp), intent(in) :: enterm
3006  real(rp), intent(inout) :: qlqout
3007  real(rp), intent(inout) :: qicout
3008  real(rp), intent(inout) :: wtw
3009  real(rp), intent(inout) :: qliq
3010  real(rp), intent(inout) :: qice
3011  real(rp), intent(inout) :: qnewlq
3012  real(rp), intent(inout) :: qnewic
3013 
3014  real(rp) :: qtot,qnew,qest,g1,wavg,conv,ratio3,oldq,ratio4,dq,pptdrg
3015 
3016  qtot=qliq+qice
3017  qnew=qnewlq+qnewic
3018 
3019 
3022  qest=0.5_rp*(qtot+qnew)
3023  g1=wtw+boterm-enterm-2._rp*g*dz*qest/1.5_rp
3024  IF(g1.LT.0.0)g1=0._rp
3025  wavg=0.5_rp*(sqrt(wtw)+sqrt(g1))
3026  conv=rate*dz/max(wavg,kf_eps) ! KF90 Eq. 9 !wig, 12-Sep-2006: added div by 0 check
3027 
3028 
3032  ratio3=qnewlq/(qnew+1.e-8_rp)
3033  ! OLDQ=QTOT
3034  qtot=qtot+0.6_rp*qnew
3035  oldq=qtot
3036  ratio4=(0.6_rp*qnewlq+qliq)/(qtot+1.e-8_rp)
3037  qtot=qtot*exp(-conv) ! KF90 Eq. 9
3038 
3039 
3041  dq=oldq-qtot
3042  qlqout=ratio4*dq
3043  qicout=(1._rp-ratio4)*dq
3044 
3045 
3047  pptdrg=0.5_rp*(oldq+qtot-0.2_rp*qnew)
3048  wtw=wtw+boterm-enterm-2._rp*g*dz*pptdrg/1.5_rp
3049  IF(abs(wtw).LT.1.e-4_rp)wtw=1.e-4_rp
3050 
3051 
3053  qliq=ratio4*qtot+ratio3*0.4_rp*qnew
3054  qice=(1._rp-ratio4)*qtot+(1._rp-ratio3)*0.4_rp*qnew
3055  qnewlq=0._rp
3056  qnewic=0._rp
3057  return
3058  end subroutine cp_kf_precipitation_oc1973
3059 
3060  !------------------------------------------------------------------------------
3064  subroutine cp_kf_precipitation_kessler( &
3065  G, DZ, BOTERM, ENTERM, &
3066  WTW, QLIQ, QICE, QNEWLQ, QNEWIC, QLQOUT, QICOUT )
3067  implicit none
3068  real(rp), intent(in) :: g
3069  real(rp), intent(in) :: dz
3070  real(rp), intent(in) :: boterm
3071  real(rp), intent(in) :: enterm
3072  real(rp), intent(inout) :: wtw
3073  real(rp), intent(inout) :: qliq
3074  real(rp), intent(inout) :: qice
3075  real(rp), intent(inout) :: qnewlq
3076  real(rp), intent(inout) :: qnewic
3077  real(rp), intent(inout) :: qlqout
3078  real(rp), intent(inout) :: qicout
3079 
3080  real(rp) :: pptdrg
3081  real(rp) :: total_liq, total_ice
3082  ! parameter module value kf_threshold
3083  real(rp) :: auto_qc, auto_qi
3084  auto_qc = kf_threshold
3085  auto_qi = kf_threshold
3086 
3087  total_liq = qliq + qnewlq
3088  total_ice = qice + qnewic
3089 
3090  ! condensate in convective updraft is converted into precipitation
3091  qlqout = max( total_liq - auto_qc, 0.0_rp )
3092  qicout = max( total_ice - auto_qi, 0.0_rp )
3093 
3094  pptdrg = max( 0.5_rp * ( total_liq + total_ice - qlqout - qicout ), 0.0_rp )
3095  wtw=wtw+boterm-enterm-2._rp*g*dz*pptdrg/1.5_rp
3096  IF(abs(wtw).LT.1.e-4_rp)wtw=1.e-4_rp
3097 
3098  qliq = max( total_liq - qlqout, 0.0_rp )
3099  qlqout = total_liq - qliq
3100 
3101  qice = max( total_ice - qicout, 0.0_rp )
3102  qicout = total_ice - qice
3103 
3104  qnewlq=0.0_rp
3105  qnewic=0.0_rp
3106 
3107  return
3108  end subroutine cp_kf_precipitation_kessler
3109 
3110  !------------------------------------------------------------------------------
3113  subroutine cp_kf_tpmix2( p,thes,qu,qliq,qice,tu,qnewlq,qnewic )
3114  use scale_const, only: &
3115  cpdry => const_cpdry, &
3116  cpvap => const_cpvap
3117  implicit none
3118  real(rp), intent(in) :: p, thes
3119  real(rp), intent(inout) :: qu, qliq, qice
3120  real(rp), intent(out) :: tu, qnewlq, qnewic
3121 
3122  real(rp) :: temp,qs,qnew,dq,qtot,rll,cpp
3123 
3124  ! parcel temperature
3125  call cp_kf_tpmix2dd( p, thes, temp, qs )
3126 
3127  dq=qs-qu
3128  IF(dq.LE.0._rp)THEN ! QS <= QU
3129  qnew=qu-qs
3130  qu=qs
3131  ELSE
3132 
3134  qnew=0._rp
3135  qtot=qliq+qice
3136 
3148  IF(qtot.GE.dq)THEN
3149  qliq=qliq-dq*qliq/(qtot+1.e-10_rp)
3150  qice=qice-dq*qice/(qtot+1.e-10_rp)
3151  qu=qs
3152  ELSE
3153  rll=xlv0-xlv1*temp
3154  cpp = cpdry + ( cpvap - cpdry ) * qu
3155  IF(qtot.LT.1.e-10_rp)THEN
3156 
3157  temp=temp+rll*(dq/(1._rp+dq))/cpp
3158  ELSE
3159 
3161  temp=temp+rll*((dq-qtot)/(1._rp+dq-qtot))/cpp
3162  qu=qu+qtot
3163  qliq=0._rp
3164  qice=0._rp
3165  ENDIF
3166  ENDIF
3167  ENDIF
3168  tu=temp
3169  qnewlq=qnew
3170  qnewic=0._rp
3171  return
3172  end subroutine cp_kf_tpmix2
3173 
3174  !------------------------------------------------------------------------------
3178  subroutine cp_kf_dtfrznew( P, QFRZ, TU, THTEU, QU, QICE )
3179  use scale_const, &
3180  pre00 => const_pre00, &
3181  tem00 => const_tem00, &
3182  rdry => const_rdry, &
3183  rvap => const_rvap, &
3184  cpdry => const_cpdry, &
3185  cpvap => const_cpvap, &
3186  epsvap => const_epsvap
3187  use scale_atmos_hydrometeor, only: &
3188  lhs, &
3189  lhf, &
3190  cv_water, &
3191  cv_ice
3192  use scale_atmos_saturation, only :&
3193  atmos_saturation_psat_liq
3194  implicit none
3195  real(rp), intent(in) :: p, qfrz
3196  real(rp), intent(inout) :: tu, thteu, qu, qice
3197 
3198  real(rp) :: rls,rlf,cpp,a,dtfrz,es,qs,dqevap,pii
3205  rls = lhs - ( cv_ice - cpvap ) * tu
3206  rlf = lhf - ( cv_ice - cv_water ) * tu
3207  cpp = cpdry + ( cpvap - cpdry ) * qu
3208 
3209 
3211  a=(cliq-bliq*dliq)/((tu-dliq)*(tu-dliq))
3212  dtfrz = rlf*qfrz/(cpp+rls*qu*a)
3213  tu = tu+dtfrz
3214  ! temporary: WRF TYPE equations are used to maintain consistency
3215  ! call ATMOS_SATURATION_psat_liq(ES,TU) !saturation vapar pressure
3216  es = aliq*exp((bliq*tu-cliq)/(tu-dliq))
3217  qs = es * epsvap / ( p - es )
3218  !
3224  dqevap = max( min(qs-qu, qice), 0.0_rp )
3225  qice = qice-dqevap
3226  qu = qu+dqevap
3227  pii=(pre00/p)**( ( rdry + rvap * qu ) / ( cpdry + cpvap * qu ) )
3228 
3229  thteu = tu*pii*exp((3374.6525_rp/tu - 2.5403_rp)*qu*(1._rp + 0.81_rp*qu))
3230  !
3231  end subroutine cp_kf_dtfrznew
3232 
3233  !------------------------------------------------------------------------------
3245  subroutine cp_kf_prof5( EQ, EE, UD )
3246  implicit none
3247  real(rp), intent(in) :: eq
3248  real(rp), intent(inout) :: ee, ud
3249 
3250  real(rp) :: sqrt2p, a1, a2, a3, p, sigma, fe
3251  real(rp) :: x, y, ey, e45, t1, t2, c1, c2
3252 
3253  DATA sqrt2p,a1,a2,a3,p,sigma,fe/2.506628_rp,0.4361836_rp,-0.1201676_rp, &
3254  0.9372980_rp,0.33267_rp,0.166666667_rp,0.202765151_rp/
3255  x = (eq - 0.5_rp)/sigma
3256  y = 6._rp*eq - 3._rp
3257  ey = exp(y*y/(-2._rp))
3258  e45 = exp(-4.5_rp)
3259  t2 = 1._rp/(1._rp + p*abs(y))
3260  t1 = 0.500498_rp
3261  c1 = a1*t1+a2*t1*t1+a3*t1*t1*t1
3262  c2 = a1*t2+a2*t2*t2+a3*t2*t2*t2
3263  IF(y.GE.0._rp)THEN
3264  ee=sigma*(0.5_rp*(sqrt2p-e45*c1-ey*c2)+sigma*(e45-ey))-e45*eq*eq/2._rp
3265  ud=sigma*(0.5_rp*(ey*c2-e45*c1)+sigma*(e45-ey))-e45*(0.5_rp+eq*eq/2._rp- &
3266  eq)
3267  ELSE
3268  ee=sigma*(0.5_rp*(ey*c2-e45*c1)+sigma*(e45-ey))-e45*eq*eq/2._rp
3269  ud=sigma*(0.5_rp*(sqrt2p-e45*c1-ey*c2)+sigma*(e45-ey))-e45*(0.5_rp+eq* &
3270  eq/2._rp-eq)
3271  ENDIF
3272  ee=ee/fe
3273  ud=ud/fe
3274  end subroutine cp_kf_prof5
3275 
3276  !------------------------------------------------------------------------------
3284  subroutine cp_kf_tpmix2dd( p, thes, ts, qs )
3285  implicit none
3286  real(rp), intent(in) :: p, thes
3287  real(rp), intent(out) :: ts, qs
3288 
3289  real(rp) :: tp,qq,bth,tth,pp,t00,t10,t01,t11,q00,q10,q01,q11
3290  integer :: iptb, ithtb
3291 
3292  ! scaling pressure and tt table index
3293  tp = ( p - plutop ) * rdpr
3294  qq = tp - aint(tp)
3295  iptb = int(tp)+1
3296  if ( iptb < 1 .or. iptb >= kfnp ) then
3297  log_warn("CP_kf_tpmix2dd",*) 'OUT OF BOUNDS (p): ', p, plutop, pbot
3298  iptb = min( max( iptb, 1 ), kfnp-1 )
3299  qq = 0.5_rp + sign(0.5_rp, iptb-2.0_rp) ! qq=0 for iptb==1, qq=1 for iptb==KFNP-1
3300  end if
3301 
3302 
3303  ! scaling the and tt table index
3304  bth = ( the0k(iptb+1) - the0k(iptb) ) * qq + the0k(iptb)
3305  tth = ( thes - bth ) * rdthk
3306  pp = tth - aint(tth)
3307  ithtb = int(tth)+1
3308  if ( ithtb < 1 .or. ithtb >= kfnt ) then
3309  log_warn("CP_kf_tpmix2dd",*) 'OUT OF BOUNDS (thes): ', thes, p, bth, bth + (kfnt-1) / rdthk
3310  ithtb = min( max( ithtb, 1 ), kfnt-1 )
3311  pp = 0.5_rp + sign(0.5_rp, ithtb-2.0_rp) ! pp=0 for ithtb==1, pp=1 for ithtb==KFNT-1
3312  end if
3313 
3314 
3315  t00=ttab(ithtb ,iptb )
3316  t10=ttab(ithtb+1,iptb )
3317  t01=ttab(ithtb ,iptb+1)
3318  t11=ttab(ithtb+1,iptb+1)
3319 
3320  q00=qstab(ithtb ,iptb )
3321  q10=qstab(ithtb+1,iptb )
3322  q01=qstab(ithtb ,iptb+1)
3323  q11=qstab(ithtb+1,iptb+1)
3324 
3325  ! parcel temperature and saturation mixing ratio
3326  ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
3327  qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
3328 
3329  return
3330  end subroutine cp_kf_tpmix2dd
3331 
3332  !------------------------------------------------------------------------------
3342  subroutine cp_kf_envirtht( P1, T1, Q1, THT1 )
3343  use scale_const, only : &
3344  p00 => const_pre00, &
3345  rdry => const_rdry, &
3346  rvap => const_rvap, &
3347  cpdry => const_cpdry, &
3348  cpvap => const_cpvap, &
3349  epsvap => const_epsvap
3350  implicit none
3351  real(rp), intent(in) :: p1, t1, q1
3352  real(rp), intent(out) :: tht1
3353 
3354  real(rp) :: ee,tlog,tdpt,tsat,tht
3355  real(rp),parameter :: c1=3374.6525_rp
3356  real(rp),parameter :: c2=2.5403_rp
3357  ee = q1 * p1 / ( epsvap + q1 )
3358  ! TLOG=ALOG(EE/ALIQ)
3359 
3360  !
3361 !! astrt=1.e-3_RP
3362 !! ainc=0.075_RP
3363 !! a1=ee/aliq
3364 !! tp=(a1-astrt)/ainc
3365 !! indlu=int(tp)+1
3366 !! value=(indlu-1)*ainc+astrt
3367 !! aintrp=(a1-value)/ainc
3368  ! tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
3369  ! change nouse lookuptable
3370  tlog = log(ee/aliq)
3371 
3372  tdpt=(cliq-dliq*tlog)/(bliq-tlog)
3373 
3374  tsat=tdpt - (0.212_rp+1.571e-3_rp*(tdpt-tem00)-4.36e-4_rp*(t1-tem00))*(t1-tdpt)
3375  ! TSAT = 2840._RP/(3.5_RP - log(ee) -4.805_RP) +55
3376 
3377  tht=t1*(p00/p1)**( ( rdry + rvap * q1 ) / ( cpdry + cpvap * q1 ) )
3378  tht1=tht*exp((c1/tsat-c2)*q1*(1._rp+0.81_rp*q1))
3379  !
3380  return
3381  end subroutine cp_kf_envirtht
3382 
3383  !------------------------------------------------------------------------------
3389  subroutine cp_kf_lutab !(SVP1,SVP2,SVP3,SVPT0)
3390  use scale_const, only :&
3391  pre00 => const_pre00, &
3392  grav => const_grav, &
3393  rdry => const_rdry, &
3394  rvap => const_rvap, &
3395  cpdry => const_cpdry, &
3396  cpvap => const_cpvap, &
3397  epsvap => const_epsvap
3398  IMPLICIT NONE
3399  integer :: kp, it, itcnt, i
3400  real(rp) :: dth = 1._rp
3401  real(rp) :: tmin = 150._rp
3402  real(rp) :: toler = 0.001_rp
3403  real(rp) :: dpr, temp, p, es, qs, pi
3404  real(rp) :: thes, tgues, thgues, thtgs
3405  real(rp) :: dt, t1, t0, f0, f1, astrt, ainc, a1
3406 
3407  ! equivalent potential temperature increment: data dth/1._RP/
3408  ! minimum starting temp: data tmin/150._RP/
3409  ! tolerance for accuracy of temperature: data toler/0.001_RP/
3410  ! top pressure (pascals)
3411  plutop=5000.0_rp
3412  ! bottom pressure (pascals)
3413  pbot=110000.0_rp
3414 
3415  ! compute parameters
3416  ! 1._over_(sat. equiv. theta increment)
3417  rdthk=1._rp/dth
3418  ! pressure increment
3419  dpr=(pbot-plutop)/real(kfnp-1)
3420  ! dpr=(pbot-plutop)/REAL(kfnp-1)
3421  ! 1._over_(pressure increment)
3422  rdpr=1._rp/dpr
3423  ! compute the spread of thes
3424  ! thespd=dth*(kfnt-1)
3425 
3426  ! calculate the starting sat. equiv. theta
3427  temp=tmin
3428  p=plutop-dpr
3429  do kp=1,kfnp
3430  p=p+dpr
3431  es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
3432  qs = epsvap * es / ( p - es )
3433  pi = (pre00/p)**( ( rdry + rvap * qs ) / ( cpdry + cpvap * qs ) )
3434  the0k(kp)=temp*pi*exp((3374.6525_rp/temp-2.5403_rp)*qs* &
3435  (1._rp+0.81_rp*qs))
3436  enddo
3437 
3438  ! compute temperatures for each sat. equiv. potential temp.
3439  p=plutop-dpr
3440  do kp=1,kfnp
3441  thes=the0k(kp)-dth
3442  p=p+dpr
3443  do it=1,kfnt
3444  ! define sat. equiv. pot. temp.
3445  thes=thes+dth
3446  ! iterate to find temperature
3447  ! find initial guess
3448  if(it.eq.1) then
3449  tgues=tmin
3450  else
3451  tgues=ttab(it-1,kp)
3452  endif
3453  es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
3454  qs = epsvap * es / ( p - es )
3455  pi = ( pre00 / p )**( ( rdry + rvap * qs ) / ( cpdry + cpvap * qs ) )
3456  thgues=tgues*pi*exp((3374.6525_rp/tgues-2.5403_rp)*qs* &
3457  (1._rp + 0.81_rp*qs))
3458  f0=thgues-thes
3459  t1=tgues-0.5_rp*f0
3460  t0=tgues
3461  itcnt=0
3462  ! iteration loop
3463  do itcnt=1,11
3464  es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
3465  qs = epsvap * es / ( p - es )
3466  pi = ( pre00 / p )**( ( rdry + rvap * qs ) / ( cpdry + cpvap * qs ) )
3467  thtgs=t1*pi*exp((3374.6525_rp/t1-2.5403_rp)*qs*(1._rp + 0.81_rp*qs))
3468  f1=thtgs-thes
3469  if(abs(f1).lt.toler)then
3470  exit
3471  endif
3472  ! itcnt=itcnt+1
3473  dt=f1*(t1-t0)/(f1-f0)
3474  t0=t1
3475  f0=f1
3476  t1=t1-dt
3477  enddo
3478  ttab(it,kp)=t1
3479  qstab(it,kp)=qs
3480  enddo
3481  enddo
3482 
3483  ! lookup table for tlog(emix/aliq)
3484  ! set up intial variable for lookup tables
3485  astrt=1.e-3_rp
3486  ainc=0.075_rp
3487  !
3488  a1=astrt-ainc
3489  do i=1,200
3490  a1=a1+ainc
3491  alu(i)=log(a1)
3492  enddo
3493  !GdCP is g/cp add for SCALE
3494  gdcp = - grav / cpdry ! inital set
3495  return
3496  end subroutine cp_kf_lutab
3497 
3498 end module scale_atmos_phy_cp_kf
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_atmos_hydrometeor::i_hr
integer, parameter, public i_hr
liquid water rain
Definition: scale_atmos_hydrometeor.F90:82
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_atmos_phy_cp_kf::atmos_phy_cp_kf_setup
subroutine, public atmos_phy_cp_kf_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, CZ, AREA, WARMRAIN_in)
Setup initial setup for Kain-Fritsch Cumulus Parameterization.
Definition: scale_atmos_phy_cp_kf.F90:179
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:133
scale_atmos_hydrometeor::i_hs
integer, parameter, public i_hs
snow
Definition: scale_atmos_hydrometeor.F90:84
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:69
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
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_phy_cp_kf::atmos_phy_cp_kf_tendency
subroutine, public atmos_phy_cp_kf_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, U, V, RHOT, TEMP, PRES, QDRY, QV_in, w0avg, FZ, KF_DTSECD, DENS_t, RHOT_t, RHOQV_t, RHOQ_t, nca, SFLX_rain, SFLX_snow, SFLX_engi, cloudtop, cloudbase, cldfrac_dp, cldfrac_sh)
ATMOS_PHY_CP_kf calculate Kain-Fritsch Cumulus Parameterization.
Definition: scale_atmos_phy_cp_kf.F90:503
scale_const::const_emelt
real(rp), parameter, public const_emelt
Definition: scale_const.F90:72
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_phy_cp_kf
module atmosphere / physics / cumulus / Kain-Fritsch
Definition: scale_atmos_phy_cp_kf.F90:45
scale_file_history
module file_history
Definition: scale_file_history.F90:15
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_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:130
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_hydrometeor.F90:83
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:44
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
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_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_hydrometeor.F90:81
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:90
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_hydrometeor::lhf
real(rp), public lhf
latent heat of fusion for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:128
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:650
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
scale_atmos_hydrometeor::lhs
real(rp), public lhs
latent heat of sublimation for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:127
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
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_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_hydrometeor::n_hyd
integer, parameter, public n_hyd
Definition: scale_atmos_hydrometeor.F90:79
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
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:134
scale_atmos_hydrometeor::cp_ice
real(rp), public cp_ice
CP for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:135
Generated by   doxygen 1.8.17