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