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