SCALE-RM
scale_atmos_phy_cp_kf.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
24 #include "inc_openmp.h"
25 module scale_atmos_phy_cp_kf
26  !------------------------------------------------------------------------------
27  !
28  !++ used modules
29  !
30  use scale_precision
31  use scale_stdio
32  use scale_prof
33  use scale_grid_index
34  use scale_const, only: &
35  tem00 => const_tem00
36 
37  use scale_tracer, only: qa
38  !------------------------------------------------------------------------------
39  implicit none
40  private
41 
42  !-----------------------------------------------------------------------------
43  !
44  !++ Public procedure
45  !
46  public :: atmos_phy_cp_kf_setup
47  public :: atmos_phy_cp_kf
48 
49  !-----------------------------------------------------------------------------
50  !++ Public pparameters & variabeles
51  real(RP), public, save :: kf_dt = 5._rp ! kf time scale [min!!!!!]
52 
53  !-----------------------------------------------------------------------------
54  !
55  !++ Private procedures
56  !
57  private :: cp_kf_trigger, cp_kf_updraft, cp_kf_downdraft, cp_kf_compensational
58  private :: calcexn
59  private :: precipitation_oc1973
60  private :: precipitation_kessler
61  private :: tpmix2, dtfrznew, prof5, tpmix2dd, envirtht
62 
63  abstract interface
64  subroutine precipitation( &
65  QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,QNEWLQ, &
66  QNEWIC,QLQOUT,QICOUT,G)
67  use scale_precision
68  real(RP), INTENT(IN ) :: g
69  real(RP), INTENT(IN ) :: dz,boterm,enterm!,RATE to be local variablebles
70  real(RP), INTENT(INOUT) :: qlqout,qicout,wtw,qliq,qice,qnewlq,qnewic
71  end subroutine precipitation
72  end interface
73 
74  !-----------------------------------------------------------------------------
75  !
76  !++ Private parameters & variables
77  !
78  ! KF subroutine look up tables VV
79  integer, private,PARAMETER :: kfnt=250,kfnp=220
80  real(RP), private,SAVE :: ttab(kfnt,kfnp),qstab(kfnt,kfnp)
81  real(RP), private,SAVE :: the0k(kfnp)
82  real(RP), private,SAVE :: alu(200)
83  real(RP), private,SAVE :: rdpr,rdthk,plutop
84  !^^^^
85  real(RP), private,SAVE :: gdcp ! GRAV/CP_dry
86  !
87  ! ALIQ saturrate watervape
88  ! BLIQ is WRF SVP2
89  ! DLIQ is WRF SVP3
90  real(RP),private,parameter :: aliq=6.112e2_rp ! Saturate pressure of water vapor [Pa]
91  real(RP),private,parameter :: bliq=17.68_rp ! emanuel 1994 (4.6.2) 17.67
92  real(RP),private,parameter :: cliq=bliq*tem00 ! convert degree to kelvin @ dew point temperature
93  real(RP),private,parameter :: dliq=29.65_rp ! 273.15 - 243.5
94  real(RP),private,parameter :: xlv1=2370._rp,xlv0=3.15e6_rp
95  ! Naming list tuning parameters
96  ! RATE is used subroutine precipitation_OC1973
97  real(RP),private, save :: rate = 0.03_rp ! ratio of cloud water and precipitation (Ogura and Cho 1973)
98  integer ,private, save :: trigger = 3 ! triger select will be modifid
99  logical ,private, save :: flag_qs = .true. ! FLAG_OS: qs is allowed or not
100  logical ,private, save :: flag_qi = .true. ! FLAG_QI: qi is allowe or not
101  real(RP),private, save :: delcape = 0.1_rp ! cape decleace rate
102  real(RP),private, save :: deeplifetime = 1800._rp ! minimum lifetimescale of deep convection
103  real(RP),private, save :: shallowlifetime = 2400._rp ! shallow convection liftime
104  real(RP),private, save :: depth_usl = 300._rp ! depth of updraft source layer [hPa]!!
105  logical ,private, save :: warmrain = .false. ! QA<=3 then ??
106  logical ,private, save :: kf_log = .false. ! KF infomation output to log file(not ERROR messeage)
107  real(RP),private, save :: kf_threshold = 1.e-3_rp ! kessler type autoconversion rate
108  integer ,private, save :: stepkf !! triger select will be modifid
109  integer ,private, save :: kf_prec = 1 ! precipitation select 1. Ogura and Cho (1973), 2. Kessler
110  procedure(precipitation), pointer,private :: p_precipitation => null()
111 
112  real(RP), private, allocatable :: lifetime (:,:) ! convectime lifetime [s]
113  integer , private, allocatable :: i_convflag(:,:) ! convection type 0:deep convection 1:shallow convection 2: no convection
114 
115  real(RP), private, allocatable :: deltaz (:,:,:) ! height interval (center level) [m]
116  real(RP), private, allocatable :: z (:,:,:) ! centerlevel real height [m]
117  real(RP) :: deltax ! delta x [m]
118 
119 
120  ! kf time controll
121  integer, private :: time_res_kf ! time step for kf
122  integer, private :: time_dstep_kf ! time interval
123  logical, private :: time_dokf ! exclude kf trigger
124 
125  ! tuning parameter
126  logical, private :: param_atmos_phy_cp_kf_wadapt = .true.
127  integer, private :: param_atmos_phy_cp_kf_w_time = 16
128 
129  !------------------------------------------------------------------------------
130 contains
131  !------------------------------------------------------------------------------
133  subroutine atmos_phy_cp_kf_setup (CP_TYPE)
134  use scale_process, only: &
135  prc_mpistop
136  use scale_time , only :&
137  time_dtsec, &
138  time_dtsec_atmos_phy_cp
139  use scale_grid_real, only: &
140  cz => real_cz
141  use scale_grid,only: &
142  dx => dx, &
143  dy => dy
144  implicit none
145 
146  character(len=*), intent(in) :: CP_TYPE
147 
148  ! tunning parameters, original parameter set is from KF2004 and NO2007
149  real(RP) :: PARAM_ATMOS_PHY_CP_kf_rate = 0.03_rp ! ratio of cloud water and precipitation (Ogura and Cho 1973)
150  integer :: PARAM_ATMOS_PHY_CP_kf_trigger = 1 ! trigger function type 1:KF2004 3:NO2007
151  logical :: PARAM_ATMOS_PHY_CP_kf_qs = .true. ! qs is allowed?
152  logical :: PARAM_ATMOS_PHY_CP_kf_qi = .true. ! qi is allowed?
153  real(DP) :: PARAM_ATMOS_PHY_CP_kf_dt = 5.0_dp ! KF convection check time interval [min]
154  real(RP) :: PARAM_ATMOS_PHY_CP_kf_dlcape = 0.1_rp ! cape decleace rate
155  real(RP) :: PARAM_ATMOS_PHY_CP_kf_dlifetime = 1800.0_rp ! minimum lifetime scale of deep convection [sec]
156  real(RP) :: PARAM_ATMOS_PHY_CP_kf_slifetime = 2400.0_rp ! lifetime of shallow convection [sec]
157  real(RP) :: PARAM_ATMOS_PHY_CP_kf_DEPTH_USL = 300.0_rp ! depth of updraft source layer [hPa]
158  integer :: PARAM_ATMOS_PHY_CP_kf_prec = 1 ! precipitation type 1:Ogura-Cho(1973) 2:Kessler
159  real(RP) :: PARAM_ATMOS_PHY_CP_kf_thres = 1.e-3_rp ! autoconversion rate for Kessler
160  logical :: PARAM_ATMOS_PHY_CP_kf_warmrain = .false. ! QQA is less equal to 3?
161  logical :: PARAM_ATMOS_PHY_CP_kf_LOG = .false. ! output log?
162 
163  namelist / param_atmos_phy_cp_kf / &
164  param_atmos_phy_cp_kf_rate, &
165  param_atmos_phy_cp_kf_trigger, &
166  param_atmos_phy_cp_kf_qs, &
167  param_atmos_phy_cp_kf_qi, &
168  param_atmos_phy_cp_kf_dt, &
169  param_atmos_phy_cp_kf_dlcape, &
170  param_atmos_phy_cp_kf_dlifetime, &
171  param_atmos_phy_cp_kf_slifetime, &
172  param_atmos_phy_cp_kf_depth_usl, &
173  param_atmos_phy_cp_kf_prec, &
174  param_atmos_phy_cp_kf_thres, &
175  param_atmos_phy_cp_kf_warmrain, &
176  param_atmos_phy_cp_kf_log, &
177  param_atmos_phy_cp_kf_wadapt, &
178  param_atmos_phy_cp_kf_w_time
179 
180  integer :: k, i, j
181  integer :: ierr
182  !---------------------------------------------------------------------------
183 
184  if( io_l ) write(io_fid_log,*)
185  if( io_l ) write(io_fid_log,*) '++++++ Module[CUMULUS] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
186  if( io_l ) write(io_fid_log,*) '+++ Kain-Fritsch scheme'
187 
188  if ( cp_type /= 'KF' ) then
189  write(*,*) 'xxx ATMOS_PHY_CP_TYPE is not KF. Check!'
190  call prc_mpistop
191  endif
192 
193  if ( abs(time_dtsec_atmos_phy_cp-time_dtsec) > 0.0_dp ) then
194  write(*,*) 'xxx TIME_DTSEC_ATMOS_PHY_CP should be same as TIME_DTSEC for KF scheme. STOP.'
195  call prc_mpistop
196  endif
197 
198  !--- read namelist
199  rewind(io_fid_conf)
200  read(io_fid_conf,nml=param_atmos_phy_cp_kf,iostat=ierr)
201  if( ierr < 0 ) then !--- missing
202  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
203  elseif( ierr > 0 ) then !--- fatal error
204  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_CP_KF. Check!'
205  call prc_mpistop
206  endif
207  if( io_l ) write(io_fid_log,nml=param_atmos_phy_cp_kf)
208 
209  call kf_lutab ! set kf look up table
210 
211  ! set kf convection check step
212  time_dstep_kf = nint( param_atmos_phy_cp_kf_dt * 60.0_dp / time_dtsec_atmos_phy_cp )
213  time_dstep_kf = max(time_dstep_kf,1) ! kf time interval step
214  time_res_kf = -1 ! initialize to keep consistent for below step
215  time_dokf = .true. ! initialize
216 
217  call cp_kf_param( param_atmos_phy_cp_kf_rate, &
218  param_atmos_phy_cp_kf_trigger, &
219  param_atmos_phy_cp_kf_qs, &
220  param_atmos_phy_cp_kf_qi, &
221  param_atmos_phy_cp_kf_dt, &
222  param_atmos_phy_cp_kf_dlcape, &
223  param_atmos_phy_cp_kf_dlifetime, &
224  param_atmos_phy_cp_kf_slifetime, &
225  param_atmos_phy_cp_kf_depth_usl, &
226  param_atmos_phy_cp_kf_prec, &
227  param_atmos_phy_cp_kf_thres, &
228  param_atmos_phy_cp_kf_warmrain, &
229  param_atmos_phy_cp_kf_log , &
230  time_dstep_kf )
231 
232  ! output parameter lists
233  if( io_l ) write(io_fid_log,*)
234  if( io_l ) write(io_fid_log,*) "*** Interval for check [step] : ", time_dstep_kf
235  if( io_l ) write(io_fid_log,*) "*** Ogura-Cho condense material convert rate : ", param_atmos_phy_cp_kf_rate
236  if( io_l ) write(io_fid_log,*) "*** Trigger function type, 1:KF2004 3:NO2007 : ", param_atmos_phy_cp_kf_trigger
237  if( io_l ) write(io_fid_log,*) "*** Exist qi? : ", param_atmos_phy_cp_kf_qi
238  if( io_l ) write(io_fid_log,*) "*** Exist qs? : ", param_atmos_phy_cp_kf_qs
239  if( io_l ) write(io_fid_log,*) "*** CAPE decrease rate : ", param_atmos_phy_cp_kf_dlcape
240  if( io_l ) write(io_fid_log,*) "*** Minimum lifetime scale of deep convection [sec] : ", param_atmos_phy_cp_kf_dlifetime
241  if( io_l ) write(io_fid_log,*) "*** Lifetime of shallow convection [sec] : ", param_atmos_phy_cp_kf_slifetime
242  if( io_l ) write(io_fid_log,*) "*** Updraft souce layer depth [hPa] : ", param_atmos_phy_cp_kf_depth_usl
243  if( io_l ) write(io_fid_log,*) "*** Precipitation type 1:Ogura-Cho(1973) 2:Kessler : ", param_atmos_phy_cp_kf_prec
244  if( io_l ) write(io_fid_log,*) "*** Kessler type precipitation's threshold : ", param_atmos_phy_cp_kf_thres
245  if( io_l ) write(io_fid_log,*) "*** Warm rain? : ", param_atmos_phy_cp_kf_warmrain
246  if( io_l ) write(io_fid_log,*) "*** Use running mean of w in adaptive timestep? : ", param_atmos_phy_cp_kf_wadapt
247  if( io_l ) write(io_fid_log,*) "*** Fixed time scale for running mean of w : ", param_atmos_phy_cp_kf_w_time
248  if( io_l ) write(io_fid_log,*) "*** Output log? : ", param_atmos_phy_cp_kf_log
249 
250  ! output variables
251  allocate( lifetime(ia,ja) )
252  allocate( i_convflag(ia,ja) )
253  lifetime(:,:) = 0.0_rp
254  i_convflag(:,:) = 2
255 
256  allocate( z(ka,ia,ja) )
257  z(:,:,:) = cz(:,:,:) ! becouse scale_atmos_phy_cp interface ,not use scale_grid
258 
259  allocate( deltaz(ka,ia,ja) )
260  ! deltaz is the interval of between model full levels(scalar point )
261  deltaz(:,:,:) = 0.0_rp
262  do j = js, je
263  do i = is, ie
264  do k = ks, ke
265  deltaz(k,i,j) = cz(k+1,i,j) - cz(k,i,j)
266  enddo
267  enddo
268  enddo
269  deltaz(ke,:,:) = 0.0_rp
270 
271  deltax = sqrt( dx*dy )
272 
273  return
274  end subroutine atmos_phy_cp_kf_setup
275 
276  !------------------------------------------------------------------------------
277  subroutine atmos_phy_cp_kf( &
278  DENS, &
279  MOMZ, &
280  MOMX, &
281  MOMY, &
282  RHOT, &
283  QTRC, &
284  DENS_t_CP, &
285  MOMZ_t_CP, &
286  MOMX_t_CP, &
287  MOMY_t_CP, &
288  RHOT_t_CP, &
289  RHOQ_t_CP, &
290  MFLX_cloudbase, &
291  SFLX_convrain, &
292  cloudtop, &
293  cloudbase, &
294  cldfrac_dp, &
295  cldfrac_sh, &
296  nca, &
297  w0avg )
298  use scale_grid_index
299  use scale_tracer, only: &
300  mp_qa
301  use scale_history, only: &
302  hist_in
303  implicit none
304 
305  real(RP), intent(in) :: DENS (ka,ia,ja)
306  real(RP), intent(in) :: MOMX (ka,ia,ja)
307  real(RP), intent(in) :: MOMY (ka,ia,ja)
308  real(RP), intent(in) :: MOMZ (ka,ia,ja)
309  real(RP), intent(in) :: RHOT (ka,ia,ja)
310  real(RP), intent(in) :: QTRC (ka,ia,ja,qa)
311  real(RP), intent(inout) :: DENS_t_CP (ka,ia,ja)
312  real(RP), intent(inout) :: MOMZ_t_CP (ka,ia,ja) ! not used
313  real(RP), intent(inout) :: MOMX_t_CP (ka,ia,ja) ! not used
314  real(RP), intent(inout) :: MOMY_t_CP (ka,ia,ja) ! not used
315  real(RP), intent(inout) :: RHOT_t_CP (ka,ia,ja)
316  real(RP), intent(inout) :: RHOQ_t_CP (ka,ia,ja,qa)
317  real(RP), intent(inout) :: MFLX_cloudbase(ia,ja) ! not used
318  real(RP), intent(inout) :: SFLX_convrain (ia,ja) ! convective rain rate [kg/m2/s]
319  real(RP), intent(inout) :: cloudtop (ia,ja) ! cloud top height [m]
320  real(RP), intent(inout) :: cloudbase (ia,ja) ! cloud base height [m]
321  real(RP), intent(inout) :: cldfrac_dp (ka,ia,ja) ! cloud fraction (deep convection)
322  real(RP), intent(inout) :: cldfrac_sh (ka,ia,ja) ! cloud fraction (shallow convection)
323  real(RP), intent(inout) :: nca (ia,ja) ! convection active time [sec]
324  real(RP), intent(inout) :: w0avg (ka,ia,ja) ! running mean of vertical velocity [m/s]
325 
326  real(RP) :: cldfrac(ka,2) ! 1 shallow , 2 deep
327 
328  !---------------------------------------------------------------------------
329 
330  if( io_l ) write(io_fid_log,*) '*** Physics step: Cumulus Parameterization(KF)'
331 
332  call kf_wmean( w0avg(:,:,:), & ! [OUT]
333  dens(:,:,:), & ! [IN]
334  momz(:,:,:) ) ! [IN]
335 
336  time_dokf = .false.
337  time_res_kf = time_res_kf + 1
338  if ( time_res_kf == time_dstep_kf ) then
339  time_dokf = .true.
340  time_res_kf = 0
341  endif
342 
343  if ( time_dokf ) then
344  if( io_l ) write(io_fid_log,*) '*** KF Convection Check '
345 
346  call prof_rapstart('CP_kf', 3)
347 
348  ! calc cumulus convection
349  call cp_kf_main( dens(:,:,:), & ! [IN]
350  momz(:,:,:), & ! [IN]
351  momx(:,:,:), & ! [IN]
352  momy(:,:,:), & ! [IN]
353  rhot(:,:,:), & ! [IN]
354  qtrc(:,:,:,:), & ! [IN]
355  w0avg(:,:,:), & ! [IN]
356  nca(:,:), & ! [INOUT]
357  dens_t_cp(:,:,:), & ! [OUT]
358  rhot_t_cp(:,:,:), & ! [OUT]
359  rhoq_t_cp(:,:,:,:), & ! [OUT]
360  sflx_convrain(:,:), & ! [OUT]
361  cldfrac_sh(:,:,:), & ! [OUT]
362  cldfrac_dp(:,:,:), & ! [OUT]
363  lifetime(:,:), & ! [OUT]
364  cloudtop(:,:), & ! [OUT]
365  cloudbase(:,:), & ! [OUT]
366  i_convflag(:,:) ) ! [OUT]
367 
368  call prof_rapend('CP_kf', 3)
369  endif
370 
371  call hist_in( lifetime(:,:), 'KF_LIFETIME', 'lifetime of KF scheme', 's' )
372  call hist_in( real(I_convflag(:,:),RP), 'KF_CONVFLAG', 'CONVECTION FLAG', '' )
373 
374  return
375  end subroutine atmos_phy_cp_kf
376 
377  !-----------------------------------------------------------------------------
379  ! WRF comment out for W0
380  !...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
381  !...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
382  !...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
383  !...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
384  !...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
385  subroutine kf_wmean( &
386  W0_avg, &
387  DENS, &
388  MOMZ )
389  use scale_time , only :&
390  time_dtsec
391  implicit none
392 
393  real(RP), intent(inout) :: W0_avg(ka,ia,ja)
394  real(RP), intent(in) :: DENS (ka,ia,ja)
395  real(RP), intent(in) :: MOMZ (ka,ia,ja)
396 
397  real(RP) :: W0
398  real(RP) :: fact_W0_avg, fact_W0
399 
400  integer :: k, i, j
401  !---------------------------------------------------------------------------
402 
403  if ( param_atmos_phy_cp_kf_wadapt ) then
404  fact_w0_avg = 2.0_rp * max(kf_dt,time_dtsec) - time_dtsec
405  fact_w0 = time_dtsec
406  else ! w_time is tuning parameter
407  fact_w0_avg = real(param_atmos_phy_cp_kf_w_time,rp)
408  fact_w0 = 1.0_rp
409  endif
410 
411  do j = js, je
412  do i = is, ie
413  do k = ks, ke
414  w0 = 0.5_rp * ( momz(k,i,j) + momz(k-1,i,j) ) / dens(k,i,j)
415 
416  w0_avg(k,i,j) = ( w0_avg(k,i,j) * fact_w0_avg &
417  + w0 * fact_w0 ) / ( fact_w0_avg + fact_w0 )
418  enddo
419  enddo
420  enddo
421 
422  return
423  end subroutine kf_wmean
424 
425  !-----------------------------------------------------------------------------
426  subroutine cp_kf_param( & ! set kf tuning parametres
427  ![IN]
428  RATE_in, &
429  TRIGGER_in, & ! INOUT
430  FLAG_QS_in, &
431  FLAG_QI_in, &
432  KF_DT_in, &
433  DELCAPE_in , &
434  DEEPLIFETIME_in, &
435  SHALLOWLIFETIME_in, &
436  DEPTH_USL_in, &
437  KF_prec_in, &
438  KF_threshold_in, &
439  WARMRAIN_in, &
440  KF_LOG_in, &
441  stepkf_in )
442  use scale_process, only: &
443  prc_mpistop
444  implicit none
445  real(RP),intent(in) :: RATE_in
446  integer, intent(inout) :: TRIGGER_in
447  logical, intent(in) :: FLAG_QS_in,FLAG_QI_in
448  real(RP),intent(in) :: KF_DT_in
449  real(RP),intent(in) :: DELCAPE_in
450  real(RP),intent(in) :: DEEPLIFETIME_in
451  real(RP),intent(in) :: SHALLOWLIFETIME_in
452  real(RP),intent(in) :: DEPTH_USL_in
453  integer, intent(in) :: KF_prec_in
454  real(RP),intent(in) :: KF_threshold_in
455  logical, intent(in) :: WARMRAIN_in
456  logical, intent(in) :: KF_LOG_in
457  integer, intent(in) :: stepkf_in
458  !
459  rate = rate_in
460  ! TRIGGER must be 1 or 3
461  if (trigger_in /= 1 .and. trigger_in /=3) then
462  if (io_l) write(io_fid_log,*) "TRIGGER must be 1 or 3 but now :",trigger_in
463  if (io_l) write(io_fid_log,*) "CHAGNGE ",trigger_in," to 3"
464  trigger_in = 3
465  end if
466  trigger = trigger_in
467  flag_qs = flag_qs_in
468  flag_qi = flag_qi_in
469  kf_dt = kf_dt_in
470  delcape = delcape_in
471  deeplifetime = deeplifetime_in
472  shallowlifetime = shallowlifetime_in
473  depth_usl = depth_usl_in
474  warmrain = warmrain_in
475  kf_prec = kf_prec_in
476  kf_threshold = kf_threshold_in
477  kf_log = kf_log_in
478  stepkf = stepkf_in
479  if (kf_prec == 1) then
480  p_precipitation => precipitation_oc1973 ! Ogura and Cho (1973)
481  elseif( kf_prec == 2) then
482  p_precipitation => precipitation_kessler ! Kessler type
483  else
484  write(*,*) 'xxx ERROR at KF namelist'
485  write(*,*) 'KF_prec must be 1 or 2'
486  call prc_mpistop
487  end if
488  return
489  end subroutine cp_kf_param
490 
491  subroutine cp_kf_main ( & ! main loop
492  ![IN]
493  dens, &
494  MOMZ, &
495  MOMX, &
496  MOMY, &
497  RHOT, &
498  QTRC, &
499  w0avg, &
500  ! [INOUT]
501  nca, &
502  ! [OUT]
503  DENS_t_CP, &
504  DT_RHOT, &
505  DT_RHOQ, &
506  rainrate_cp, &
507  cldfrac_sh, &
508  cldfrac_dp, &
509  timecp, &
510  cloudtop, &
511  zlcl, &
512  I_convflag )
513  use scale_precision
514  use scale_grid_index
515  use scale_tracer
516  use scale_const, only: &
517  grav => const_grav, &
518  r => const_rdry
519  use scale_time , only :&
520  dt => time_dtsec_atmos_phy_cp
521  use scale_grid_real, only: &
522  fz => real_fz
523  use scale_atmos_thermodyn, only: &
524  thermodyn_temp_pres => atmos_thermodyn_temp_pres, &
525  thermodyn_rhoe => atmos_thermodyn_rhoe, &
526  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e, &
527  thermodyn_qd => atmos_thermodyn_qd, &
528  thermodyn_pott => atmos_thermodyn_pott
529  use scale_atmos_saturation ,only :&
530  saturation_psat_liq => atmos_saturation_psat_liq
531  implicit none
532  ![IN]
533  real(RP),intent(in) :: DENS(ka,ia,ja) ! density [kg/m**3]
534  real(RP),intent(in) :: MOMZ(ka,ia,ja) ! momentum
535  real(RP),intent(in) :: MOMX(ka,ia,ja) ! momentum
536  real(RP),intent(in) :: MOMY(ka,ia,ja) ! momentum
537  real(RP),intent(in) :: RHOT(ka,ia,ja) ! density*PT
538  real(RP),intent(in) :: QTRC(ka,ia,ja,qa) ! raito of water elements
539  real(RP),intent(in) :: w0avg(ka,ia,ja) ! running mean vertical wind velocity [m/s]
540  ! [INOUT]
541  real(RP),intent(inout) :: nca(ia,ja) ! num of step convection active [step]
542  real(RP),intent(out) :: DENS_t_CP(ka,ia,ja) ! dens/dt
543  real(RP),intent(out) :: DT_RHOT(ka,ia,ja) ! dens*PT/dt
544  real(RP),intent(out) :: DT_RHOQ(ka,ia,ja,qa) ! dens*q/dt
545  real(RP),intent(out) :: rainrate_cp(ia,ja) ! rain PPTFLX(prcp_flux)/deltax**2*dt ! convective rain
546  real(RP),intent(out) :: cldfrac_sh(ka,ia,ja) ! cloud fraction
547  real(RP),intent(out) :: cldfrac_dp(ka,ia,ja) ! cloud fraction
548  real(RP),intent(out) :: timecp(ia,ja) ! timescale of cumulus parameterization
549  real(RP),intent(out) :: cloudtop(ia,ja) ! cloud top height
550  real(RP),intent(out) :: zlcl(ia,ja) ! hight of lcl cloud bottom height[m]
551  integer, intent(out) :: I_convflag(ia,ja) ! convection type
555 
556  ! [Internal Work]
557  real(RP) :: u (ka) ! x-direction wind velocity [m/s]
558  real(RP) :: v (ka) ! y-direction wind velocity [m/s]
559  real(RP) :: temp (ka) ! temperature [K]
560  real(RP) :: pres (ka) ! pressure [Pa]
561  real(RP) :: qv (ka) ! water vapor mixing ratio[kg/kg]
562  real(RP) :: QDRY (ka) ! ratio of dry air
563  real(RP) :: qes (ka) ! saturate vapor [kg/kg]
564  real(RP) :: PSAT (ka) ! saturation vaper pressure
565  real(RP) :: QSAT (ka) ! saturate water vaper mixing ratio [kg/kg]
566  real(RP) :: rh (ka) ! saturate vapor [%]
567  real(RP) :: deltap(ka) ! delta Pressure [Pa]
568 
569  real(RP) :: q_hyd(ka,mp_qa) ! water mixing ratio [kg/kg]
570  real(RP) :: dens_nw(ka) ! density [kg/m**3]
571  integer :: nic
572  real(RP) :: time_advec ! advection timescale
573  real(RP) :: umf(ka) ! Updraft Mass Flux
574  real(RP) :: umflcl ! Updraft Mass Flux
575  real(RP) :: upent(ka) ! Updraft Mass flux ENTrainment
576  real(RP) :: updet(ka) ! Updraft Mass flux DETrainment
577  ! updraft value
578  real(RP) :: temp_u(ka) ! updraft temperature [K]
579  real(RP) :: tempv(ka) ! vertual temperature [K]
580  real(RP) :: qv_u(ka) ! updraft qv
581  real(RP) :: cape ! cape
582  ! water
583  real(RP) :: qc(ka) ! cloud water
584  real(RP) :: qi(ka) ! cloud ice
585  real(RP) :: qvdet(ka) ! vapor detrainment
586  real(RP) :: qcdet(ka) ! liquit detrainment
587  real(RP) :: qidet(ka) ! ice detrainment
588  real(RP) :: totalprcp ! total precipitation
589  real(RP) :: flux_qr(ka) ! rain flux
590  real(RP) :: flux_qs(ka) ! snow flux
591  ! upward theta
592  real(RP) :: theta_eu(ka) ! updraft equivalent PT
593  real(RP) :: theta_ee(ka) ! environment equivalent PT
594  ! convection type
595  ! index valuables
596  integer :: k_lcl ! index of LCL layer
597  integer :: k_top ! index of cloud top hight
598  integer :: k_ml ! index of melt layer (temp < tem00)
599  integer :: k_lc,k_let,k_pbl ! indexs
600  real(RP) :: zmix ! usl(up source layer) layer depth [m]
601  real(RP) :: presmix ! usl layer depth [Pa]
602  real(RP) :: umfnewdold(ka) ! umfnew/umfold
603  real(RP) :: dpthmx ! max depth of pressure
604  ! move internalwark
605  real(RP) :: ems(ka) ! box weight[kg]
606  real(RP) :: emsd(ka) ! 1/ems
607  ! [OUT]@downdraft
608  real(RP) :: wspd(3) ! horizontal wind speed 1 k_lcl, 2 k_z5,3 k_top
609  real(RP) :: dmf(ka) ! downdraft massflux
610  real(RP) :: downent(ka) ! downdraft entrainment
611  real(RP) :: downdet(ka) ! downdraft detrainment
612  real(RP) :: theta_d(ka) ! downdraft theta
613  real(RP) :: qv_d(ka) ! water vapor of downdraft
614  real(RP) :: prcp_flux ! precpitation flux
615  real(RP) :: tder ! temperature change from evap
616  real(RP) :: CPR ! all precipitation before consider cloud bottom evaporation
617  integer :: k_lfs ! LFS layer index
618  ![OUT] @ compensasional subsidence
619  real(RP) :: temp_g(ka) ! temporarly temperature -> after iteration then new variable
620  real(RP) :: qv_g(ka) ! temporarly vapor -> after iteration then new variable
621  real(RP) :: qc_nw(ka) ! new cloud water mixing ratio [kg/kg]
622  real(RP) :: qi_nw(ka) ! new cloud ice mixing ratio [kg/kg]
623  real(RP) :: qr_nw(ka) ! new rain water mixing ratio [kg/kg]
624  real(RP) :: qs_nw(ka) ! new snow water mixing ratio [kg/kg]
625  ! new variable
626  real(RP) :: rhot_nw(ka) ! rho*PT
627  real(RP) :: qtrc_nw(ka,qa) ! qv,qc,qr,qi,qs (qg not change)
628  real(RP) :: pott_nw(ka) ! new PT
629  !
630  real(RP) :: cldfrac_KF(ka,2) ! cloud fraction
631  !
632  real(RP) :: qd(ka) ! dry mixing ratio
633  ! do loop variable
634  integer :: k, i, j, iq
635 
636 
637  do j = js, je
638  do i = is, ie
639 
640  nca(i,j) = nca(i,j) - real(TIME_DSTEP_KF,RP) * dt
641 
642  ! check convection
643  if ( nca(i,j) .ge. 0.5_dp * dt ) cycle
644 
645 
646  ! convert variables
647 
648  ! calculate u(x-directin velocity ), v(y-direction velocity)
649  do k = ks, ke
650  u(k) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
651  v(k) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
652  enddo
653 
654  do k = ks, ke
655  call thermodyn_temp_pres( temp(k), & ! [OUT]
656  pres(k), & ! [OUT]
657  dens(k,i,j), & ! [IN]
658  rhot(k,i,j), & ! [IN]
659  qtrc(k,i,j,:) ) ! [IN]
660 
661  ! calculate water vaper and relative humidity
662  call thermodyn_qd( qdry(k), qtrc(k,i,j,:) )
663 
664  call saturation_psat_liq( psat(k), temp(k) )
665 
666  qsat(k) = 0.622_rp * psat(k) / ( pres(k) - ( 1.0_rp-0.622_rp ) * psat(k) )
667  qv(k) = qtrc(k,i,j,i_qv) / qdry(k)
668  qv(k) = max( 0.000001_rp, min( qsat(k), qv(k) ) ) ! conpare QSAT and QV, guess lower limit
669  rh(k) = qv(k) / qsat(k)
670  enddo
671 
672  ! calculate delta P by hydrostatic balance
673  ! deltap is the pressure interval between half levels(face levels) @ SCALE
674  do k = ks, ke
675  deltap(k) = dens(k,i,j) * grav * ( fz(k+1,i,j) - fz(k,i,j) ) ! rho*g*dz
676  enddo
677 
678 
679  dens_t_cp(ks:ke,i,j) = 0.0_rp
680  dt_rhot(ks:ke,i,j) = 0.0_rp
681  do iq = 1, qqs
682  dt_rhoq(ks:ke,i,j,iq) = 0.0_rp
683  end do
684  cldfrac_kf(ks:ke,:) = 0.0_rp
685  rainrate_cp(i,j) = 0.0_rp
686  timecp(i,j) = 0.0_rp
687  cloudtop(i,j) = 0.0_rp
688  zlcl(i,j) = 0.0_rp
689 
690  do iq = 1, mp_qa
691  do k = ks, ke
692  q_hyd(k,iq) = qtrc(k,i,j,i_mp2all(iq)) / qdry(k)
693  enddo
694  enddo
695 
696 
697  call cp_kf_trigger ( &
698  ! [IN]
699  deltaz(:,i,j), z(:,i,j), & ! deltaz and height[m]
700  qv(:), & ! water vapor mixing ratio [kg/kg]
701  qes(:), & ! saturation water vapor mixing ratio
702  pres(:), & ! pressure [Pa]
703  deltap(:), & ! interval of pressure [Pa]
704  temp(:), & ! temperature [K]
705  w0avg(:,i,j), & ! average w
706  ! [OUT]
707  i_convflag(i,j), & ! convection flag
708  cloudtop(i,j), & ! cloud top height [m]
709  temp_u(:), & ! updraft temperature
710  tempv(:), & ! virtual temperature
711  ! water values
712  qv_u(:), & ! updraft water vapor mixing ratio
713  qc(:), & ! cloud water mixing ratio
714  qi(:), & ! ice mixing ratio
715  qvdet(:), & ! updraft detrain vapor
716  qcdet(:), & ! updraft detrain water
717  qidet(:), & ! updraft detrain ice
718  flux_qr(:), & ! rain flux
719  flux_qs(:), & ! snow flux
720  totalprcp, & ! total precipitation(rain and snow ) fall
721  ! thermodynamics values
722  theta_eu(:), & ! updraft equivalent PT
723  theta_ee(:), & ! environment equivalent PT
724  cape, & ! cape
725  ! massflux values
726  umf(:), & ! updraft mass flux
727  umflcl, & ! updraft mass flux@lcl
728  upent(:), & ! updraft entrainment
729  updet(:), & ! updraft detrainment
730  ! layer nambers
731  k_lcl, & ! index of LCL layer
732  k_lc, & ! index of USL layer botom
733  k_pbl, & ! index of USL layer top
734  k_top, & ! index of cloud top
735  k_let, & ! index of below detrain layer
736  k_ml, & ! index of temprerure < tem00
737  ! pbl layer values
738  presmix, & ! USL layer pressure [Pa]
739  dpthmx, & ! USL layer depth ( [Pa])
740  zlcl(i,j), & ! lcl hight
741  zmix, & ! USL layer depth
742  umfnewdold(:) )
743  !------------------------------------------------------------------------
744  ! calc ems(box weight[kg])
745  !------------------------------------------------------------------------
746  if(i_convflag(i,j) /= 2) then
747  ems(k_top+1:ke) = 0._rp
748  emsd(k_top+1:ke) = 0._rp
749  do k = ks, k_top
750  ems(k) = deltap(k)*deltax*deltax/grav
751  emsd(k) = 1._rp/ems(k)
752  umfnewdold(k) = 1._rp/umfnewdold(k)
753  end do
754  end if
755  !------------------------------------------------------------------------
756  call cp_kf_downdraft ( &
757  ! [IN]
758  deltaz(:,i,j), z(:,i,j) ,& ! deltaz and height [m]
759  u(:) ,& ! u-wind m/s
760  v(:) ,& ! v-wind m/s
761  zlcl(i,j) ,& ! lcl height [m]
762  rh(:) ,& ! relative humidity
763  deltap(:) ,& ! interval of pressure [Pa]
764  pres(:) ,& ! pressure [Pa]
765  qv(:) ,& ! water vapor mixing ratio [kg/kg]
766  ems(:) ,& ! ems(box weight[kg])
767  emsd(:) ,& ! 1/ems
768  theta_ee(:) ,& ! environment theta_E
769  umf(:) ,& ! updraft mass flux
770  totalprcp ,& ! total precipitation
771  flux_qs(:) ,& ! snow flux
772  tempv(:) ,& ! vertual temperature [K]
773  i_convflag(i,j) ,& ! convection type
774  ! layer index
775  k_lcl ,& ! index of LCL layer
776  k_ml ,& ! index of melt layer
777  k_top ,& ! index of cloud top
778  k_pbl ,& ! index of USL layer top
779  k_let ,& ! index of below detrain layer
780  k_lc ,& ! index of USL layer botom
781  ! [out]
782  wspd(:) ,& ! wind speed 1 k_lcl, 2 k_z5,3 k_top
783  dmf(:) ,& ! downward mass flux
784  downent(:) ,& ! downdraft entrainment
785  downdet(:) ,& ! downdraft detrainment
786  theta_d(:) ,& ! downdraft theta
787  qv_d(:) ,& ! downdraft water vaper mixng ratio
788  prcp_flux ,& ! precipitateion
789  k_lfs ,& ! index of LFS layer
790  cpr ,& ! all precipitation before consider cloud bottom evaporation
791  tder)
792  !------------------------------------------------------------------------
793  call cp_kf_compensational ( &
794  ![IN]
795  deltaz(:,i,j),z(:,i,j) ,& ! deltaz and height [m]
796  pres(:), & ! pressure [Pa]
797  deltap(:), & ! deltap [Pa]
798  temp(:), & ! temperature [K]
799  qv(:), & ! water vapor mixing ratio
800  ! form kf_trigger
801  ems(:), & ! box weight[kg]
802  emsd(:), & ! 1/ems
803  presmix, & ! usl layer pressure depth
804  zmix, & ! usl layer pressure depth
805  dpthmx, & ! usl layer depth
806  cape, & ! CAPE
807  temp_u(:), & ! updraft temperature
808  qvdet(:), & ! detrainment water vaper mixing ratio
809  umflcl, & ! umf @LCL
810  umf(:), & ! upward mass flux (updraft)
811  upent(:), & ! updraft entrainment
812  updet(:), & ! updraft detrainment
813  qcdet(:), & ! updraft detrainment qc
814  qidet(:), & ! updraft detrainment qi
815  qc(:), & ! cloud water mixing ratio
816  qi(:), & ! cloud ice mixing ratio
817  flux_qr(:), & ! rain flux
818  flux_qs(:), & ! snow flux
819  umfnewdold(:), & ! umfnew/umfold ratio
820  ! from kf_downdraft
821  wspd(:), & ! wind speed 1. 2 is used
822  dmf(:), & ! downward mass flux (downdraft)
823  downent(:), & ! downdraft entrainment
824  downdet(:), & ! downdraft detrainment
825  qv_d(:), & ! downdraft detrainment qv
826  theta_d(:), & ! potential temperature @downdraft
827  prcp_flux, & ! precipitation
828  tder, & ! evapolation
829  cpr, & ! all precipitation before consider cloud bottom evaporation
830  i_convflag(i,j), & ! intent inout convective flag
831  ! layer index
832  ! from kf_triger
833  k_top, & ! index of cloud top hight
834  k_lcl, & ! index of LCL layer
835  k_lc , & ! index of LCL layer botoom
836  k_pbl, & ! index of USL layer top
837  k_ml, & ! index of melt layer (temp < tem00)
838  ! from kf_downdraft
839  k_lfs, & ! index of LFS layer
840  ![OUT] new valuavles after timestep
841  temp_g(:), & ! update temperature [K]
842  qv_g(:), & ! update water vapor mixing ratio
843  qc_nw(:), & ! update cloud water mixing ratio
844  qi_nw(:), & ! update cloud ice mixing ratio
845  qr_nw(:), & ! update rain water mixing ratio
846  qs_nw(:), & ! update snow water mixing ratio
847  rainrate_cp(i,j), & ! rain PPTFLX(prcp_flux)/deltax**2*dt ! convective rain
848  nic, & ! (timescale of cumulus parameterization)/dt integer
849  cldfrac_kf, & ! cloud fraction
850  timecp(i,j), & ! timescale of cumulus parameterization
851  time_advec) ! advection timescale
852 
853  !------------------------------------------------------------------------
854  ! compute tendencys
855  !------------------------------------------------------------------------
856 
857  if(i_convflag(i,j) == 2) then ! no convection
858  rainrate_cp(i,j) = 0.0_rp
859  cldfrac_kf(ks:ke,:) = 0.0_rp
860  timecp(i,j) = 0.0_rp
861  cloudtop(i,j) = 0.0_rp
862  zlcl(i,j) = 0.0_rp
863  else ! convection allowed I_convflag=0 or 1
864  ! chek
865  !
866  !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
867  !
868  !...IF THE ADVECTIVE TIME PERIOD (time_advec) IS LESS THAN SPECIFIED MINIMUM
869  !...TIMECP, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
870  !
871  if (i_convflag(i,j) == 0) then ! deep
872  if (time_advec < timecp(i,j)) nic=nint(time_advec/dt)
873  nca(i,j) = real(nic,rp)*dt ! convection feed back act this time span
874  elseif (i_convflag(i,j) == 1) then ! shallow
875  timecp(i,j) = shallowlifetime
876  nca(i,j) = kf_dt*60._rp ! convection feed back act this time span
877  end if
878  ! update qd
879  if (qa > 3) then ! assume TOMITA08
880  q_hyd(ks:ke,i_qc-1) = ( qc_nw(ks:ke) + q_hyd(ks:ke,i_qc-1) )
881  q_hyd(ks:ke,i_qr-1) = ( qr_nw(ks:ke) + q_hyd(ks:ke,i_qr-1) )
882  q_hyd(ks:ke,i_qi-1) = ( qi_nw(ks:ke) + q_hyd(ks:ke,i_qi-1) )
883  q_hyd(ks:ke,i_qs-1) = ( qs_nw(ks:ke) + q_hyd(ks:ke,i_qs-1) )
884  qd(ks:ke) = 1._rp / ( 1._rp + qv_g(ks:ke) + sum(q_hyd(ks:ke,:),2)) ! new qdry
885  qtrc_nw(ks:ke,i_qv) = qv_g(ks:ke)*qd(ks:ke)
886  ! new qtrc
887  do iq = 1,mp_qa
888  qtrc_nw(ks:ke,i_mp2all(iq)) = ( q_hyd(ks:ke,iq) )*qd(ks:ke)
889  end do
890  ! new density
891  dens_nw(ks:ke) = dens(ks:ke,i,j)
892  do iq = qqs,qqe
893  dens_nw(ks:ke) = dens_nw(ks:ke) + ( qtrc_nw(ks:ke,iq) - qtrc(ks:ke,i,j,iq) )*dens(ks:ke,i,j)
894  end do
895  ! dens_(n+1) = dens_(n) + tend_q*dens_(n)
896  else ! kessular
897  q_hyd(ks:ke,i_qc-1) = ( qc_nw(ks:ke) + q_hyd(ks:ke,i_qc-1) )
898  q_hyd(ks:ke,i_qr-1) = ( qr_nw(ks:ke) + q_hyd(ks:ke,i_qr-1) )
899  !
900  qd(ks:ke) = 1._rp / ( 1._rp + qv_g(ks:ke) + sum(q_hyd(ks:ke,:),2))
901  ! new qtrc
902  do iq = 1,mp_qa
903  qtrc_nw(ks:ke,i_mp2all(iq)) = ( q_hyd(ks:ke,iq) )*qd(ks:ke)
904  end do
905  ! new density
906  dens_nw(ks:ke) = dens(ks:ke,i,j)
907  do iq = qqs,qqe
908  dens_nw(ks:ke) = dens_nw(ks:ke) + ( qtrc_nw(ks:ke,iq) - qtrc(ks:ke,i,j,iq) )*dens(ks:ke,i,j)
909  end do
910  ! dens_(n+1) = dens_(n) + tend_q*dens_(n)
911  end if
912  ! calc new potential temperature
913  call thermodyn_pott(&
914  pott_nw(:), &
915  temp_g(:),pres(:),qtrc_nw(:,:))
916  ! update rhot
917  rhot_nw(ks:ke) = dens_nw(ks:ke)*pott_nw(ks:ke)
918  ! calc tendency
919  dens_t_cp(ks:ke,i,j) = (dens_nw(ks:ke) - dens(ks:ke,i,j))/timecp(i,j)
920  dt_rhot(ks:ke,i,j) = (rhot_nw(ks:ke) - rhot(ks:ke,i,j))/timecp(i,j)
921  do iq = qqs,qqe
922  dt_rhoq(ks:ke,i,j,iq) = (dens_nw(ks:ke)*qtrc_nw(ks:ke,iq) - dens(ks:ke,i,j)*qtrc(ks:ke,i,j,iq))/timecp(i,j)
923  end do
924  ! if noconvection then nca is same value before call. nca only modifyed convectioned
925  end if
926 
927  cldfrac_sh(ks:ke,i,j) = cldfrac_kf(ks:ke,1)
928  cldfrac_dp(ks:ke,i,j) = cldfrac_kf(ks:ke,2)
929 
930  end do
931  end do
932 
933  return
934  end subroutine cp_kf_main
935  !------------------------------------------------------------------------------
936  ! calc trigger function and upward mass flux
937  !------------------------------------------------------------------------------
938  subroutine cp_kf_trigger ( & !> 1d model ! triger function
939  ! [IN]
940  dz_kf,z_kf ,& ! deltaz and height [m]
941  qv ,& ! water vapor mixing ratio [kg/kg]
942  qes ,& ! saturated vapor [kg/kg]
943  pres ,& ! pressure [Pa]
944  deltap ,& ! interval of pressure [Pa]
945  temp ,& ! temperature [K]
946  w0avg ,& ! running mean w
947  ! [OUT]
948  I_convflag ,& ! convection flag
949  cloudtop ,& ! cloud top height [m]
950  temp_u ,& ! updraft temperature
951  tempv ,& ! virtual temperature
952  ! water variable
953  qv_u ,& ! updraft water vapor mixing ratio
954  qc ,& ! cloud water mixing ratio
955  qi ,& ! ice mixingratio
956  qvdet ,& ! updraft detrain vapor
957  qcdet ,& ! updraft detrain water
958  qidet ,& ! updraft detrain ice
959  flux_qr ,& ! rain flux
960  flux_qs ,& ! snow flux
961  totalprcp ,& ! total precipitation (rain and snow fall)
962  ! thermodynamics variables
963  theta_eu ,& ! updraft equivalent theta
964  theta_ee ,& ! environment equivalent theta
965  cape ,& ! cape
966  ! massflux variables
967  umf ,& ! updraft mass flux
968  umflcl ,& ! updraft mass flux@LCL
969  upent ,& ! updraft entrainment
970  updet ,& ! updraft detrainment
971  ! layer indexs
972  k_lcl ,& ! index of LCL layer
973  k_lc ,& ! index of USL layer botom
974  k_pbl ,& ! index of USL layer top
975  k_top ,& ! index of cloud top
976  k_let ,& ! index of below detrain layer
977  k_ml ,& ! index of temprerure < tem00 (melting layer)
978  ! pbl layer variables
979  presmix ,& ! USL layer pressure
980  dpthmx ,& ! USL layer depth (pressure [Pa])
981  zlcl ,& ! LCL layer height[m]
982  zmix ,& ! usl layer depth [m]
983  umfnewdold ) ! ratio of umf/umfold(see updraft)
984  use scale_precision
985  use scale_grid_index
986  use scale_const,only :&
987  cp => const_cpdry , &
988  pre00 => const_pre00, &
989  tem00 => const_tem00, &
990  grav => const_grav
991  use scale_process, only: &
992  prc_mpistop
993  implicit none
994  ! [IN]
995  real(RP), intent(in) :: dz_kf(ka),z_kf(ka) ! delta z and height [m]
996  real(RP), intent(in) :: qv(ka) ! water vapor
997  real(RP), intent(in) :: qes(ka) ! saturation water vapor
998  real(RP), intent(in) :: pres(ka) ! pressure [Pa]
999  real(RP), intent(in) :: deltap(ka) ! delta pressure [Pa]
1000  real(RP), intent(in) :: temp(ka) ! temperature
1001  real(RP), intent(in) :: w0avg(ka) ! running mean w
1002  ! [OUT]
1003  ! mass flux
1004  real(RP), intent(out) :: umf(ka) ! upward mass flux
1005  real(RP), intent(out) :: umflcl ! upward mass flux @lcl
1006  real(RP), intent(out) :: upent(ka) ! upward mass flux entrainment
1007  real(RP), intent(out) :: updet(ka) ! upward mass flux detrainment
1008  ! updraft variable
1009  real(RP), intent(out) :: temp_u(ka) ! updraft temperature
1010  real(RP), intent(out) :: tempv(ka) ! vertual temperature
1011  ! real(RP) :: tempvq_u(KA)
1012  real(RP), intent(out) :: qv_u(ka) ! updraft qv
1013  real(RP), intent(out) :: totalprcp ! total prcpitation (rain+snow)
1014  real(RP), intent(out) :: cape ! CAPE
1015  real(RP), intent(out) :: cloudtop ! cloud top height
1016  real(RP) :: cloudhight ! cloud depth (cloud top - cloud base)
1017  ! water
1018  real(RP), intent(out) :: qc(ka) ! cloud water mixing ratio
1019  real(RP), intent(out) :: qi(ka) ! cloud ice mixing ratio
1020  real(RP) :: qrout(ka) ! rain
1021  real(RP) :: qsout(ka) ! snow
1022  real(RP), intent(out) :: qvdet(ka) ! detrainment water vapor
1023  real(RP), intent(out) :: qcdet(ka) ! detrainment cloud water
1024  real(RP), intent(out) :: qidet(ka) ! detrainment cloud ice
1025  real(RP), intent(out) :: flux_qr(ka) !
1026  real(RP), intent(out) :: flux_qs(ka) !
1027  ! upward theta
1028  real(RP), intent(out) :: theta_eu(ka)
1029  real(RP), intent(out) :: theta_ee(ka)
1030  ! convection type
1031  integer, intent(out) :: I_convflag
1035  ! index valuables
1036  integer, intent(out) :: k_lcl ! index of LCL layer
1037  integer, intent(out) :: k_top ! index of cloud top hight
1038  integer, intent(out) :: k_ml ! index of melt layer (temp < tem00)
1039  real(RP), intent(out) :: zlcl ! hight of lcl
1040  integer, intent(out) :: k_lc,k_let,k_pbl ! indexs
1041  real(RP), intent(out) :: zmix ! usl layer depth [m]
1042  real(RP), intent(out) :: presmix ! usl layer depth [Pa]
1043  real(RP), intent(out) :: umfnewdold(ka) !! umfnew/umfold
1044  real(RP) :: fbfrc
1045  real(RP), intent(out) :: dpthmx
1046  !---------------------------------------------------------------------------
1047 
1048  integer, parameter :: itr_max = 10000
1049 
1051  real(RP) :: pres300 ! pressure sfc-300hpa
1053  integer :: kk,nk ! kk is "do loop index" and nk is counter of "do loop "
1054  integer :: k_llfc ! k of LFC height haight of lowest
1056  integer :: n_uslcheck ! usl chek layer number
1057  real(RP) :: pres15 ! temporaly valuables pressure 15 hpa interval
1058  ! calculate mix tempreature (usl has 50mb or so)
1059  ! usl layer variable VV
1060  real(RP) :: theta_mix ! theta in usl layer
1061  real(RP) :: temp_mix ! temperature in usl layer
1062  real(RP) :: qv_mix ! water vaper mixing ratio in usl layer
1063  real(RP) :: emix ! saturate water vaper mixing ratio in usl layer
1064  !
1065  real(RP) :: temp_dew ! dew point temperature
1066  real(RP) :: templog ! log emix/aliq
1067  real(RP) :: deltaz ! dz used for interpolation of lcl
1068  real(RP) :: temp_env ! environment temperature
1069  real(RP) :: tempv_env ! environment temperature vertual
1070  real(RP) :: qvenv ! environment q vapor
1071  real(RP) :: w_cz ! velocity Kain(2004) c(z)
1072  real(RP) :: w_g ! velocity Kain(2004) wg
1073  real(RP) :: w_lcl ! vertical velocity @lcl
1074  real(RP) :: temp_lcl ! temperature @lcl
1075  real(RP) :: tempv_lcl ! vertal temperature @lcl
1076  real(RP) :: pres_lcl ! pressure @ lcl
1077  real(RP) :: dtvv ! delta Tvv Kain(2004)
1078  real(RP) :: dtrh ! trigger rh
1079  real(RP) :: radius ! cloud radius
1080  real(RP) :: dptt ! pressure thickness
1081  real(RP) :: dumfdp ! temp var
1082  real(RP) :: d_min
1083  ! calc in subroutin kf_updraft
1084  real(RP) :: umfnew,umfold ! from updraft
1085  integer :: k_lclm1 ! k_lcl -1
1086  integer :: k_start ! tempraly val
1087  integer :: k_check(ka) ! check layer index (because of 15mb interbal)
1088  integer :: n_check ! num of check
1089  integer :: n_layers ! num of USL layer
1090  integer :: nchm ! used shallow convection layer index
1091  integer :: nuchm ! (used shallow convection check)
1092  real(RP) :: CHMAX ! max cloud height used in shallow convection
1093  integer :: NNN ! internal work
1094  real(RP) :: CLDHGT(ka) ! used for shallow convection
1095  real(RP) :: dpthmin
1096  ! trigger 3 variables
1097  real(RP) :: rh_lcl
1098  real(RP) :: U00
1099  real(RP) :: DQSSDT
1100  real(RP) :: qs_lcl
1101 
1102  integer :: itr
1103  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1104  !! start cood
1105  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1106  pres300 = pres(ks) - depth_usl*100._rp ! pressure @ surface - 300 mb. maybe 700mb or so default depth_usl is 300hPa
1107  tempv(:) = temp(:)*(1._rp + 0.608_rp*qv(:)) ! vertual temperature
1108  ! search above 300 hPa index to "k_llfc"
1109  do kk = ks, ke
1110  if (pres(kk) >= pres300) k_llfc = kk
1111  end do
1112  ! usl(updraft sourcer layer) has interval for 15mb
1113  n_check = ks ! first layer
1114  k_check(ks) = ks ! first layer
1115  pres15 = pres(ks) - 15.e2_rp
1116 ! k_check = KS
1117  ! calc 15 hpa interval Num of layer(n_check) and index of layer(k_check)
1118  do kk = ks+1, k_llfc
1119  if(pres(kk) < pres15 ) then
1120  n_check = n_check + 1
1121  k_check(n_check) = kk
1122  pres15 = pres15 - 15.e2_rp
1123  end if
1124  end do
1125  ! main routine
1126  ! set initial
1127  fbfrc = 0._rp ! set inital 0
1128  dpthmin = 50.e2_rp ! set initialy USL layer depth (50 hPa)
1129  i_convflag = 2 ! no convection set initialy
1130  nuchm = ks-1 ! initialize (used shallow convection check) !like "0"
1131  n_uslcheck = ks-1 ! initialize index of usl check like "0"
1132  nchm = ks-1 ! initializelike "0"
1133  cldhgt(:) = 0._rp ! cloud hight is all 0
1134  do itr = 1, itr_max ! usl
1135  n_uslcheck = n_uslcheck + 1 ! nu in WRF
1136  ! calc shallow convection after all layer checked
1137  ! NOTE:
1138  ! This 'if block' is only used
1139  if (n_uslcheck > n_check) then ! if over potentially usl layer
1140  if (i_convflag == 1) then ! if schallow convection then
1141  chmax = 0._rp ! initialize Cloud Height Max
1142  nchm = ks-1 ! initialize index of Cloud Height Max "like 0"
1143  do kk = ks, n_check
1144  nnn = k_check(kk) ! index of checked layer (15 hpa interval layer)
1145  if (cldhgt(nnn) > chmax) then ! get max cloud height in shallow convection
1146  nchm = nnn
1147  nuchm = kk
1148  chmax = cldhgt(nnn)
1149  end if
1150  end do
1151  n_uslcheck = nuchm - 1
1152  fbfrc = 1._rp ! shallow convection is no precipitation
1153  cycle ! usl ! recalc updraft for shallow convection
1154  else ! not shallow convection then
1155  i_convflag = 2 ! no convection
1156  return
1157  end if ! end if schallow convection then
1158  end if ! end if over potentially usl layer
1159  k_lc = k_check(n_uslcheck)
1160  n_layers = 0 ! number of USL layer contain discretization layers
1161  dpthmx = 0._rp ! depth of USL (pressure [Pa])
1162  ! nk = k_lc - KS !< check k_lc -KS is not surface or bottom of ground
1163  nk = k_lc - 1
1164  ! calculate above k_lc layers depth of pressure
1165  ! usl layer depth is nessesary 50mb or more
1166  if ( nk + 1 < ks ) then
1167  if(kf_log) then
1168  if(io_l) write(io_fid_log,*) 'would go off bottom: cu_kf',pres(ks),k_lc,nk+1,n_uslcheck !,nk i,j
1169  end if
1170  else
1171  do ! serach USL layer index. USL layer is nessesally more 50hPa
1172  nk = nk +1
1173  if ( nk > ke ) then
1174  if(kf_log) then
1175  if(io_l) write(io_fid_log,*) 'would go off top: cu_kf'!,nk i,j
1176  end if
1177  exit
1178  end if
1179  dpthmx = dpthmx + deltap(nk) ! depth of pressure
1180  n_layers = n_layers + 1
1181  if (dpthmx > dpthmin) then
1182  exit
1183  end if
1184  end do
1185  end if
1186  if (dpthmx < dpthmin) then ! if dpthmx(USL layer depth) < dptmin(minimum USL layerdepth)
1187  i_convflag = 2
1188  return ! chenge at 3d version
1189  end if
1190  k_pbl = k_lc + n_layers -1
1191  ! initialize
1192  theta_mix = 0._rp; temp_mix = 0._rp; qv_mix = 0._rp; zmix = 0._rp;presmix = 0._rp
1193  ! clc mix variable (USL layer average valuable)
1194  do kk = k_lc, k_pbl
1195  temp_mix = temp_mix + deltap(kk)*temp(kk)
1196  qv_mix = qv_mix + deltap(kk)*qv(kk)
1197  zmix = zmix + deltap(kk)*z_kf(kk)
1198  presmix = presmix + deltap(kk)*pres(kk)
1199  end do
1200  ! mix tmp calculate
1201  temp_mix = temp_mix/dpthmx
1202  qv_mix = qv_mix/dpthmx
1203  presmix = presmix/dpthmx
1204  zmix = zmix/dpthmx
1205  emix = qv_mix*presmix/(0.622_rp + qv_mix) ! water vapor pressure
1206  ! calculate dewpoint and LCL temperature not use look up table
1207  ! LCL is Lifted condensation level
1208  ! this dewpoint formulation is Bolton's formuration (see Emanuel 1994 4.6.2)
1209  templog = log(emix/aliq)
1210  temp_dew = (cliq - dliq*templog)/(bliq - templog) ! dew point temperature
1211  ! temperature @ lcl setting need dewpoit
1212  ! this algolizm proposed by Davies-Jones (1983)
1213  temp_lcl = temp_dew - (0.212_rp + 1.571e-3_rp*(temp_dew - tem00) &
1214  - 4.36e-4_rp*(temp_mix - tem00))*(temp_mix - temp_dew) ! LCL temperature
1215  temp_lcl = min(temp_lcl,temp_mix)
1216  tempv_lcl = temp_lcl*(1._rp + 0.608_rp*qv_mix) ! LCL vertual temperature
1217  zlcl = zmix + (temp_lcl - temp_mix)/gdcp ! height of LCL
1218  ! index of z@lcl
1219  ! z(k_lclm1) < zlcl < z(k_lcl)
1220  do kk = k_lc, ke
1221  k_lcl = kk
1222  if( zlcl <= z_kf(kk) ) exit
1223  end do
1224  k_lclm1 = k_lcl - 1
1225  if (zlcl > z_kf(ke)) then
1226  i_convflag = 2 ! no convection
1227  return
1228  end if
1229  !! interpolate environment temperature and vapor at LCL
1230  deltaz = ( zlcl - z_kf(k_lclm1) )/( z_kf(k_lcl)- z_kf(k_lclm1 ) )
1231  temp_env = temp(k_lclm1) + ( temp(k_lcl) - temp(k_lclm1) )*deltaz
1232  qvenv = qv(k_lclm1) + ( qv(k_lcl) - qv(k_lclm1) )*deltaz
1233  tempv_env = temp_env*( 1._rp + 0.608_rp*qvenv )
1234  ! w_cz set review Kain(2004) eq.2
1235  ! w_g set
1236  ! dtlcl setting
1237  ! wg is
1238  if (zlcl < 2.e3_rp) then ! Kain 2004 eq.2
1239 ! w_cz = 0.02_RP*zlcl/2.e3_RP!
1240  w_cz = 1.e-5_rp*zlcl !
1241  else
1242  w_cz = 0.02_rp
1243  end if
1244  !! wg is iapproximate running mean grid resolved vertical velocity at the lcl(m/s)
1245  !! need w0avg
1246  !!wg = wg-c(z)
1247  w_g = (w0avg(k_lclm1) + (w0avg(k_lcl) - w0avg(k_lclm1))*deltaz )*deltax/25.e3_rp - w_cz ! need w0avg
1248  if ( w_g < 1.e-4_rp) then ! too small wg
1249  dtvv = 0._rp
1250  else
1251  dtvv = 4.64_rp*w_g**0.33_rp ! Kain (2004) eq.1 **1/3
1252  end if
1253  !! triggers will be make
1254  dtrh = 0._rp
1255  !! in WRF has 3 type of trigger function
1256  if ( trigger == 2) then ! new function none
1257  else if ( trigger == 3) then ! NO2007 trigger function
1258  !...for ETA model, give parcel an extra temperature perturbation based
1259  !...the threshold RH for condensation (U00)...
1260  ! as described in Narita and Ohmori (2007), 12th Mesoscale Conf.
1261  !
1262  !...for now, just assume U00=0.75...
1263  !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
1264  u00 = 0.75_rp
1265  if(u00 < 1._rp) then
1266  qs_lcl=qes(k_lclm1)+(qes(k_lcl)-qes(k_lclm1))*deltaz
1267  rh_lcl = qvenv/qs_lcl
1268  dqssdt = qv_mix*(cliq-bliq*dliq)/((temp_lcl-dliq)*(temp_lcl-dliq))
1269  if(rh_lcl >= 0.75_rp .and. rh_lcl <=0.95_rp)then
1270  dtrh = 0.25*(rh_lcl-0.75)*qv_mix/dqssdt
1271  elseif(rh_lcl > 0.95_rp) then
1272  dtrh = (1._rp/rh_lcl-1._rp)*qv_mix/dqssdt
1273  else
1274  dtrh = 0._rp
1275  end if
1276  end if
1277  end if
1278 !@================================================================================ ! check
1279  if (temp_lcl + dtvv + dtrh < temp_env) then ! kf triggerfucn dtrh is used @ NHM trigger func
1280  ! parcel is not bouyant
1281  ! cycle and check one more up layer(15 hPa )
1282  cycle
1283  else ! parcel is bouyant ,determin updraft
1284  ! theta_lclm1 is need
1285 
1286  ! calc updraft theta_E
1287  call envirtht(presmix,temp_mix,qv_mix,theta_eu(k_lclm1))
1288  !
1289  if (dtvv + dtrh > 1.e-4_rp ) then
1290  w_lcl = 1._rp + 0.5_rp*sqrt(2._rp*grav*(dtvv + dtrh)*500._rp/tempv_env)! Kain(2004) eq. 3??
1291  w_lcl = min(w_lcl,3._rp)
1292  else
1293  w_lcl = 1._rp
1294  end if
1295  k_let = k_lcl ! add k_let
1296  ! interpolate pressure
1297  pres_lcl = pres(k_lclm1) + (pres(k_lcl) - pres(k_lclm1))*deltaz
1298  ! denslcl = pres_lcl/(R*tempv_lcl)
1299  ! temp_lcl is already caliculated
1300  if (w_g < 0 ) then
1301  radius = 1000._rp
1302  elseif (w_g > 0.1_rp) then
1303  radius = 2000._rp
1304  else
1305  radius = 1000._rp + 1000._rp*w_g*10._rp ! wg/0.1 -> w_g*10
1306  end if
1307 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1308 !@ Compute updraft properties
1309 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1310  call cp_kf_updraft ( &
1311  ! [OUT]
1312  umf(:) ,& ! upward mass flux
1313  umflcl ,& ! updraft mass flux@LCL
1314  upent(:) ,& ! updraft entrainment
1315  updet(:) ,& ! updraft detrainment
1316  umfnewdold(:) ,& ! ratio of umfnew/umfold
1317  umfnew,umfold ,& !
1318  temp_u(:) ,& ! updraft temperature internal work??
1319  theta_eu(:) ,& ! updraft theta_E
1320  theta_ee(:) ,& ! environment theta_E inout
1321  cloudhight ,& ! cloud depth
1322  totalprcp ,& ! total amount of precipitation and snow
1323  cloudtop ,& ! cloud top hight
1324  qv_u(:) ,& ! updraft qv
1325  qc(:) ,& ! cloud water
1326  qi(:) ,& ! cloud ice
1327  qrout(:) ,& ! rain
1328  qsout(:) ,& ! snow
1329  qvdet(:) ,& ! detrainment of qv
1330  qcdet(:) ,& ! detrainment of qc
1331  qidet(:) ,& ! detrainment of qc
1332  cape ,& ! cape !
1333  flux_qr(:) ,& ! flux of qr
1334  flux_qs(:) ,& ! flux of qi
1335  k_top ,& ! index of cloud top
1336  k_let ,& ! index of below detrain layer
1337  ! [IN]
1338  dz_kf(:),z_kf(:) ,&
1339  w_lcl ,& ! vertical velocity of LCL layer
1340  temp_lcl ,& ! temperature of LCL layer
1341  tempv_lcl ,& ! vertural temperature of LCL layer
1342  pres_lcl ,& ! pressure of LCL layer
1343  qv_mix ,& ! vapor of USL layer
1344  qv(:) ,& ! vapor
1345  temp(:) ,& ! temperature
1346  tempv_env ,& ! lcl temprature
1347  zlcl ,& ! z[m]@LCL will be intent in
1348  pres(:) ,& ! preassure
1349  deltap(:) ,& ! dp
1350  radius ,& ! cloud radius
1351  dpthmx ,& ! pressure of depth of USL layer
1352  k_lcl ,& ! index of LCL layer
1353  k_pbl & ! index of USL layer
1354  )
1355 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1356 !@ Compute updraft properties
1357 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1358 
1359  !! temporaly cloud hight
1360  !! for shallow convection
1361  cldhgt(k_lc) = cloudhight
1362  !! minimum cloud hight for deep convection
1363  !! Kain (2004) eq.7
1364  if(temp_lcl > 293._rp) then
1365  d_min = 4.e3_rp
1366  elseif (temp_lcl < 293._rp .and. temp_lcl > 273._rp) then
1367  d_min = 2.e3_rp + 100._rp*(temp_lcl - tem00)
1368  else
1369  d_min = 2.e3_rp
1370  end if
1371 
1372  !! convection type check
1373  if ( k_top <= k_lcl .or. &
1374  k_top <= k_pbl .or. &
1375  k_let+1 <= k_pbl ) then ! no convection allowd
1376  cloudhight = 0._rp
1377  cldhgt(k_lc) = cloudhight ! 0
1378  do kk = k_lclm1,k_top
1379  umf(kk) = 0._rp
1380  upent(kk) = 0._rp
1381  updet(kk) = 0._rp
1382  qcdet(kk) = 0._rp
1383  qidet(kk) = 0._rp
1384  flux_qr(kk) = 0._rp
1385  flux_qs(kk) = 0._rp
1386  end do
1387  elseif (cloudhight > d_min .and. cape > 1._rp) then ! shallow convection
1388  i_convflag = 0 ! deep convection
1389  exit ! exit usl loop
1390  else ! shallow convection
1391  !! shallow convection after determin
1392  !! remember this layer (by virture of non-zero CLDHGT) as potential shallo-cloudlayer
1393  !! detern shallow convection is after checked all potentialy USL layer
1394  !! when no deep convection is found (this is deterned first if block in do loop)
1395  i_convflag = 1
1396  if(n_uslcheck == nuchm) then !
1397  exit ! exit usl loop
1398  else
1399  do kk = k_lclm1,k_top
1400  umf(kk) = 0._rp
1401  upent(kk) = 0._rp
1402  updet(kk) = 0._rp
1403  qcdet(kk) = 0._rp
1404  qidet(kk) = 0._rp
1405  flux_qr(kk) = 0._rp
1406  flux_qs(kk) = 0._rp
1407  end do
1408  end if ! if(n_uslcheck == NUCHM) then
1409  end if ! convection type
1410  end if ! triggeer
1411  end do ! usl
1412  if ( itr .ge. itr_max ) then
1413  write(*,*) 'xxx iteration max count was reached in the USL loop in the KF scheme'
1414  call prc_mpistop
1415  end if
1416 
1417 
1418  if (i_convflag == 1) then ! shallow convection
1419  k_start = max(k_pbl,k_lcl)
1420  k_let = k_start
1421  end if
1422  !
1423  !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1424  ! THIS LEVEL...
1425  !
1426  if (k_let == k_top) then
1427  updet(k_top) = umf(k_top) + updet(k_top) - upent(k_top)
1428  qcdet(k_top) = qc(k_top)*updet(k_top)*umfnew/umfold
1429  qidet(k_top) = qi(k_top)*updet(k_top)*umfnew/umfold
1430  upent(k_top) = 0._rp
1431  umf(k_top) = 0._rp
1432  else
1433  !! begin total detrainment at the level above the k_let
1434  dptt = 0._rp
1435  do kk = k_let+1,k_top
1436  dptt = dptt + deltap(kk)
1437  end do
1438  dumfdp = umf(k_let)/dptt
1439  !!
1440  !!... adjust mass flux profiles, detrainment rates, and precipitation fall
1441  !!... rates to reflect the linear decrease in mass flux between the k_let and
1442  !!
1443  !!...entrainment is allowed at every level except for K_TOP, so disallow
1444  !!...entrainment at k_TOP and adjust entrainment rates between k_LET and k_TOP
1445  !!...so the the dilution factor due to entrainment is not changed but
1446  !!...the actual entrainment rate will change due due forced total
1447  !!...detrainment in this layer...
1448  do kk = k_let+1,k_top
1449  if(kk == k_top) then
1450  updet(kk) = umf(kk-1)
1451  upent(kk) = 0._rp
1452  qcdet(kk) = updet(kk)*qc(kk)*umfnewdold(kk)
1453  qidet(kk) = updet(kk)*qi(kk)*umfnewdold(kk)
1454  else
1455  umf(kk) = umf(kk-1) - deltap(kk)*dumfdp
1456  upent(kk) = umf(kk)*(1._rp - 1._rp/umfnewdold(kk))
1457  updet(kk) = umf(kk-1) - umf(kk) + upent(kk)
1458  qcdet(kk) = updet(kk)*qc(kk)*umfnewdold(kk)
1459  qidet(kk) = updet(kk)*qi(kk)*umfnewdold(kk)
1460  end if
1461  if (kk >= k_let+2) then
1462  totalprcp = totalprcp - flux_qr(kk) - flux_qs(kk)
1463  flux_qr(kk) = umf(kk-1)*qrout(kk)
1464  flux_qs(kk) = umf(kk-1)*qsout(kk)
1465  totalprcp = totalprcp + flux_qr(kk) + flux_qs(kk)
1466  end if
1467  !
1468  end do
1469 
1470  end if ! let== ktop
1471 
1472  !! initialize some arrays below cloud base and above cloud top
1473  !! below cloud base
1474  !! melt layer setting
1475  do kk = ks,k_top
1476  if(temp(kk) > tem00) k_ml = kk !! melt layer
1477  end do
1478  !
1479  do kk = ks,k_lclm1
1480  !!
1481  if(kk >= k_lc) then
1482  !
1483  if (kk == k_lc) then ! if bototm of USL(pbl)
1484  upent(kk) = umflcl*deltap(kk)/dpthmx
1485  umf(kk) = umflcl*deltap(kk)/dpthmx
1486  elseif (kk <= k_pbl) then ! if in USL(pbl)
1487  upent(kk) = umflcl*deltap(kk)/dpthmx
1488  umf(kk) = umf(kk-1) + upent(kk) !! assume no detrain
1489  else
1490  upent(kk) = 0._rp
1491  umf(kk) = umflcl
1492  end if
1493  !
1494  temp_u(kk) = temp_mix + (z_kf(kk) - zmix)*gdcp ! layner assumption of temp_u
1495  qv_u(kk) = qv_mix
1496  ! wu(kk) = w_lcl no use @wrf
1497  else ! below USL layer no updraft
1498  temp_u(kk) = 0._rp
1499  qv_u(kk) = 0._rp
1500  umf(kk) = 0._rp
1501  upent(kk) = 0._rp
1502  end if
1503  !
1504  updet(kk) = 0._rp
1505  qvdet(kk) = 0._rp
1506  qc(kk) = 0._rp
1507  qi(kk) = 0._rp
1508  qrout(kk) = 0._rp
1509  qsout(kk) = 0._rp
1510  flux_qr(kk) = 0._rp
1511  flux_qs(kk) = 0._rp
1512  qcdet(kk) = 0._rp
1513  qidet(kk) = 0._rp
1514  !!calc theta_E environment
1515  call envirtht(pres(kk),temp(kk),qv(kk),theta_ee(kk))
1516  !!
1517  end do
1518 
1519  ! define variables above cloud top
1520  do kk = k_top+1,ke
1521  umf(kk) = 0._rp
1522  upent(kk) = 0._rp
1523  updet(kk) = 0._rp
1524  qvdet(kk) = 0._rp
1525  qc(kk) = 0._rp
1526  qi(kk) = 0._rp
1527  qrout(kk) = 0._rp
1528  qsout(kk) = 0._rp
1529  flux_qr(kk) = 0._rp
1530  flux_qs(kk) = 0._rp
1531  qcdet(kk) = 0._rp
1532  qidet(kk) = 0._rp
1533  end do
1534  do kk = k_top+2,ke
1535  temp_u(kk) = 0._rp
1536  qv_u(kk) = 0._rp
1537  end do
1538 
1539  return
1540  end subroutine cp_kf_trigger
1541  !!-------------------------------------------------------------------------------
1542  !! updraft
1543  !!-------------------------------------------------------------------------------
1544  subroutine cp_kf_updraft (& ! 1d model updraft
1545  ![OUT]
1546  umf ,& ! upward mass flux(updraft)
1547  umflcl ,& ! upward mass flux@LCL
1548  upent ,& ! updraft entrainment
1549  updet ,& ! updraft detrainment
1550  umfnewdold ,& ! UMF side entrainment/detrainment calculation before(old) or after(new)
1551  umfnew,umfold ,& ! updraft
1552  temp_u ,& ! updraft temperature internal work??
1553  theta_eu ,& ! updraft theta_E
1554  theta_ee ,& ! thera_e @env
1555  cloudhight ,& ! cloud depth [m]
1556  totalprcp ,& ! total amount of precipitation and snow
1557  cloudtop ,& ! cloud top hight [m]
1558  qv_u ,& ! updraft qv
1559  qc ,& ! cloud water
1560  qi ,& ! cloud ice
1561  qrout ,& ! rain
1562  qsout ,& ! snow
1563  qvdet ,& ! detrainment of qv
1564  qcdet ,& ! detrainment of qc
1565  qidet ,& ! detrainment of qc
1566  cape ,& ! cape
1567  flux_qr ,& ! flux of qr
1568  flux_qs ,& ! flux of qi
1569  k_top ,& ! index of cloud top layer
1570  k_let ,& ! index of below detrain layer
1571  ! [IN]
1572  dz_kf,z_kf ,& ! delataz and z height
1573  w_lcl ,& ! vertical velocity of LCL layer
1574  temp_lcl ,& ! temperature of LCL layer
1575  tempv_lcl ,& ! vertural temperature of LCL layer
1576  pres_lcl ,& ! pressure of LCL layer
1577  qv_mix ,& ! water vapor of USL layer
1578  qv ,& ! water vapor
1579  temp ,& ! temperature
1580  tempv_env ,& ! temperature @ environmetnt
1581  zlcl ,& ! z[m]@LCL will be intent in
1582  pres ,& ! preassure [hPa]
1583  deltap ,& ! dp [pa]
1584  radius ,& ! cloud radius [m]
1585  dpthmx ,& ! pressure of depth of USL layer
1586  k_lcl ,& ! index of LCL layer
1587  k_pbl ) ! index of USL layer
1588  use scale_precision
1589  use scale_grid_index
1590  use scale_const,only :&
1591  cp => const_cpdry , &
1592  pre00 => const_pre00 , &
1593  grav => const_grav , &
1594  r => const_rdry
1595  implicit none
1596  ![OUT]
1597  real(RP),intent(out) :: umf(ka) ! upward mass flux
1598  real(RP),intent(out) :: umflcl ! upward mass flux @lcl
1599  real(RP),intent(out) :: upent(ka) ! updraft entrainment
1600  real(RP),intent(out) :: updet(ka) ! updraft detrainment
1601  real(RP),intent(out) :: umfnewdold(ka) ! ratio of below umfnew/umfold
1602  real(RP),intent(out) :: umfnew,umfold ! UMF side entrainment/detrainment calculation before(old) or after(new)
1603  real(RP),intent(out) :: temp_u(ka) ! updraft temperature internal work
1604  real(RP),intent(out) :: theta_ee(ka)
1605  real(RP),intent(inout) :: theta_eu(ka)
1606  real(RP),intent(out) :: cloudhight ! cloud depth
1607  real(RP),intent(out) :: totalprcp ! precipitation
1608  real(RP),intent(out) :: cloudtop ! cloud top height
1609  real(RP),intent(out) :: qv_u(ka) ! updraft water vape mixing ratio
1610  real(RP),intent(out) :: qc(ka) ! liquit QLIQ in WRF
1611  real(RP),intent(out) :: qi(ka) ! QICE in WRF
1612  real(RP),intent(out) :: qrout(ka) ! QLQOUT in WRF
1613  real(RP),intent(out) :: qsout(ka) ! QICOUT in WRF
1614  real(RP),intent(out) :: qvdet(ka) ! will be intent out
1615  real(RP),intent(out) :: qcdet(ka) ! detrainment of qc
1616  real(RP),intent(out) :: qidet(ka) ! detrainment of qc
1617  real(RP),intent(out) :: cape ! CAPE
1618  real(RP),intent(out) :: flux_qr(ka) ! ratio of generation of liquit fallout(RAIN)
1619  real(RP),intent(out) :: flux_qs(ka) ! ratio of generation of ice fallout(SNOW)
1620  integer, intent(out) :: k_top ! top of convection layer index
1621  integer, intent(inout) :: k_let ! top of convection layer not detrain only layer
1622  ! [IN]
1623  real(RP),intent(in) :: dz_kf(ka),z_kf(ka)
1624  ! lcl layer valiable
1625  real(RP),intent(in) :: w_lcl
1626  real(RP),intent(in) :: temp_lcl
1627  real(RP),intent(in) :: tempv_lcl
1628  real(RP),intent(in) :: pres_lcl
1629  real(RP),intent(in) :: qv_mix
1630  real(RP),intent(in) :: qv(ka)
1631  real(RP),intent(in) :: temp(ka)
1632  real(RP),intent(in) :: tempv_env ! environment vertual temperature
1633  real(RP),intent(in) :: zlcl ! z[m]@LCL will be intent in
1634  real(RP),intent(in) :: pres(ka)
1635  real(RP),intent(in) :: deltap(ka)
1636  real(RP),intent(in) :: radius
1637  real(RP),intent(in) :: dpthmx ! pressure of depth of USL layer
1638  integer ,intent(in) :: k_lcl
1639  integer ,intent(in) :: k_pbl
1640  ![INTERNAL WORK]
1641  integer :: kk,kkp1 ! kk : do loop ,kkp1: kk+1
1642  integer :: k_lclm1 ! k_lcl -1
1643  real(RP) :: tempv_u(ka) ! updraft vertual temperature internalwork
1644  real(RP) :: tempvq_u(ka) ! temperature vertial for updraft ,qv, qc, qi
1645  real(RP) :: denslcl ! density @LCL
1646  real(RP) :: ee1,ud1, ee2,ud2 ! entrainment and detrainment calc valiables
1647  real(RP) :: f_eq(ka) ! prof5 variable
1648  real(RP) :: f_mix1,f_mix2 ! factor of mixed fraction
1649  real(RP) :: REI,DILBE ! REI KF(1990) Eq.1 , DILBE is tempvar for calc cape
1650  real(RP) :: qcnew,qinew ! qcnew is qc new var , qinew is qi newver
1651  real(RP) :: qfrz
1652  real(RP) :: f_frozen1
1653  real(RP) :: temptmp ! temporaly temperature
1654  real(RP) :: temptmp_ice ! temporaly temperature for ice or liquid face calculateion
1655  real(RP) :: tempv(ka) ! virtual temperature [K]
1656  real(RP) :: wtw ! w**2
1657  real(RP) :: boeff ! bouyancy effect
1658  real(RP) :: boterm ! bouyancy term for calc vertical vilocity
1659  real(RP) :: dztmp
1660  real(RP) :: entterm ! entrainment term for calc vertical vilocity
1661  real(RP) :: theta_tmp ! tmporaly temperature
1662  real(RP) :: wu(ka) ! vertical velocity of updraft
1663  real(RP) :: qvtmp, qctmp, qitmp
1664  real(RP) :: temp_u95, temp_u10 ! temporaly Temperature value use determin Mixed Fraction
1665  real(RP),parameter :: temp_frzT = 268.16_rp
1666  real(RP),parameter :: temp_frzB = 248.16_rp
1670  tempv = temp*(1._rp + 0.608_rp*qv)
1671  ! initial updraft mass flux
1672  umfnewdold(:) = 1._rp
1673  k_lclm1 = k_lcl - 1
1674  wtw = w_lcl*w_lcl
1675  denslcl = pres_lcl/(r*tempv_lcl)
1676  umf(k_lclm1) = denslcl*1.e-2_rp*deltax*deltax ! (A0)
1677  umflcl = umf(k_lclm1)
1678  umfold = umflcl
1679  umfnew = umfold
1680  ! initial vas setup
1681  temp_u(k_lclm1) = temp_lcl
1682  tempv_u(k_lclm1) = tempv_lcl
1683  qv_u(k_lclm1) = qv_mix
1684  ! initialize
1685  ! only qv exist other water variables is none
1686  f_eq(k_lclm1) = 0._rp
1687  upent = 0._rp
1688  qc(k_lclm1) = 0._rp
1689  qi(k_lclm1) = 0._rp
1690  qrout(k_lclm1) = 0._rp
1691  qsout(k_lclm1) = 0._rp
1692  qcdet(k_lclm1) = 0._rp
1693  qidet(k_lclm1) = 0._rp
1694  flux_qr = 0._rp
1695  flux_qs = 0._rp
1696  totalprcp = 0._rp
1697  !
1698  ee1 = 1._rp
1699  ud1 = 0._rp
1700  rei = 0._rp
1701  dilbe = 0._rp
1702  cape = 0._rp
1703  !! temptmp is used during clculation of the linearglaciation
1704  !! process. it is initially set to the temperature at which freezing
1705  !! is specified to begin. within the glaciation
1706  !! interval, it is set equal to the updraft temp at the
1707  !! previous model level...
1708  temptmp_ice = temp_frzt
1709 
1712  do kk = k_lclm1,ke-1 ! up_main original(wrf cood K is k_lclm1)
1713  kkp1 = kk + 1
1714  ! temporaly use below layer valuables
1715  f_frozen1 = 0._rp ! frozen rate initialize this variable 0 (water)to 1(ice)
1716  temp_u(kkp1) = temp(kkp1) ! up parcel temperature assumed env temperature
1717  theta_eu(kkp1) = theta_eu(kk) ! up parcel theta_E
1718  qv_u(kkp1) = qv_u(kk) ! up parcel vapor
1719  qc(kkp1) = qc(kk) ! up parcel water
1720  qi(kkp1) = qi(kk) ! up parcel ice
1721  !!# calqulate satulation and cleate q and qic
1722  !!# and calc updraft pacel temperature
1723  !!# it is use determinnent of frozn or not
1724  !!#
1725  call tpmix2(pres(kkp1),theta_eu(kkp1),temp_u(kkp1),qv_u(kkp1),qc(kkp1),qi(kkp1),qcnew,qinew)
1726  !!# check to see if updraft temperature is avove the temperature at which
1727  !!# glaciation is assumed to initiate. if it is, calculate the
1728  !!# fraction of remainning liquid water to freeze... temp_frzT is the
1729  !!# temperature at which freezing begins, temp_frzB isthe temperature below which all
1730  !!# liquid water is frozen at each level...
1731  if (temp_u(kkp1) <= temp_frzt ) then ! if frozen temperature then ice is made
1732  if (temp_u(kkp1) > temp_frzb) then
1733  if (temptmp_ice > temp_frzt) temptmp_ice = temp_frzt
1734  !!# liner assumption determin frozen ratio
1735  f_frozen1 = (temptmp_ice - temp_u(kkp1))/(temptmp_ice - temp_frzb)
1736  else
1737  f_frozen1 = 1._rp ! all water is frozen
1738  endif
1739  temptmp_ice = temp_u(kkp1)
1740  !!# calc how much ice is a layer ?
1741  qfrz = (qc(kkp1) + qcnew)*f_frozen1 ! all ice
1742  qinew = qinew + qcnew*f_frozen1 ! new create ice
1743  qcnew = qcnew - qcnew*f_frozen1 ! new create liquit
1744  qi(kkp1) = qi(kkp1) + qc(kkp1)*f_frozen1 ! ice old + convert liquit to ice
1745  qc(kkp1) = qc(kkp1) - qc(kkp1)*f_frozen1 ! liquit old - convet liquit to ice
1746  !!# calculate effect of freezing
1747  !!# and determin new create frozen
1748  call dtfrznew(temp_u(kkp1), pres(kkp1),theta_eu(kkp1),qv_u(kkp1),qfrz,qi(kkp1) )
1749  end if
1750  tempv_u(kkp1) = temp_u(kkp1)*(1._rp + 0.608_rp*qv_u(kkp1)) ! updraft vertual temperature
1751  ! calc bouyancy term for verticl velocity
1752  if (kk == k_lclm1) then !! lcl layer exist between kk and kk+1 layer then use interporate value
1753  boeff = (tempv_lcl + tempv_u(kkp1))/(tempv_env + tempv(kkp1)) - 1._rp
1754  boterm = 2._rp*(z_kf(kkp1) - zlcl)*grav*boeff/1.5_rp
1755  dztmp = z_kf(kkp1) - zlcl
1756  else
1757  boeff = (tempv_u(kk) + tempv_u(kkp1))/(tempv(kk) + tempv(kkp1)) - 1._rp
1758  boterm = 2._rp*(dz_kf(kk) )*grav*boeff/1.5_rp
1759  dztmp = dz_kf(kk)
1760  end if
1761  !!# entrainment term
1762  entterm = 2._rp*rei*wtw/umfold
1763  !! calc precp(qr) , snow(qr) and vertical virocity
1764  call p_precipitation(&
1765  qc(kkp1), qi(kkp1), wtw, dztmp, boterm, entterm, qcnew, qinew, qrout(kkp1), &
1766  qsout(kkp1),grav )
1767  !! precipitation_OC1973 or precipitation_Kessler
1768  !! in Ogura and Cho (1973), 'RATE' is nessesary but this rate is parameter in subroutine "precipitation_OC1973"
1769  !! RATE is tuning parameter
1770  !! read Kain 1990
1771  if (wtw < 1.e-3_rp ) then ! if no vertical velocity
1772  exit ! end calc updraft
1773  else
1774  wu(kkp1) = sqrt(wtw)
1775  end if
1776  !!# calc tehta_e in environment to entrain into updraft
1777  call envirtht(pres(kkp1),temp(kkp1),qv(kkp1),theta_ee(kkp1))
1778  ! rei is the rate of environment inflow
1779  rei = umflcl*deltap(kkp1)*0.03_rp/radius !!# Kain 1990 eq.1 ;Kain 2004 eq.5
1780  !!--------------------------------------------------
1781  !! calc cape
1782  !!--------------------------------------------------
1783  tempvq_u(kkp1) = temp_u(kkp1)*(1._rp + 0.608_rp*qv_u(kkp1) - qc(kkp1) - qi(kkp1))
1784  if (kk == k_lclm1) then!! lcl layer exist between kk and kk+1 then use interporate value
1785  dilbe = ((tempv_lcl + tempvq_u(kkp1))/(tempv_env + tempv(kkp1)) - 1._rp)*dztmp
1786  else
1787  dilbe = ((tempvq_u(kk) + tempvq_u(kkp1))/(tempv(kk) + tempv(kkp1)) - 1._rp)*dztmp
1788  end if
1789  if(dilbe > 0._rp) cape = cape + dilbe*grav
1790  !!--------------------------------------------------
1791  !! if cloud parcels are virtually colder then the entrainment , minimal
1792  !! entrainment 0.5*rei is imposed...
1793  !! read Kain 2004
1794  !!--------------------------------------------------
1795  !! calc entrainment/detrainment
1796  if(tempvq_u(kkp1) <= tempv(kkp1)) then ! if entrain and detrain
1797  ! original KF90 no entrainment allow
1798  ee2 = 0.5_rp ! Kain (2004) eq.4
1799  ud2 = 1._rp
1800  f_eq(kkp1) = 0._rp
1801  else
1802  k_let = kkp1
1803  temptmp = tempvq_u(kkp1)
1804  !! determine teh critical mixed fraction of updraft and environmental air ...
1805  !! if few mix moisture air
1806  f_mix1 = 0.95_rp
1807  f_mix2 = 1._rp - f_mix1
1808  theta_tmp = f_mix1*theta_ee(kkp1) + f_mix2*theta_eu(kkp1)
1809  qvtmp = f_mix1*qv(kkp1) + f_mix2*qv_u(kkp1)
1810  qctmp = f_mix2*qc(kkp1)
1811  qitmp = f_mix2*qi(kkp1)
1812  !! need only temptmp because calc bouyancy
1813  call tpmix2(pres(kkp1),theta_tmp,temptmp,qvtmp,qctmp,qitmp,qcnew,qinew)
1814  !! qinew and qcnew is damy valuavle(not use )
1815  temp_u95 = temptmp*(1._rp + 0.608_rp*qvtmp - qctmp - qitmp)
1816  ! TU95 in old coad
1817  if ( temp_u95 > tempv(kkp1)) then ! few mix but bouyant then ! if95
1818  ee2 = 1._rp ! rate of entrain is 1 -> all entrain
1819  ud2 = 0._rp
1820  f_eq(kkp1) = 1._rp
1821  else
1822  f_mix1 = 0.10_rp
1823  f_mix2 = 1._rp - f_mix1
1824  theta_tmp = f_mix1*theta_ee(kkp1) + f_mix2*theta_eu(kkp1)
1825  qvtmp = f_mix1*qv(kkp1) + f_mix2*qv_u(kkp1)
1826  qctmp = f_mix2*qc(kkp1)
1827  qitmp = f_mix2*qi(kkp1)
1828  !! need only temptmp because calc bouyancy
1829  call tpmix2(pres(kkp1),theta_tmp,temptmp,qvtmp,qctmp,qitmp,qcnew,qinew)
1830  !! qinew and qcnew is damy valuavle(not use )
1831  temp_u10 = temptmp*(1._rp + 0.608_rp*qvtmp - qctmp - qitmp)
1832  if (abs(temp_u10 - tempvq_u(kkp1)) < 1.e-3_rp ) then !if10%
1833  ee2 = 1._rp ! all entrain
1834  ud2 = 0._rp
1835  f_eq(kkp1) = 1._rp
1836  else
1837  f_eq(kkp1) = (tempv(kkp1) - tempvq_u(kkp1))*f_mix1 &
1838  & /(temp_u10 - tempvq_u(kkp1))
1839  f_eq(kkp1) = max(0._rp,f_eq(kkp1) )
1840  f_eq(kkp1) = min(1._rp,f_eq(kkp1) )
1841  if (f_eq(kkp1) == 1._rp) then ! f_eq
1842  ee2 = 1._rp
1843  ud2 = 0._rp
1844  elseif (f_eq(kkp1) == 0._rp) then
1845  ee2 = 0._rp
1846  ud2 = 1._rp
1847  else
1848  !!# subroutine prof5 integrates over the gaussian dist to determine
1849  !!# the fractional entrainment and detrainment rates...
1850  call prof5(f_eq(kkp1),ee2,ud2)
1851  end if ! end of f_iq
1852  end if ! end of if10%
1853  end if ! end of if95%
1854  end if! end of entrain/detrain
1855  ee2 = max(ee2,0.5_rp)
1856  ud2 = 1.5_rp*ud2
1857  !
1858  upent(kkp1) = 0.5*rei*(ee1 + ee2)
1859  updet(kkp1) = 0.5*rei*(ud1 + ud2)
1860  !! if the calculated updraft detrainment rate is greater then the total
1861  !! updraft mass flux(umf), all cloud mass detrains
1862  if (umf(kk) - updet(kkp1) < 10._rp) then ! why 10.???
1863  !! if the calculated detrained mass flux is greater than totla upd mass
1864  !! flux, impose total detrainment of updraft mass at the previous model level..
1865  !! first, correct cape calculation if neede
1866  if (dilbe > 0._rp) then
1867  cape = cape - dilbe*grav
1868  end if
1869  k_let = kk ! then k_let = k_top
1870  exit
1871  else
1872  ee1 = ee2
1873  ud1 = ud2
1874  umfold = umf(kk) - updet(kkp1)
1875  umfnew = umfold + upent(kkp1)
1876  umf(kkp1) = umfnew
1877  umfnewdold(kkp1) = umfnew/umfold
1878  qcdet(kkp1) = qc(kkp1)*updet(kkp1)
1879  qidet(kkp1) = qi(kkp1)*updet(kkp1)
1880  qvdet(kkp1) = qv_u(kkp1)
1881  !! below layer updraft qv and entrain q /new updraft massflux
1882  qv_u(kkp1) = (umfold*qv_u(kkp1) + upent(kkp1)*qv(kkp1))/umfnew
1883  theta_eu(kkp1) = (umfold*theta_eu(kkp1) + upent(kkp1)*theta_ee(kkp1))/umfnew
1884  ! in WRF
1885  ! PPTLIQ(KKP11) = qlqout(KKP11)*UMF(KKP1)
1886  qc(kkp1) = qc(kkp1)*umfold/umfnew
1887  qi(kkp1) = qi(kkp1)*umfold/umfnew
1888  ! flux_qr is ratio of generation of liquit fallout(RAIN)
1889  ! flux_qi is ratio of generation of ice fallout(SNOW)
1890  flux_qr(kkp1) = qrout(kkp1)*umf(kk)
1891  flux_qs(kkp1) = qsout(kkp1)*umf(kk)
1892  totalprcp = totalprcp + flux_qr(kkp1) + flux_qs(kkp1) ! total prcp vars
1893  !! if below cloud base then
1894  !! updraft entrainment is umf@LCL*dp/depth
1895  if(kkp1 <= k_pbl ) upent(kkp1) = upent(kkp1) + umflcl*deltap(kkp1)/dpthmx
1896  end if
1897  end do ! this loop is exit at w=0.
1898  k_top = kk ! vertical coodinate index @ w=0
1899  ! cloud height
1900  cloudhight = z_kf(k_top) - zlcl
1901  cloudtop = z_kf(k_top)
1902  return
1903  end subroutine cp_kf_updraft
1904  !!------------------------------------------------------------------------------
1905  !! Down draft
1906  !!------------------------------------------------------------------------------
1907  subroutine cp_kf_downdraft ( &
1908  ![IN]
1909  dz_kf,z_kf ,& !
1910  u ,& ! x-direction wind
1911  v ,& ! meridional wind
1912  zlcl ,& ! lcl_hight [m]
1913  rh ,& ! relative humidity initial make kf_init
1914  deltap ,& ! delta P
1915  pres ,& ! pressure [Pa]
1916  qv ,& ! watervapor
1917  ems ,& !
1918  emsd ,& !
1919  theta_ee ,& ! environment equivalent theta
1920  umf ,& ! upward mass flux
1921  totalprcp ,& ! total precipitation
1922  flux_qs ,& !
1923  tempv ,&
1924  I_convflag ,& ! convection flag
1925  ! index
1926  k_lcl ,&
1927  k_ml ,&
1928  k_top ,&
1929  k_pbl ,&
1930  k_let ,&
1931  k_lc ,&
1932  ![OUT]
1933  wspd ,& ! wind speed 1 k_lcl, 2 k_z5,3 k_top
1934  dmf ,& ! down
1935  downent ,&
1936  downdet ,&
1937  theta_d ,&
1938  qv_d ,&
1939  prcp_flux ,&
1940  k_lfs ,&
1941  CPR ,&
1942  tder)
1943  use scale_precision
1944  use scale_grid_index
1945  use scale_const,only :&
1946  cp => const_cpdry , &
1947  pre00 => const_pre00 , &
1948  emelt => const_emelt , &
1949  grav => const_grav , &
1950  r => const_rdry
1951  use scale_atmos_saturation ,only :&
1952  atmos_saturation_psat_liq
1953  use scale_process, only: &
1954  prc_mpistop
1955  implicit none
1956  ! [IN]
1957  ! from init
1958  real(RP),intent(in) :: dz_kf(ka),z_kf(ka)
1959  real(RP),intent(in) :: u(ka) ! x velocity
1960  real(RP),intent(in) :: v(ka) ! y velocity
1961  real(RP),intent(in) :: rh(ka) ! ! initial make kf_init
1962  real(RP),intent(in) :: deltap(ka) ! delta pressure
1963  real(RP),intent(in) :: pres(ka) ! pressure
1964  real(RP),intent(in) :: qv(ka) ! watervapor
1965  real(RP),intent(in) :: ems(ka) ! dp*(deltax)**2/g
1966  real(RP),intent(in) :: emsd(ka) ! 1._RP/ems
1967  ! from kf_triger
1968  real(RP),intent(in) :: zlcl ! lcl_hight
1969  real(RP),intent(in) :: umf(ka) ! upward mass flux
1970  real(RP),intent(in) :: totalprcp ! in !total precipitation
1971  real(RP),intent(in) :: flux_qs(ka) ! snow fall
1972  real(RP),intent(in) :: tempv(ka) ! virtual env temperature internal work
1973  real(RP),intent(in) :: theta_ee(ka) ! in ! environment equivalent theta
1974  integer :: k_z5 ! internalwork 500hPa layer index
1975  ! from kf_triger
1976  integer,intent(in) :: I_convflag ! index of convection
1977  integer,intent(in) :: k_lcl ! index of lcl layer
1978  integer,intent(in) :: k_top ! index of cloud top layer
1979  integer,intent(in) :: k_pbl ! index of USL layer top
1980  integer,intent(in) :: k_let ! index of updraft only detrainment
1981  integer,intent(in) :: k_lc ! index of USL layer bottom
1982  integer,intent(in) :: k_ml ! index of melt layer
1983  ! [OUT]
1984  real(RP),intent(out) :: wspd(3) ! wind speed 1 k_lcl, 2 k_z5,3 k_top
1985  real(RP),intent(out) :: dmf(ka) ! downdraft massflux
1986  real(RP),intent(out) :: downent(ka) ! downdraft entrainment
1987  real(RP),intent(out) :: downdet(ka) ! downdraft detrainment
1988  real(RP),intent(out) :: theta_d(ka) ! downdraft theta
1989  real(RP),intent(out) :: qv_d(ka) ! vapor of downdraft
1990  real(RP),intent(out) :: prcp_flux ! out! precpitation flux
1991  real(RP),intent(out) :: tder ! temperature change from evap
1992  real(RP),intent(out) :: CPR ! all precipitation before consider cloud bottom evaporation
1993  integer,intent(out) :: k_lfs ! LFS layer index
1994  ! [INTERNAL WORK]
1995  integer :: kk, kkp1 ! do loop variable
1996  real(RP) :: shsign ! frip variable tmp use
1997  real(RP) :: vws ! vertical wind shear
1998  real(RP) :: pef ! precipitation efficiency
1999  real(RP) :: pef_wind ! precipitation efficiency from wind shear
2000  real(RP) :: pef_cloud ! precipitation efficiency from cloud hight
2001  real(RP) :: cbh ! cloud base hight
2002  real(RP) :: rcbh ! E @ Zhang and Fritsch 1986 eq.13
2003  !
2004  real(RP) :: dens_d ! density of downdraft @ LFS
2005  real(RP) :: rhbar ! rerative humidity RH, bar is mean of
2006  real(RP) :: rhh ! decreace rhh
2007  !! humidity adjust ment tempolaly variables VV
2008  real(RP) :: dssdt
2009  real(RP) :: dtmp
2010  real(RP) :: qsrh
2011  real(RP) :: RL
2012  real(RP) :: T1rh
2013  !! ^^^^^^^^^^^^^
2014  real(RP) :: dpthtmp ! temporaly depth of downdraft source layer
2015  real(RP) :: dpthdet ! downdraft detrainment depth (Pa)
2016  real(RP) :: temp_d(ka) ! downdraft temperature
2017  real(RP) :: tempv_d(ka) ! downdraft virtual temperature
2018  real(RP) :: theta_ed(ka) ! downdraft equivalent theta
2019  real(RP) :: qvsd(ka) ! saturate watervapor in downdraft
2020  real(RP) :: qvs_tmp ! saturate qv temp
2021  real(RP) :: exn(ka) ! exner function
2022  real(RP) :: f_dmf ! factor of dmf Kain(2004) eq.11
2023  real(RP) :: dtempmlt ! melting effect of temperatture
2024  real(RP) :: es ! saturate vapor pressure
2025  real(RP) :: ddinc ! downdraft ??atode
2026  real(RP) :: prcpmlt ! melt precp internal
2027  integer :: k_dstart ! 150mb below of downdraft start (near USL top)
2028  integer :: k_lfstmp ! downdraft start layer
2029  integer :: k_ldt ! index top of downdrraft detraiment layer
2030  integer :: k_ldb ! index botom of downdraft detrainmengt layer
2031  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
2032  !! start cood
2033  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
2034 
2035  ! if no convection
2036  if (i_convflag == 2) return ! if 3d, may be change
2037  ! detremin
2038  do kk = ks, ke
2039  if (pres(kk) >= pres(ks)*0.5_rp) k_z5 = kk
2040  end do
2041  ! calc wind speed at LCL and cloud top and
2042  wspd(1) = sqrt(u(k_lcl)*u(k_lcl) + v(k_lcl)*v(k_lcl))
2043  wspd(2) = sqrt(u(k_z5)*u(k_z5) + v(k_z5)*v(k_z5))
2044  wspd(3) = sqrt(u(k_top)*u(k_top) + v(k_top)*v(k_top))
2045  ! calc wind shear and precipitation efficiency
2046  ! precipitation efficiency is calc two way
2047  ! 1. wind shear 2.cloud base hight
2048  ! we use mean of them
2049  if(wspd(3) > wspd(1)) then
2050  shsign = 1._rp
2051  else
2052  shsign = -1._rp
2053  end if
2054  vws = (u(k_top) - u(k_lcl))*(u(k_top) - u(k_lcl)) &
2055  + (v(k_top) - v(k_lcl))*(v(k_top) - v(k_lcl))
2056 
2057  vws = 1.e3_rp*shsign*sqrt(vws)/(z_kf(k_top) - z_kf(k_lcl))
2058  !! PEF(precipitation efficiency ) from Fritsch and Chappel 1980
2059  ! wind shear type
2060  pef_wind = 1.591_rp + vws*(-0.639_rp + vws*(9.53e-2_rp - vws*4.96e-3_rp) )
2061  ! 0.2 < pef_wind< 0.9
2062  pef_wind = max(pef_wind,0.2_rp)
2063  pef_wind = min(pef_wind,0.9_rp)
2064 
2065  !! 2. cloud base height
2066  ! refarence Zhang and Fritsch 1986
2067  cbh = (zlcl - z_kf(ks))*3.281e-3_rp ! convert m-> feet that's Amaerican
2068  if( cbh < 3._rp) then
2069  rcbh = 0.02_rp
2070  else
2071  rcbh = 0.96729352_rp + cbh*(-0.70034167_rp + cbh*(0.162179896_rp + cbh*(- &
2072  1.2569798e-2_rp + cbh*(4.2772e-4_rp - cbh*5.44e-6_rp))))
2073  end if
2074  if(cbh > 0.25) rcbh = 2.4_rp
2075  pef_cloud = 1._rp/(1._rp + rcbh)
2076  pef_cloud = min(pef_cloud,0.9_rp)
2077  ! pef is mean of wind shear type and cloud base heigt type
2078  pef = 0.5_rp*(pef_wind + pef_cloud)
2079  !!--------------------------------------------------
2080  !! compute downdraft propeties
2081  !!--------------------------------------------------
2082  tder = 0._rp ! initialize evapolation valuavle
2083  if(i_convflag == 1) then ! down devap
2084  k_lfs = ks ! shallow convection
2085  else
2086  !! start downdraft about 150 hPa above cloud base
2087  k_dstart = k_pbl + 1 ! index of usl layer top +1
2088  k_lfstmp = k_let - 1
2089  do kk = k_dstart+1, ke
2090  dpthtmp = pres(k_dstart) - pres(kk)
2091  if(dpthtmp > 150.e2_rp) then ! if dpth > 150hpa then
2092  k_lfstmp = kk
2093  exit
2094  end if
2095  end do
2096  k_lfstmp = min(k_lfstmp, k_let - 1) ! k_let is equivalent temperature layer index then
2097  k_lfs = k_lfstmp
2098 
2099  !... if LFS in not at least 50 mb above cloud base (implying that the
2100  !... level of equil temp, let ,is just above cloud base) do not allow downdraft
2101  ! dondraft layer
2102  ! downdraft is saturated
2103  if(pres(k_dstart) - pres(k_lfs) > 50.e2_rp) then ! LFS > 50mb(minimum of downdraft source layer)
2104  theta_ed(k_lfs) = theta_ee(k_lfs)
2105  qv_d(k_lfs) = qv(k_lfs)
2106  !! call tpmix2dd to find wet-bulb temperature and qv
2107  call tpmix2dd(pres(k_lfs),theta_ed(k_lfs),temp_d(k_lfs),qvs_tmp)
2108  call calcexn(pres(k_lfs),qvs_tmp,exn(k_lfs))
2109  !! exn(kk) = (PRE00/pres(k_lfs))**(0.2854*(1._RP - 0.28_RP*qv_d(k_lfs)))
2110  !! theta_d(k_lfs) = temp_d(k_lfs)*(PRE00/pres(k_lfs))**(0.2854*(1._RP -
2111  !! 0.28*qv_d(k_lfs)))
2112  theta_d(k_lfs) = temp_d(k_lfs)*exn(k_lfs)
2113  !! take a first guess at hte initial downdraft mass flux
2114  tempv_d(k_lfs) = temp_d(k_lfs)*(1._rp + 0.608_rp*qvs_tmp)
2115  dens_d = pres(k_lfs)/(r*tempv_d(k_lfs))
2116  dmf(k_lfs) = -(1._rp - pef)*1.e-2_rp*deltax*deltax*dens_d ! AU0 = 1.e-2*dx**2
2117  downent(k_lfs) = dmf(k_lfs)
2118  downdet(k_lfs) = 0._rp
2119  rhbar = rh(k_lfs)*deltap(k_lfs)
2120  dpthtmp = deltap(k_lfs)
2121  !! calc downdraft entrainment and (detrainment=0.)
2122  !! and dmf
2123  !! downdraft theta and q
2124  do kk = k_lfs-1,k_dstart,-1
2125  kkp1 = kk + 1
2126  downent(kk) = downent(k_lfs)*ems(kk)/ems(k_lfs)
2127  downdet(kk) = 0._rp
2128  dmf(kk) = dmf(kkp1) + downent(kk)
2129  theta_ed(kk) = ( theta_ed(kkp1)*dmf(kkp1) + theta_ee(kk)*downent(kk) )/dmf(kk)
2130  qv_d(kk) = ( qv_d(kkp1)*dmf(kkp1) + qv(kk)*downent(kk) )/dmf(kk)
2131  dpthtmp = dpthtmp + deltap(kk)
2132  rhbar = rhbar + rh(kk)*deltap(kk) ! rh average@ downdraft layer
2133  end do
2134  rhbar = rhbar/dpthtmp
2135  f_dmf = 2._rp*(1._rp - rhbar) ! Kain 2004 eq.11
2136  !
2137  !! calc melting effect
2138  !! first, cocalc total frozen precipitation generated..
2139  prcpmlt = 0._rp
2140  do kk = k_lcl,k_top
2141  prcpmlt = prcpmlt + flux_qs(kk)
2142  end do
2143  if (k_lc < k_ml ) then ! if below melt level layer then
2144  !! RLF is EMELT
2145  !! dtempmlt = RLF*prcpmlt/(CP*umf(k_lcl))
2146  dtempmlt = emelt*prcpmlt/(cp*umf(k_lcl))
2147  else
2148  dtempmlt = 0._rp
2149  end if
2150  call tpmix2dd(pres(k_dstart),theta_ed(k_dstart),temp_d(k_dstart),qvs_tmp)
2151  temp_d(k_dstart) = temp_d(k_dstart) - dtempmlt
2152  !! use check theis subroutine is this
2153  call atmos_saturation_psat_liq(es,temp_d(k_dstart)) !saturation vapar pressure
2154  qvs_tmp = 0.622_rp*es/(pres(k_dstart) - es )
2155  !! Bolton 1980 pseudoequivalent potential temperature
2156  theta_ed(k_dstart) = temp_d(k_dstart)*(pre00/pres(k_dstart))**(0.2854_rp*(1._rp - 0.28_rp*qvs_tmp))* &
2157  & exp((3374.6525_rp/temp_d(k_dstart)-2.5403_rp)*qvs_tmp*(1._rp + 0.81_rp*qvs_tmp))
2158  k_ldt = min(k_lfs-1, k_dstart-1)
2159  dpthdet = 0._rp
2160  do kk = k_ldt,ks,-1 ! start calc downdraft detrain layer index
2161  !!
2162  dpthdet = dpthdet + deltap(kk)
2163  theta_ed(kk) = theta_ed(k_dstart)
2164  qv_d(kk) = qv_d(k_dstart)
2165  call tpmix2dd(pres(kk),theta_ed(kk),temp_d(kk),qvs_tmp)
2166  qvsd(kk) = qvs_tmp
2167  !... specify RH decrease of 20%/km indowndraft
2168  rhh = 1._rp - 2.e-4_rp*(z_kf(k_dstart) -z_kf(kk) ) ! 0.2/1000.
2169  !!
2170  !!... adjust downdraft temp,qv to secified RH :
2171  !!
2172  !! Kain(2004)
2173  if(rhh < 1._rp) then
2174  !!
2175  dssdt = (cliq - bliq*dliq)/((temp_d(kk) - dliq)*(temp_d(kk) - dliq))
2176  rl = xlv0 - xlv1*temp_d(kk)
2177  dtmp = rl*qvs_tmp*(1._rp - rhh )/(cp + rl*rhh*qvs_tmp*dssdt )
2178  t1rh = temp_d(kk) + dtmp
2179  call atmos_saturation_psat_liq(es,t1rh) !saturation vapar pressure
2180  es = rhh*es
2181  qsrh = 0.622_rp*es/(pres(kk) - es)
2182  if(qsrh < qv_d(kk) ) then
2183  qsrh = qv_d(kk)
2184  t1rh = temp_d(kk) + (qvs_tmp - qsrh)*rl/cp
2185  end if
2186 
2187  temp_d(kk) = t1rh
2188  qvs_tmp = qsrh
2189  qvsd(kk) = qvs_tmp
2190  !!
2191  end if
2192  !!
2193  tempv_d(kk) = temp_d(kk)*( 1._rp + 0.608_rp*qvsd(kk) )
2194  if(tempv_d(kk) > tempv(kk) .or. kk == ks) then
2195  k_ldb = kk
2196  exit
2197  end if
2198  !!
2199  end do ! end calc downdraft detrain layer index
2200  if((pres(k_ldb)-pres(k_lfs)) > 50.e2_rp ) then ! minimum downdraft depth !
2201  do kk = k_ldt,k_ldb,-1
2202  kkp1 = kk + 1
2203  downdet(kk) = -dmf(k_dstart)*deltap(kk)/dpthdet
2204  downent(kk) = 0._rp
2205  dmf(kk) = dmf(kkp1) + downdet(kk)
2206  tder = tder + (qvsd(kk) - qv_d(kk))*downdet(kk)
2207  qv_d(kk) = qvsd(kk)
2208  call calcexn(pres(kk),qv_d(kk),exn(kk))
2209  theta_d(kk) = temp_d(kk)*exn(kk)
2210  end do
2211  end if
2212  end if ! LFS>50mb
2213  end if ! down devap
2214  !!
2215  !!... if downdraft does not evaporate any water for specified relative
2216  !!...humidity, no downdraft is allowed ..
2217  !!and shallow convection does not have downdraft
2218  if (tder < 1._rp) then ! dmf modify
2219  prcp_flux = totalprcp
2220  cpr = totalprcp
2221  tder = 0._rp
2222  k_ldb = k_lfs
2223  do kk = ks, k_top
2224  dmf(kk) = 0._rp
2225  downdet(kk) = 0._rp
2226  downent(kk) = 0._rp
2227  theta_d(kk) = 0._rp
2228  temp_d(kk) = 0._rp
2229  qv_d(kk) = 0._rp
2230  end do
2231  else
2232  ddinc = -f_dmf*umf(k_lcl)/dmf(k_dstart) ! downdraft keisuu Kain 2004 eq.11
2233  if(tder*ddinc > totalprcp) then
2234  ddinc = totalprcp/tder ! rate of prcp/evap
2235  end if
2236  tder = tder*ddinc
2237  do kk = k_ldb,k_lfs
2238  dmf(kk) = dmf(kk)*ddinc
2239  downent(kk) = downent(kk)*ddinc
2240  downdet(kk) = downdet(kk)*ddinc
2241  end do
2242  cpr = totalprcp
2243  prcp_flux = totalprcp - tder ! precipitation - evapolate water
2244  pef = prcp_flux/totalprcp
2245  !
2246  !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
2247  ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
2248  ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
2249  !
2250  ! DO NK=LC,LFS
2251  ! UMF(NK)=UMF(NK)*UPDINC
2252  ! UDR(NK)=UDR(NK)*UPDINC
2253  ! UER(NK)=UER(NK)*UPDINC
2254  ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
2255  ! PPTICE(NK)=PPTICE(NK)*UPDINC
2256  ! DETLQ(NK)=DETLQ(NK)*UPDINC
2257  ! DETIC(NK)=DETIC(NK)*UPDINC
2258  ! ENDDO
2259  !
2260  !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
2261  !...DOWNDRAFT...
2262  !
2263  if (k_ldb > ks) then
2264  do kk = ks,k_ldb-1
2265  dmf(kk) = 0._rp
2266  downdet(kk) = 0._rp
2267  downent(kk) = 0._rp
2268  theta_d(kk) = 0._rp
2269  temp_d(kk) = 0._rp
2270  qv_d(kk) = 0._rp
2271  end do
2272  end if
2273  ! no dmf is above LFS layer
2274  do kk = k_lfs+1,ke
2275  dmf(kk) = 0._rp
2276  downdet(kk) = 0._rp
2277  downent(kk) = 0._rp
2278  theta_d(kk) = 0._rp
2279  temp_d(kk) = 0._rp
2280  qv_d(kk) = 0._rp
2281  end do
2282  !! no temperature and qv avave downdraft detrainment startlayer (k_ldt)
2283  do kk = k_ldt+1,k_lfs-1
2284  temp_d(kk) = 0._rp
2285  qv_d(kk) = 0._rp
2286  theta_d(kk) = 0._rp
2287  end do
2288  end if ! dmf modify
2289 
2290  return
2291  end subroutine cp_kf_downdraft
2292  !!------------------------------------------------------------------------------
2293  !! compute properties for compensational subsidence
2294  !! return new valuavles
2295  !!------------------------------------------------------------------------------
2296  !!compute properties for compensational subsidence
2297  !!and determin dt variables
2298  subroutine cp_kf_compensational (&
2299  !![in]
2300  dz_kf,z_kf ,&
2301  pres ,& ! pressure
2302  deltap ,& ! deltap
2303  temp_bf ,& ! temperature
2304  qv ,& ! water vapor
2305  ems ,&
2306  emsd ,&
2307  !! trigger
2308  presmix ,& ! usl layer pressure depth
2309  zmix ,& ! usl layer pressure depth
2310  dpthmx ,& ! usl layer depth
2311  cape ,& !cape
2312  temp_u ,& ! updraft temperature
2313  qvdet ,& ! updraft water vapor
2314  umflcl ,& ! umf @LCL
2315  umf ,& ! UMF
2316  upent ,& ! updraft entrainment
2317  updet ,& ! updraft detrainment
2318  qcdet ,& ! updraft detrainment qc
2319  qidet ,& ! updraft detrainment qi
2320  qc ,& ! cloud water
2321  qi ,& ! cloud ice
2322  flux_qr ,& ! rain flux
2323  flux_qs ,& ! snow flux
2324  umfnewdold ,& ! 1/(umfnew/umfold ratio) !
2325  !!downdraft
2326  wspd ,& ! wind speed 1. 2 is used
2327  dmf ,& ! DMF
2328  downent ,& ! downdraft entrainment
2329  downdet ,& ! downdraft detrainment
2330  qv_d ,& ! downdraft detrainment qv
2331  theta_d ,& ! potential temperature @downdraft
2332  prcp_flux ,& ! precipitation
2333  tder ,& ! evapolation
2334  cpr ,& ! all precipitation before consider cloud bottom evaporation
2335  !
2336  i_convflag ,& ! intent inout convective flag
2337  !triger
2338  k_top ,&
2339  k_lcl_bf ,&
2340  k_lc ,&
2341  k_pbl ,&
2342  k_ml ,&
2343  !downdraft
2344  k_lfs ,&
2345  !![OUT] new valuavles after timestep
2346  temp_g ,&
2347  qv_g ,&
2348  qc_nw ,&
2349  qi_nw ,&
2350  qr_nw ,&
2351  qs_nw ,&
2352  rainrate_cp ,&
2353  nic ,&
2354  cldfrac_kf ,&
2355  timecp ,&
2356  time_advec)
2357  use scale_precision
2358  use scale_grid_index
2359  use scale_const,only :&
2360  cp => const_cpdry , &
2361  pre00 => const_pre00 , &
2362  grav => const_grav , &
2363  emelt => const_emelt
2364  use scale_atmos_saturation ,only :&
2365  atmos_saturation_psat_liq
2366  use scale_time , only :&
2367  dt => time_dtsec_atmos_phy_cp
2368  use scale_process, only: &
2369  prc_mpistop
2370  implicit none
2371  ! in
2372  ! form init
2373  real(RP),intent(in) :: dz_kf(ka),z_kf(ka) ! delta Z, and haight [m] full point
2374  real(RP),intent(in) :: pres(ka) ! pressure [Pa]
2375  real(RP),intent(in) :: deltap(ka) ! delta pressure
2376  real(RP),intent(in) :: temp_bf(ka) ! temperature berore
2377  real(RP),intent(in) :: qv(ka) ! cloud vaper mixing ratio
2378  real(RP),intent(in) :: ems(ka)
2379  real(RP),intent(in) :: emsd(ka)
2380  ! from trigger
2381  real(RP),intent(in) :: presmix
2382  real(RP),intent(in) :: zmix
2383  real(RP),intent(in) :: dpthmx ! usl layer depth
2384  real(RP),intent(in) :: cape
2385  real(RP),intent(in) :: temp_u(ka)
2386  real(RP),intent(in) :: qvdet(ka)
2387  real(RP),intent(in) :: umflcl ! umf@LCL internal work(umf(k_lcl - 1)
2388  real(RP),intent(inout) :: umf(ka) ! UMF
2389  real(RP),intent(inout) :: upent(ka) ! updraft entrainment
2390  real(RP),intent(inout) :: updet(ka) ! updraft detrainment
2391  real(RP),intent(inout) :: qcdet(ka) ! updraft detrainment qc
2392  real(RP),intent(inout) :: qidet(ka) ! updraft detrainment qi
2393  real(RP),intent(in) :: qc(ka) ! cloud water mixing ratio
2394  real(RP),intent(in) :: qi(ka) ! cloud ice mixing ratio
2395  real(RP),intent(in) :: flux_qr(ka)
2396  real(RP),intent(in) :: flux_qs(ka)
2397  real(RP),intent(in) :: umfnewdold(ka) ! umfnew/umfold ratio
2398  ! from dwndraft
2399  real(RP),intent(in) :: wspd(3) ! wind speed 1. 2 is used
2400  real(RP),intent(inout) :: dmf(ka) ! DMF
2401  real(RP),intent(inout) :: downent(ka) ! downdraft entrainment
2402  real(RP),intent(inout) :: downdet(ka) ! downdraft detrainment
2403  real(RP),intent(in) :: qv_d(ka) ! downdraft detrainment qv
2404  real(RP),intent(in) :: theta_d(ka)
2405  real(RP),intent(inout) :: prcp_flux
2406  real(RP),intent(inout) :: tder ! evapolation
2407  real(RP),intent(in) :: cpr
2408  integer,intent(inout) :: I_convflag ! intent inout
2409  ! from trigger
2410  integer,intent(in) :: k_top
2411  integer,intent(inout) :: k_lcl_bf
2412  integer,intent(in) :: k_lc
2413  integer,intent(in) :: k_pbl
2414  integer,intent(in) :: k_ml
2415  integer,intent(in) :: k_lfs
2416  ! out
2417  real(RP),intent(out) :: temp_g(ka) ! temporarly temperature -> after iteration then new variable
2418  real(RP),intent(out) :: qv_g(ka)
2419  real(RP),intent(out) :: qc_nw(ka)
2420  real(RP),intent(out) :: qi_nw(ka)
2421  real(RP),intent(out) :: qr_nw(ka)
2422  real(RP),intent(out) :: qs_nw(ka)
2423  integer,intent (out) :: nic ! num of step convection active [step]
2424  ! dt(cumulus parametrization deltat)/cumulus parameterization time scale
2425  real(RP),intent(out) :: rainrate_cp !PPTFLX(prcp_flux)/deltax**2/dt convective rain rate
2426  real(RP),intent(out) :: cldfrac_KF(ka,2) ! 1.shallow and 2.deep
2427  real(RP),intent(out) :: timecp ! timescale of cumulus parameterization
2428  real(RP),intent(out) :: time_advec ! advection timescale
2429  ![INTERNAL WORK]
2430  ! iteratevar
2431  real(RP) :: umf2(ka) ! UMF
2432  real(RP) :: upent2(ka) ! updraft entrainment
2433  real(RP) :: updet2(ka) ! updraft detrainment
2434  real(RP) :: qcdet2(ka) ! updraft detrainment qc
2435  real(RP) :: qidet2(ka) ! updraft detrainment qi
2436  real(RP) :: dmf2(ka) ! DMF
2437  real(RP) :: downent2(ka) ! downdraft entrainment
2438  real(RP) :: downdet2(ka) ! downdraft detrainment
2439  real(RP) :: prcp_flux2 ! precpitation flux
2440  real(RP) :: tder2 ! evaporation
2441  !
2442  real(RP) :: tkemax !tkemax tuning prameter
2443  real(RP) :: z_lcl ! lcl layer hight
2444  real(RP) :: theta(ka) ! theta is not same as SCALE theta. This theta is only assume qv
2445  real(RP) :: theta_u(ka) ! theta in updraft
2446  real(RP) :: theta_eu(ka) ! equivalent PT updraft
2447  real(RP) :: theta_eg(ka) ! equivalent PT environment
2448  real(RP) :: exn(ka) ! exner function
2449  real(RP) :: qv_env ! environment qv
2450  real(RP) :: qv_mix ! USL layer mean qv
2451  real(RP) :: qv_gu(ka) ! updraft qv (used calc CAPE)
2452  real(RP) :: temp_env ! temporarly environment temperature lcl layer
2453  real(RP) :: tempv_env ! temporarly environment virtual temperature lcl layer
2454  real(RP) :: temp_lcl ! LCL temperature used calcurate CAPE
2455  real(RP) :: tempv_lcl ! temporarly environment virtual temperature lcl layer
2456  real(RP) :: temp_mix ! temporarly environment temperature USL layer mean
2457  real(RP) :: temp_gu(ka) ! temporarly updraft temperature
2458  real(RP) :: tempv_g(ka) ! temporaly virtual
2459  real(RP) :: tempvq_u(ka) ! temporaly virtual
2460  real(RP) :: es ! saturate vapor pressure
2461  real(RP) :: qvss ! saturate vapor pressure mixingratio
2462  ! calc dew point temperature etc.
2463  real(RP) :: DQ,TDPT,DSSDT,emix,RL,TLOG
2464  !
2465  real(RP) :: vconv
2466  real(RP) :: dzz ! lcl layer depth used calcurate CAPE for interpolate lcl layer
2467  real(RP) :: deltaz ! lcl layer depth used calcurate CAPE for interpolate lcl layer
2468  real(RP) :: dilbe ! used calculate CAPE
2469  integer :: k_lcl,k_lclm1 ! LCL and LCL-1 layer index used calucrat CAPE
2470  ! inter porlatevariable
2471  ! new variables
2472  real(RP) :: theta_nw(ka) ! new PT (used iteration
2473  real(RP) :: theta_g(ka) ! new PT
2474  real(RP) :: qv_nw(ka) ! tempolaly new qv, becouse itaration
2475  real(RP) :: dpth ! pressure depth of cloud used check
2476  real(RP) :: cape_g ! new cape clculate after timestep
2477  real(RP) :: dcape ! deltacape compared 10% of original
2478  real(RP) :: fxm(ka) ! flux factor
2479  real(RP) :: f_cape ,f_capeold ! cape ratio new/old
2480  real(RP) :: stab ! 0.95 stablevariable
2481  real(RP) :: dtt_tmp
2482  real(RP) :: dtt ! deltat of cumulus prameterization determined by layer depth/omega
2483  real(RP) :: deltat ! deltat of cumulus parameterization same as dtt but reduced
2484  integer :: istop ! effor flag 1 mass blance check
2485  integer :: ncount ! countor for iteration
2486  integer :: ntimecount ! timecount do loop index
2487  integer :: nstep ! max step of ntimecount
2488  integer :: noiter ! iteration flag
2489  ! temporaly
2490  integer :: kk,kkp1 ! index of do loop
2491  integer :: k_lmax ! max of k_lcl,k_lfs
2492  real(RP) :: tma,tmb,tmm ! eff
2493  real(RP) :: acoeff,bcoeff ! eff
2494  real(RP) :: evac ! shallow convection TKE factor
2495  real(RP) :: ainc,ainctmp, aincmx,aincold ! factors ainctmp is tmpvariable; aincmx is max of ainc
2496  real(RP) :: omg(ka) ! pressure velocity
2497  real(RP) :: topomg ! cloud top omg calc by updraft
2498  real(RP) :: fbfrc ! precpitation to be fedback 0.0 -1.0 shallo-> 1.0(no rain) deep->0.0
2499  real(RP) :: frc2
2500  real(RP) :: dfda
2501  real(RP) :: domg_dp(ka) ! d omega/dp
2502  real(RP) :: absomgtc,absomg
2503  real(RP) :: f_dp
2504  ! fluxs
2505  real(RP) :: theta_fxin(ka) ! cmpensational subsidence flux form
2506  real(RP) :: theta_fxout(ka) ! cmpensational subsidence flux form
2507  real(RP) :: qv_fxin(ka) ! cmpensational subsidence flux form
2508  real(RP) :: qv_fxout(ka) ! cmpensational subsidence flux form
2509  real(RP) :: qc_fxin(ka) ! cmpensational subsidence flux form
2510  real(RP) :: qc_fxout(ka) ! cmpensational subsidence flux form
2511  real(RP) :: qi_fxin(ka) ! cmpensational subsidence flux form
2512  real(RP) :: qi_fxout(ka) ! cmpensational subsidence flux form
2513  real(RP) :: qr_fxin(ka) ! cmpensational subsidence flux form
2514  real(RP) :: qr_fxout(ka) ! cmpensational subsidence flux form
2515  real(RP) :: qs_fxin(ka) ! cmpensational subsidence flux form
2516  real(RP) :: qs_fxout(ka) ! cmpensational subsidence flux form
2517  real(RP) :: rainfb(ka) ! rain fall
2518  real(RP) :: snowfb(ka) ! snow fall
2519  ! evariable moisture budeget
2520  real(RP) :: err ! error (tmp var)
2521  real(RP) :: qinit ! init water condensation (only qv)
2522  real(RP) :: qfinl ! fineal water condensation (producted by cumulus parameterization)
2523  ! warm rain
2524  real(RP) :: cpm ! tempolaly variable
2525  real(RP) :: UMF_tmp ! tempolaly variable
2526  real(RP) :: xcldfrac ! tempolaly variable
2527  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
2528  !! start cood
2529  !!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
2530 
2531  !! initialize
2532  if(i_convflag == 2) return ! noconvection
2533  ! start cood
2534  k_lcl = k_lcl_bf
2535  fbfrc = 0._rp
2536  if(i_convflag == 1) fbfrc = 1._rp
2537  !! set timecp(TIMEC)
2538  ! time scale of adjustment
2539  vconv = 0.5_rp*(wspd(1) + wspd(2)) ! k_lcl + k_z5
2540  timecp = deltax/vconv
2541  time_advec = timecp ! advection time sclale
2542  ! 30 minutes < timecp < 60 minutes
2543  timecp = max(deeplifetime, timecp)
2544  timecp = min(3600._rp, timecp)
2545  if(i_convflag == 1) timecp = shallowlifetime ! shallow convection timescale is 40 minutes
2546  nic = nint(timecp/dt)
2547  ! timecp is KF timescale
2548  ! must be bigger than "cumplus parameterization timestep " given by namelist
2549  timecp = real( nic,rp )*dt! determin timecp not change below
2550  !
2551  ! maximam of ainc calculate
2552  aincmx = 1000._rp
2553  k_lmax = max(k_lcl,k_lfs)
2554  do kk = k_lc,k_lmax
2555  if((upent(kk)-downent(kk)) > 1.e-3_rp) then
2556  ainctmp = ems(kk)/( (upent(kk) -downent(kk))*timecp )
2557  aincmx = min(aincmx,ainctmp)
2558  end if
2559  end do
2560  ainc = 1._rp
2561  if(aincmx < ainc ) ainc = aincmx
2562  !! for interpolation save original variable
2563  tder2 = tder
2564  prcp_flux2 = prcp_flux
2565  do kk = ks,k_top
2566  umf2(kk) = umf(kk)
2567  upent2(kk) = upent(kk)
2568  updet2(kk) = updet(kk)
2569  qcdet2(kk) = qcdet(kk)
2570  qidet2(kk) = qidet(kk)
2571  dmf2(kk) = dmf(kk)
2572  downent2(kk) = downent(kk)
2573  downdet2(kk) = downdet(kk)
2574  end do
2575  f_cape = 1._rp
2576  stab = 1.05_rp - delcape ! default 0.95
2577  noiter = 0 ! noiter=1 then stop iteration
2578  istop = 0
2579  if (i_convflag == 1) then ! shallow convection
2580  ! refarence Kain 2004
2581  tkemax = 5._rp
2582  evac = 0.50_rp*tkemax*1.e-1_rp
2583  ainc = evac*dpthmx*deltax*deltax/( umflcl*grav*timecp)
2584  tder = tder2*ainc
2585  prcp_flux = prcp_flux2*ainc
2586  do kk = ks,k_top
2587  umf(kk) = umf2(kk)*ainc
2588  upent(kk) = upent2(kk)*ainc
2589  updet(kk) = updet2(kk)*ainc
2590  qcdet(kk) = qcdet2(kk)*ainc
2591  qidet(kk) = qidet2(kk)*ainc
2592  dmf(kk) = dmf2(kk)*ainc
2593  downent(kk) = downent2(kk)*ainc
2594  downdet(kk) = downdet2(kk)*ainc
2595  end do
2596  end if
2597  !! theta set up by Emanuel Atomospheric convection, 1994 111p
2598  !! original KF theta is calced apploximatly.
2599  do kk = ks,k_top
2600  call calcexn(pres(kk),qv(kk),exn(kk))
2601  theta(kk) = temp_bf(kk)*exn(kk)
2602  call calcexn(pres(kk),qvdet(kk),exn(kk))
2603  theta_u(kk) = temp_u(kk)*exn(kk)
2604  end do
2605  temp_g(k_top+1:ke) = temp_bf(k_top+1:ke)
2606  qv_g(k_top+1:ke) = qv(k_top+1:ke)
2607  omg(:) = 0._rp ! initialize
2608  !!
2609  !!@@@ start main loop @@@
2610  do ncount = 1,10 ! itaration
2611  dtt = timecp
2612  do kk = ks,k_top
2613  domg_dp(kk) = -(upent(kk) - downent(kk) - updet(kk) - downdet(kk))*emsd(kk)
2614  if(kk > ks)then
2615  omg(kk) = omg(kk-1) - deltap(kk-1)*domg_dp(kk-1)
2616  absomg = abs(omg(kk))
2617  absomgtc = absomg*timecp
2618  f_dp = 0.75*deltap(kk-1)
2619  if(absomgtc > f_dp)THEN
2620  dtt_tmp = f_dp/abs(omg(kk))
2621  dtt=min(dtt,dtt_tmp) ! it is use determin nstep
2622  end if
2623  end if
2624  end do
2625  !! theta_nw and qv_nw has valus only in 1 to k_top
2626  do kk = ks, k_top
2627  theta_nw(kk) = theta(kk)
2628  qv_nw(kk) = qv(kk)
2629  fxm(kk) = omg(kk)*deltax*deltax/grav ! fluxmass
2630  end do
2631  nstep = nint(timecp/dtt + 1) ! how many time step forwad
2632  deltat = timecp/real(nstep,RP) ! deltat*nstep = timecp
2633  do ntimecount = 1, nstep
2634  ! tempolaly time forward start iteration
2635  !!... assign theta and q variables at the top and bottom of each layer based
2636  !!... on sign of omega
2637  do kk = ks, k_top ! initialize
2638  theta_fxin(kk) = 0._rp
2639  theta_fxout(kk) = 0._rp
2640  qv_fxin(kk) = 0._rp
2641  qv_fxout(kk) = 0._rp
2642  end do
2643  do kk = ks+1,k_top ! calc flux variable
2644  if(omg(kk) <= 0._rp) then
2645  theta_fxin(kk) = -fxm(kk)*theta_nw(kk-1)
2646  qv_fxin(kk) = -fxm(kk)*qv_nw(kk-1)
2647  theta_fxout(kk - 1) = theta_fxout(kk-1) + theta_fxin(kk)
2648  qv_fxout(kk - 1) = qv_fxout(kk-1) + qv_fxin(kk)
2649  else
2650  theta_fxout(kk) = fxm(kk)*theta_nw(kk)
2651  qv_fxout(kk) = fxm(kk)*qv_nw(kk)
2652  theta_fxin(kk - 1) = theta_fxin(kk-1) + theta_fxout(kk)
2653  qv_fxin(kk - 1) = qv_fxin(kk-1) + qv_fxout(kk)
2654  end if
2655  end do
2656  !! update the theta and qv variables at each level
2657  !! only theta and qv calc cape below and if cape > 10%of old cape then iterate
2658  do kk = ks, k_top
2659  theta_nw(kk) = theta_nw(kk) &
2660  + (theta_fxin(kk) - theta_fxout(kk) &
2661  + updet(kk)*theta_u(kk) + downdet(kk)*theta_d(kk) &
2662  - ( upent(kk) - downent(kk) )*theta(kk) ) *deltat*emsd(kk)
2663  qv_nw(kk) = qv_nw(kk) &
2664  + (qv_fxin(kk) - qv_fxout(kk) &
2665  + updet(kk)*qvdet(kk) + downdet(kk)*qv_d(kk) &
2666  - ( upent(kk) - downent(kk) )*qv(kk) )*deltat*emsd(kk)
2667  end do
2668  end do ! ntimecount
2669 
2670  do kk = ks, k_top
2671  theta_g(kk) = theta_nw(kk)
2672  qv_g(kk) = qv_nw(kk)
2673  end do
2674  ! Check to see if mixing ratio dips below zero anywere ; if so, borrow
2675  ! moisture forom adjacent layers to bring it back up abobe zero
2676  do kk = ks,k_top
2677  if(qv_g(kk) < 0._rp ) then ! negative moisture
2678  if(kk == ks) then
2679  write(*,*) "error qv<0 @ Kain-Fritsch cumulus parameterization"
2680  write(*,*) "@sub scale_atmos_phy_cp_kf",__file__,__line__
2681  call prc_mpistop
2682  end if
2683  kkp1 = kk + 1
2684  if(kk == k_top) then
2685  kkp1 = k_lcl
2686  end if
2687  tma = qv_g(kkp1)*ems(kkp1)
2688  tmb = qv_g(kk-1)*ems(kk-1)
2689  tmm = (qv_g(kk) - 1.e-9_rp )*ems(kk)
2690  bcoeff = -tmm/((tma*tma)/tmb + tmb)
2691  acoeff = bcoeff*tma/tmb
2692  tmb = tmb*(1._rp - bcoeff)
2693  tma = tma*(1._rp - acoeff)
2694  qv_g(kk) = 1.e-9_rp
2695  qv_g(kkp1) = tma*emsd(kkp1)
2696  qv_g(kk-1) = tmb*emsd(kk-1)
2697  end if
2698  end do
2699  ! clculate top layer omega and conpare determ omg
2700  topomg = (updet(k_top) - upent(k_top))*deltap(k_top)*emsd(k_top)
2701  if( abs(topomg - omg(k_top)) > 1.e-3_rp) then ! not same omega velocity error
2702  istop = 1
2703  write(*,*) "xxxERROR@KF omega is not consistent",ncount
2704  write(*,*) "omega error",abs(topomg - omg(k_top)),k_top,topomg,omg(k_top)
2705  call prc_mpistop
2706  end if
2707  ! convert theta to T
2708  do kk = ks,k_top
2709  call calcexn(pres(kk),qv_g(kk),exn(kk))
2710  temp_g(kk) = theta_g(kk)/exn(kk)
2711  tempv_g(kk) = temp_g(kk)*(1._rp + 0.608_rp*qv_g(kk))
2712  end do
2713  !------------------------------------------------------------------------------
2714  ! compute new cloud and change in available bouyant energy(CAPE)
2715  !------------------------------------------------------------------------------
2716  if(i_convflag == 1) then
2717  exit ! if shallow convection , calc Cape is not used
2718  end if
2719  temp_mix = 0._rp
2720  qv_mix = 0._rp
2721  do kk = k_lc, k_pbl
2722  temp_mix = temp_mix + deltap(kk)*temp_g(kk)
2723  qv_mix = qv_mix + deltap(kk)*qv_g(kk)
2724  ! presmix = presmix + deltap(kk)
2725  end do
2726  temp_mix = temp_mix/dpthmx
2727  qv_mix = qv_mix/dpthmx
2728  ! calc saturate water vapor pressure
2729  call atmos_saturation_psat_liq(es,temp_mix)
2730  qvss = 0.622_rp*es/(presmix -es) ! saturate watervapor
2731  !!
2732  !!... Remove supersaturation for diagnostic purposes, if necessary..
2733  !!
2734  if (qv_mix > qvss) then ! saturate then
2735  rl = xlv0 -xlv1*temp_mix
2736  cpm = cp*(1._rp + 0.887_rp*qv_mix)
2737  dssdt = qvss*(cliq-bliq*dliq)/( (temp_mix-dliq)**2)
2738  dq = (qv_mix -qvss)/(1._rp + rl*dssdt/cpm)
2739  temp_mix = temp_mix + rl/cp*dq
2740  qv_mix = qv_mix - dq
2741  temp_lcl = temp_mix
2742  else !same as detern trigger
2743  qv_mix = max(qv_mix,0._rp)
2744  emix = qv_mix*presmix/(0.6222_rp + qv_mix)
2745  tlog = log(emix/aliq)
2746  ! dew point temperature Bolton(1980)
2747  tdpt = (cliq - dliq*tlog)/(bliq - tlog)
2748  temp_lcl = tdpt - (0.212_rp + 1.57e-3_rp*(tdpt - tem00) - 4.36e-4_rp*(temp_mix - tem00))*(temp_mix - tdpt)
2749  temp_lcl = min(temp_lcl,temp_mix)
2750  end if
2751  tempv_lcl = temp_lcl*(1._rp + 0.608_rp*qv_mix)
2752  z_lcl = zmix + (temp_lcl - temp_mix)/gdcp
2753  do kk = k_lc, ke
2754  k_lcl = kk
2755  if( z_lcl <= z_kf(kk) ) exit
2756  end do
2757  ! estimate environmental tempeature and mixing ratio
2758  ! interpolate environment temperature and vapor at LCL
2759  k_lclm1 = k_lcl - 1
2760  deltaz = ( z_lcl - z_kf(k_lclm1) )/( z_kf(k_lcl)- z_kf(k_lclm1 ) )
2761  temp_env = temp_g(k_lclm1) + ( temp_g(k_lcl) - temp_g(k_lclm1) )*deltaz
2762  qv_env = qv_g(k_lclm1) + ( qv_g(k_lcl) - qv_g(k_lclm1) )*deltaz
2763  tempv_env = temp_env*( 1._rp + 0.608_rp*qv_env )
2764  !! pres_lcl=pres(k_lcl-1)+(pres(k_lcl-1)-pres(k_lcl-1))*deltaz
2765  theta_eu(k_lclm1)=temp_mix*(pre00/presmix)**(0.2854_rp*(1._rp - 0.28_rp*qv_mix))* &
2766  exp((3374.6525_rp/temp_lcl-2.5403_rp)*qv_mix*(1._rp + 0.81_rp*qv_mix))
2767  !!
2768  !!...COMPUTE ADJUSTED ABE(ABEG).(CAPE)
2769  !!
2770  cape_g = 0._rp ! cape "_g" add because
2771  do kk=k_lclm1,k_top-1 ! LTOPM1
2772  kkp1=kk+1
2773  theta_eu(kkp1) = theta_eu(kk)
2774  call tpmix2dd(pres(kkp1),theta_eu(kkp1),temp_gu(kkp1),qv_gu(kkp1)) ! get temp_gu and qv_gu
2775  tempvq_u(kkp1) = temp_gu(kkp1)*(1._rp + 0.608_rp*qv_gu(kkp1) - qc(kkp1)- qi(kkp1))
2776  if(kk == k_lclm1) then ! interporate
2777  dzz = z_kf(k_lcl) - z_lcl
2778  dilbe = ((tempv_lcl + tempvq_u(kkp1))/(tempv_env + tempv_g(kkp1)) - 1._rp)*dzz
2779  else
2780  dzz = dz_kf(kk)
2781  dilbe = ((tempvq_u(kk) + tempvq_u(kkp1))/(tempv_g(kk) + tempv_g(kkp1)) - 1._rp)*dzz
2782  end if
2783  if(dilbe > 0._rp) cape_g = cape_g + dilbe*grav
2784  !!
2785  !!...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
2786  !!
2787  call envirtht(pres(kkp1),temp_g(kkp1),qv_g(kkp1),theta_eg(kkp1)) ! calc get theta_eg
2788  !! theta_eg(environment theta_E )
2789  ! theta_eu(kkp1) = theta_eu(kkp1)*(1._RP/umfnewdold(kkp1)) + theta_eg(kkp1)*(1._RP - (1._RP/umfnewdold(kkp1)))
2790  theta_eu(kkp1) = theta_eu(kkp1)*(umfnewdold(kkp1)) + theta_eg(kkp1)*(1._rp - (umfnewdold(kkp1)))
2791  end do
2792  if (noiter == 1) exit ! noiteration
2793  dcape = max(cape - cape_g,cape*0.1_rp) ! delta cape
2794  f_cape = cape_g/cape ! ratio of cape new/old
2795  if(f_cape > 1._rp .and. i_convflag == 0) then ! if deep convection and cape is inclease this loop
2796  i_convflag = 2
2797  return
2798  end if
2799  if(ncount /= 1) then
2800  if(abs(ainc - aincold) < 1.e-4_rp) then ! IN cycle not change facter then exit iter next loop
2801  noiter = 1 ! exit this loop in nex step
2802  ainc = aincold
2803  cycle ! iter
2804  end if
2805  dfda = (f_cape - f_capeold)/(ainc - aincold)
2806  if (dfda > 0._rp ) then
2807  noiter = 1 ! exit this loop @next loop step
2808  ainc = aincold
2809  cycle ! iter
2810  end if
2811  end if
2812  aincold = ainc
2813  f_capeold = f_cape
2814  !! loop exit
2815  !! aincmx is upper limit of massflux factor
2816  !! if massflux factor 'ainc' is near "aincmax" then exit
2817  !! but need CAPE is less than 90% of original
2818  if (ainc/aincmx > 0.999 .and. f_cape > 1.05-stab ) then
2819  exit
2820  end if
2821  !! loop exit 1. or 2.
2822  !! 1. NEW cape is less than 10% of oliginal cape
2823  !! 2. ncount = 10
2824  if( (f_cape <= 1.05-stab .and. f_cape >= 0.95-stab) .or. ncount == 10) then
2825  exit
2826  else ! no exit
2827  if(ncount > 10) exit ! sayfty ?? ncount musn't over 10...
2828  if(f_cape == 0._rp) then ! f_cape is 0 -> new cape is 0 : too reducted
2829  ainc = ainc*0.5_rp
2830  else
2831  if(dcape < 1.e-4) then ! too small cape then exit at next loop(iter loop) step
2832  noiter = 1
2833  ainc = aincold
2834  cycle
2835  else ! calculate new factor ainc
2836  ainc = ainc*stab*cape/dcape
2837  end if
2838  end if
2839  ainc = min(aincmx,ainc) ! ainc must be less than aincmx
2840  !! if ainc becomes very small, effects of convection will be minimal so just ignore it
2841  if (ainc < 0.05) then !! noconvection
2842  i_convflag = 2
2843  return
2844  end if
2845  !!update valuables use factar ainc
2846  tder = tder2*ainc
2847  prcp_flux = prcp_flux2*ainc
2848  do kk = ks,k_top
2849  umf(kk) = umf2(kk)*ainc
2850  upent(kk) = upent2(kk)*ainc
2851  updet(kk) = updet2(kk)*ainc
2852  qcdet(kk) = qcdet2(kk)*ainc
2853  qidet(kk) = qidet2(kk)*ainc
2854  dmf(kk) = dmf2(kk)*ainc
2855  downent(kk) = downent2(kk)*ainc
2856  downdet(kk) = downdet2(kk)*ainc
2857  end do
2858  ! go back up for another iteration
2859  end if
2860 
2861  end do ! iter(ncount)
2862  ! get the cloud fraction
2863  cldfrac_kf(:,:) = 0._rp
2864  if (i_convflag == 1) then
2865  do kk = k_lcl-1, k_top
2866  umf_tmp = umf(kk)/(deltax*deltax)
2867  xcldfrac = 0.07*log(1._rp+(500._rp*umf_tmp))
2868  xcldfrac = max(1.e-2_rp,xcldfrac)
2869  cldfrac_kf(kk,1) = min(2.e-2_rp,xcldfrac) ! shallow
2870  end do
2871  else
2872  do kk = k_lcl-1, k_top
2873  umf_tmp = umf(kk)/(deltax*deltax)
2874  xcldfrac = 0.14*log(1._rp+(500._rp*umf_tmp))
2875  xcldfrac = max(1.e-2_rp,xcldfrac)
2876  cldfrac_kf(kk,2) = min(6.e-1_rp,xcldfrac) ! deep
2877  end do
2878  end if
2879  !!
2880  !!...compute hydrometeor tendencies as is done for T,qv...
2881  !!
2882  !! ...frc2 is the fraction of total condensate
2883  !!... generated that goes into precipitation
2884  !!
2885  !! Redistribute hydrometerors according to the final mass-flux variables
2886  !! dtheta/dt and dqv/dt is
2887  if (cpr > 0._rp) then
2888  frc2 = prcp_flux/(cpr*ainc)
2889  else
2890  frc2 = 0._rp
2891  end if
2892  !! no qc qi qr qs inputted in KF scheme
2893  qc_nw(ks:ke) = 0._rp
2894  qi_nw(ks:ke) = 0._rp
2895  qr_nw(ks:ke) = 0._rp
2896  qs_nw(ks:ke) = 0._rp
2897  do kk = ks,k_top
2898  rainfb(kk) = flux_qr(kk)*ainc*fbfrc*frc2
2899  snowfb(kk) = flux_qs(kk)*ainc*fbfrc*frc2
2900  end do
2901 
2902  do ntimecount = 1, nstep ! same as T, QV
2903  do kk = ks, k_top ! initialize
2904  qc_fxin(kk) = 0._rp
2905  qc_fxout(kk) = 0._rp
2906  qi_fxin(kk) = 0._rp
2907  qi_fxout(kk) = 0._rp
2908  qr_fxin(kk) = 0._rp
2909  qr_fxout(kk) = 0._rp
2910  qs_fxin(kk) = 0._rp
2911  qs_fxout(kk) = 0._rp
2912  end do
2913  do kk = ks+1,k_top ! calc flux variable
2914  if(omg(kk) <= 0._rp) then
2915  qc_fxin(kk) = -fxm(kk)*qc_nw(kk-1)
2916  qi_fxin(kk) = -fxm(kk)*qi_nw(kk-1)
2917  qr_fxin(kk) = -fxm(kk)*qr_nw(kk-1)
2918  qs_fxin(kk) = -fxm(kk)*qs_nw(kk-1)
2919  !
2920  qc_fxout(kk-1) = qc_fxout(kk-1) + qc_fxin(kk)
2921  qi_fxout(kk-1) = qi_fxout(kk-1) + qi_fxin(kk)
2922  qr_fxout(kk-1) = qr_fxout(kk-1) + qr_fxin(kk)
2923  qs_fxout(kk-1) = qs_fxout(kk-1) + qs_fxin(kk)
2924  else
2925  qc_fxout(kk) = fxm(kk)*qc_nw(kk)
2926  qi_fxout(kk) = fxm(kk)*qi_nw(kk)
2927  qr_fxout(kk) = fxm(kk)*qr_nw(kk)
2928  qs_fxout(kk) = fxm(kk)*qs_nw(kk)
2929  !
2930  qc_fxin(kk-1) = qc_fxin(kk-1) + qc_fxout(kk)
2931  qi_fxin(kk-1) = qi_fxin(kk-1) + qi_fxout(kk)
2932  qr_fxin(kk-1) = qr_fxin(kk-1) + qr_fxout(kk)
2933  qs_fxin(kk-1) = qs_fxin(kk-1) + qs_fxout(kk)
2934  end if
2935  end do
2936 
2937  do kk = ks, k_top
2938  qc_nw(kk) = qc_nw(kk) + (qc_fxin(kk) - qc_fxout(kk) + qcdet(kk) )*deltat*emsd(kk)
2939  qi_nw(kk) = qi_nw(kk) + (qi_fxin(kk) - qi_fxout(kk) + qidet(kk) )*deltat*emsd(kk)
2940  qr_nw(kk) = qr_nw(kk) + (qr_fxin(kk) - qr_fxout(kk) + rainfb(kk) )*deltat*emsd(kk)
2941  qs_nw(kk) = qs_nw(kk) + (qs_fxin(kk) - qs_fxout(kk) + snowfb(kk) )*deltat*emsd(kk)
2942  end do
2943  ! in original qlg= qlpa but it is not nessesary
2944  end do
2945 
2946  !! cumulus parameterization rain (rainrate_cp)and rain rate (rainratecp) is detern
2947  rainrate_cp = prcp_flux*(1._rp - fbfrc)/(deltax*deltax) ! if shallow convection then fbfrc = 1. -> noprcpitation
2948  !! evaluate moisuture budget
2949  qinit = 0._rp ! initial qv
2950  qfinl = 0._rp ! final water subsidence qv,qi,qr,qs...
2951  dpth = 0._rp !
2952  do kk = ks,k_top
2953  dpth = dpth + deltap(kk)
2954  qinit = qinit + qv(kk)*ems(kk)
2955  qfinl = qfinl + qv_g(kk)*ems(kk) ! qv
2956  qfinl = qfinl + (qc_nw(kk) + qi_nw(kk) + qr_nw(kk) + qs_nw(kk))*ems(kk)
2957  end do
2958  qfinl = qfinl + prcp_flux*timecp*(1._rp - fbfrc)
2959  err = (qfinl - qinit )*100._rp/qinit
2960  if (abs(err) > 0.05_rp .and. istop == 0) then
2961  ! write error message
2962  ! moisture budget error
2963  istop = 1
2964  write(*,*) "XXXX ERROR@KF,MOISTURE"
2965  call prc_mpistop
2966  end if
2967  !! feed back to resolvable scale tendencies
2968  !! if the advective time period(time_advec) is less than specified minimum
2969  !! timec, allow feed back to occur only during time_advec
2970  !! calc tendency of Theta_,qv, qc, qi, qr, qs
2971  !! temp_g ,qv_g, qc_nw,qi_nw,qr_nw,qs_nw
2972  !! calc new theta (SCALE theta)
2973  !! then return tendency
2974  !! warm rain
2975  if (warmrain) then ! assume Kessuler scheem?
2976  !!
2977  !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2978  !!
2979  do kk = ks,ke
2980  cpm = cp*(1._rp + 0.887_rp*qv_g(kk))
2981  temp_g(kk) = temp_g(kk) - (qi_nw(kk) + qs_nw(kk))*emelt/cpm
2982  qc_nw(kk) = qc_nw(kk) + qi_nw(kk)
2983  qr_nw(kk) = qr_nw(kk) + qs_nw(kk)
2984  qi_nw(kk) = 0._rp
2985  qs_nw(kk) = 0._rp
2986  end do
2987  return
2988  elseif(.not. flag_qs ) then
2989  !!
2990  !!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
2991  !!...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2992  !!
2993  do kk = ks,ke
2994  cpm = cp*(1._rp + 0.887*qv_g(kk))
2995  if(kk < k_ml) then
2996  temp_g(kk) = temp_g(kk) - (qi_nw(kk) + qs_nw(kk))*emelt/cpm
2997  elseif(kk > k_ml) then! kk == k_ml no melt
2998  temp_g(kk) = temp_g(kk) + (qi_nw(kk) + qs_nw(kk))*emelt/cpm
2999  end if
3000  qc_nw(kk) = qc_nw(kk) + qi_nw(kk)
3001  qr_nw(kk) = qr_nw(kk) + qs_nw(kk)
3002  qi_nw(kk) = 0._rp
3003  qs_nw(kk) = 0._rp
3004  end do
3005  return
3006  !!
3007  !!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
3008  !!...OF HYDROMETEORS DIRECTLY...
3009  !!
3010  elseif( flag_qs ) then
3011  if(.not.flag_qi) then
3012  do kk = ks, ke
3013  qs_nw(kk) = qs_nw(kk) + qi_nw(kk)
3014  end do
3015  end if
3016  return
3017  else ! not allow
3018  write(*,*) "xxx ERROR@KF,NOTallow namelist"
3019  write(*,*) "!!!!!Moiture setting error in KF check namelist"
3020  call prc_mpistop
3021  end if
3022  return
3023  end subroutine cp_kf_compensational
3024  !! calc exner function for potential temperature(contain water vapor)
3025  subroutine calcexn( pres,qv,exn) ! emanuel 1994 111pp ! check potential temperature definition
3026  use scale_const,only :&
3027  cp_dry => const_cpdry , &
3028  pre00 => const_pre00 , &
3029  r_dry => const_rdry
3030  implicit none
3031  real(RP),intent(in) :: pres
3032  real(RP),intent(in) :: qv
3033  real(RP),intent(out) :: exn
3034  exn = (pre00/pres)**(0.2854_rp*(1._rp - 0.28_rp*qv ))
3035  !! exect
3036  ! exn = (PRE00/pres)**(( R_dry + qv*R_vap)/(Cp_dry + qv*Cp_vap))
3037  return
3038  end subroutine calcexn
3039 !@!------------------------------------------------------------------------------
3040 !@! kf subroutine from WRF cood
3041 !@!------------------------------------------------------------------------------
3042  subroutine precipitation_oc1973( &
3043  QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,QNEWLQ, &
3044  QNEWIC,QLQOUT,QICOUT,G)
3045  !-----------------------------------------------------------------------
3046  implicit none
3047  !-----------------------------------------------------------------------
3048  ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
3049  ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
3050  ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
3051  ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
3052  ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
3053 
3054  real(RP), intent(in) :: G
3055  real(RP), intent(in) :: DZ,BOTERM,ENTERM!,RATE to be local variablebles
3056  real(RP), intent(inout) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
3057 
3058  real(RP) :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
3059 
3060  qtot=qliq+qice
3061  qnew=qnewlq+qnewic
3062  !
3063  ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
3064  ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
3065  ! LEVELS...
3066  !
3067  qest=0.5_rp*(qtot+qnew)
3068  g1=wtw+boterm-enterm-2._rp*g*dz*qest/1.5_rp
3069  IF(g1.LT.0.0)g1=0._rp
3070  wavg=0.5_rp*(sqrt(wtw)+sqrt(g1))
3071  conv=rate*dz/wavg ! KF90 Eq. 9
3072 
3073  !
3074  ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
3075  ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
3076  ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
3077  ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
3078  !
3079  ratio3=qnewlq/(qnew+1.e-8_rp)
3080  ! OLDQ=QTOT
3081  qtot=qtot+0.6_rp*qnew
3082  oldq=qtot
3083  ratio4=(0.6_rp*qnewlq+qliq)/(qtot+1.e-8_rp)
3084  qtot=qtot*exp(-conv) ! KF90 Eq. 9
3085  !
3086  ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
3087  ! PARCEL AT THIS LEVEL...
3088  !
3089  dq=oldq-qtot
3090  qlqout=ratio4*dq
3091  qicout=(1._rp-ratio4)*dq
3092  !
3093  ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
3094  ! LATE VERTICAL VELOCITY
3095  !
3096  pptdrg=0.5_rp*(oldq+qtot-0.2_rp*qnew)
3097  wtw=wtw+boterm-enterm-2._rp*g*dz*pptdrg/1.5_rp
3098  IF(abs(wtw).LT.1.e-4_rp)wtw=1.e-4_rp
3099  !
3100  ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
3101  ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
3102  !
3103  qliq=ratio4*qtot+ratio3*0.4_rp*qnew
3104  qice=(1._rp-ratio4)*qtot+(1._rp-ratio3)*0.4_rp*qnew
3105  qnewlq=0._rp
3106  qnewic=0._rp
3107  return
3108  end subroutine precipitation_oc1973
3109 
3110  ! Kessler type auto conversion
3111  subroutine precipitation_kessler( &
3112  QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,QNEWLQ, &
3113  QNEWIC,QLQOUT,QICOUT,G)
3114  implicit none
3115  ! kesseler type
3116  ! auto conversion
3117  real(RP), intent(in ) :: G
3118  real(RP), intent(in ) :: DZ,BOTERM,ENTERM!,RATE to be local variablebles
3119  real(RP), intent(inout) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
3120  real(RP) :: pptdrg
3121  real(RP) :: total_liq, total_ice
3122  ! parameter module value kf_threshold
3123  real(RP) :: auto_qc, auto_qi
3124  auto_qc = kf_threshold
3125  auto_qi = kf_threshold
3126 
3127  total_liq = qliq + qnewlq
3128  total_ice = qice + qnewic
3129 
3130  ! condensate in convective updraft is converted into precipitation
3131  qlqout = max( total_liq - auto_qc, 0.0_rp )
3132  qicout = max( total_ice - auto_qi, 0.0_rp )
3133 
3134  pptdrg = max( 0.5_rp * ( total_liq + total_ice - qlqout - qicout ), 0.0_rp )
3135  wtw=wtw+boterm-enterm-2._rp*g*dz*pptdrg/1.5_rp
3136  IF(abs(wtw).LT.1.e-4_rp)wtw=1.e-4_rp
3137 
3138  qliq = max( total_liq - qlqout, 0.0_rp )
3139  qlqout = total_liq - qliq
3140 
3141  qice = max( total_ice - qicout, 0.0_rp )
3142  qicout = total_ice - qice
3143 
3144  qnewlq=0.0_rp
3145  qnewic=0.0_rp
3146 
3147  return
3148  end subroutine precipitation_kessler
3149 
3150  ! ----------------------------------------------------------------------
3151  subroutine tpmix2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic )!,XLV1,XLV0)
3152  !
3153  ! Lookup table variables:
3154  ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
3155  ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3156  ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3157  ! REAL, SAVE, DIMENSION(1:200) :: ALU
3158  ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3159  ! End of Lookup table variables:
3160  !-----------------------------------------------------------------------
3161  implicit none
3162  !-----------------------------------------------------------------------
3163  real(RP), intent(in) :: P,THES!,XLV1,XLV0
3164  real(RP), intent(out) :: QNEWLQ,QNEWIC
3165  real(RP), intent(inout) :: TU,QU,QLIQ,QICE
3166 
3167  real(RP) :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
3168  real(RP) :: TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
3169  integer :: IPTB,ITHTB
3170  !-----------------------------------------------------------------------
3171 
3172  !c******** LOOKUP TABLE VARIABLES... ****************************
3173  ! parameter(kfnt=250,kfnp=220)
3174  !c
3175  ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
3176  ! * alu(200),rdpr,rdthk,plutop
3177  !C***************************************************************
3178  !c
3179  !c***********************************************************************
3180  !c scaling pressure and tt table index
3181  !c***********************************************************************
3182  !c
3183  tp=(p-plutop)*rdpr
3184  qq=tp-aint(tp)
3185  iptb=int(tp)+1
3186 
3187  !
3188  !***********************************************************************
3189  ! base and scaling factor for the
3190  !***********************************************************************
3191  !
3192  ! scaling the and tt table index
3193  bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
3194  tth=(thes-bth)*rdthk
3195  pp =tth-aint(tth)
3196  ithtb=int(tth)+1
3197  ! IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
3198  IF(iptb.GE.kfnp .OR. iptb.LE.1 .OR. ithtb.GE.250 .OR. ithtb.LE.1)THEN
3199  ! modify
3200  write(*,*) '**** OUT OF BOUNDS *********',iptb,ithtb,p,thes
3201  ! flush(98)
3202  ENDIF
3203  !
3204  t00=ttab(ithtb ,iptb )
3205  t10=ttab(ithtb+1,iptb )
3206  t01=ttab(ithtb ,iptb+1)
3207  t11=ttab(ithtb+1,iptb+1)
3208  !
3209  q00=qstab(ithtb ,iptb )
3210  q10=qstab(ithtb+1,iptb )
3211  q01=qstab(ithtb ,iptb+1)
3212  q11=qstab(ithtb+1,iptb+1)
3213  !
3214  !***********************************************************************
3215  ! parcel temperature
3216  !***********************************************************************
3217  !
3218  temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
3219  !
3220  qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
3221  !
3222  dq=qs-qu
3223  IF(dq.LE.0._rp)THEN
3224  qnew=qu-qs
3225  qu=qs
3226  ELSE
3227  !
3228  ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
3229  ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
3230  !
3231  qnew=0._rp
3232  qtot=qliq+qice
3233  !
3234  ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
3235  ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
3236  ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
3237  ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
3238  ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
3239  !
3240  !...subsaturated values only occur in calculations involving various mixtures of
3241  !...updraft and environmental air for estimation of entrainment and detrainment.
3242  !...For these purposes, assume that reasonable estimates can be given using
3243  !...liquid water saturation calculations only - i.e., ignore the effect of the
3244  !...ice phase in this process only...will not affect conservative properties...
3245  !
3246  IF(qtot.GE.dq)THEN
3247  qliq=qliq-dq*qliq/(qtot+1.e-10_rp)
3248  qice=qice-dq*qice/(qtot+1.e-10_rp)
3249  qu=qs
3250  ELSE
3251  rll=xlv0-xlv1*temp
3252  cpp=1004.5_rp*(1._rp+0.89_rp*qu)
3253  IF(qtot.LT.1.e-10_rp)THEN
3254  !
3255  !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
3256  temp=temp+rll*(dq/(1._rp+dq))/cpp
3257  ELSE
3258  !
3259  !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
3260  ! THE TEMPERATURE IS GIVEN BY:
3261  !
3262  temp=temp+rll*((dq-qtot)/(1._rp+dq-qtot))/cpp
3263  qu=qu+qtot
3264  qtot=0._rp
3265  qliq=0._rp
3266  qice=0._rp
3267  ENDIF
3268  ENDIF
3269  ENDIF
3270  tu=temp
3271  qnewlq=qnew
3272  qnewic=0._rp
3273  return
3274  end subroutine tpmix2
3275  !******************************************************************************
3276  subroutine dtfrznew(TU,P,THTEU,QU,QFRZ,QICE)!,ALIQ,BLIQ,CLIQ,DLIQ)
3277  !-----------------------------------------------------------------------
3278  use scale_precision
3279  use scale_atmos_saturation ,only :&
3280  atmos_saturation_psat_liq
3281  implicit none
3282  !-----------------------------------------------------------------------
3283  real(RP), intent(in) :: P,QFRZ!to module variable,ALIQ,BLIQ,CLIQ,DLIQ
3284  real(RP), intent(inout) :: TU,THTEU,QU,QICE
3285 
3286  real(RP) :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
3287  !-----------------------------------------------------------------------
3288  !
3289  !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
3290  !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
3291  !...TTFRZ TO TBFRZ...
3292  !...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
3293  !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
3294  !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
3295  !
3296 
3297  rlc=2.5e6_rp-2369.276_rp*(tu-273.16_rp)
3298  ! RLC=2.5E6_RP-2369.276_RP*(TU-273.15_RP) ! 273.16 -> 273.15 ??
3299  rls=2833922._rp-259.532_rp*(tu-273.16_rp)
3300  ! RLS=2833922._RP-259.532_RP*(TU-273.15_RP) ! 273.16 -> 273.15 ??
3301  rlf=rls-rlc
3302  cpp=1004.5_rp*(1._rp+0.89_rp*qu)
3303  !
3304  ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
3305  ! FOR SATURATION VAPOR PRESSURE...
3306  !
3307  a=(cliq-bliq*dliq)/((tu-dliq)*(tu-dliq))
3308  dtfrz = rlf*qfrz/(cpp+rls*qu*a)
3309  tu = tu+dtfrz
3310  ! call ATMOS_SATURATION_psat_liq(ES,TU) !saturation vapar pressure
3311 
3312  es = aliq*exp((bliq*tu-cliq)/(tu-dliq))
3313  qs = es*0.622_rp/(p-es)
3314  !
3315  !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
3316  !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
3317  !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
3318  !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
3319  !...TEMPERATURE TO THE SATURATION VARIABLE...
3320  !
3321  dqevap = qs-qu
3322  qice = qice-dqevap
3323  qu = qu+dqevap
3324  pii=(1.e5_rp/p)**(0.2854_rp*(1._rp-0.28_rp*qu))
3325  ! Bolton 1980
3326  ! Emanuel 1994 132p eq(4.7.9) pseudoequivalent PT
3327  thteu = tu*pii*exp((3374.6525_rp/tu - 2.5403_rp)*qu*(1._rp + 0.81_rp*qu))
3328  !
3329  end subroutine dtfrznew
3330  ! --------------------------------------------------------------------------------
3331  subroutine prof5(EQ,EE,UD)
3332  !
3333  !***********************************************************************
3334  !***** GAUSSIAN TYPE MIXING PROFILE....******************************
3335  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3336  ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
3337  ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
3338  ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
3339  ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
3340  ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
3341  ! JACK KAIN
3342  ! 7/6/89
3343  ! Solves for KF90 Eq. 2
3344  !
3345  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3346  !-----------------------------------------------------------------------
3347  implicit none
3348  !-----------------------------------------------------------------------
3349  real(RP), intent(in) :: EQ
3350  real(RP), intent(inout) :: EE,UD
3351 
3352  real(RP) :: SQRT2P,A1,A2,A3,P,SIGMA,FE
3353  real(RP) :: X,Y,EY,E45,T1,T2,C1,C2
3354 
3355  DATA sqrt2p,a1,a2,a3,p,sigma,fe/2.506628_rp,0.4361836_rp,-0.1201676_rp, &
3356  0.9372980_rp,0.33267_rp,0.166666667_rp,0.202765151_rp/
3357  x = (eq - 0.5_rp)/sigma
3358  y = 6._rp*eq - 3._rp
3359  ey = exp(y*y/(-2._rp))
3360  e45 = exp(-4.5_rp)
3361  t2 = 1._rp/(1._rp + p*abs(y))
3362  t1 = 0.500498_rp
3363  c1 = a1*t1+a2*t1*t1+a3*t1*t1*t1
3364  c2 = a1*t2+a2*t2*t2+a3*t2*t2*t2
3365  IF(y.GE.0._rp)THEN
3366  ee=sigma*(0.5_rp*(sqrt2p-e45*c1-ey*c2)+sigma*(e45-ey))-e45*eq*eq/2._rp
3367  ud=sigma*(0.5_rp*(ey*c2-e45*c1)+sigma*(e45-ey))-e45*(0.5_rp+eq*eq/2._rp- &
3368  eq)
3369  ELSE
3370  ee=sigma*(0.5_rp*(ey*c2-e45*c1)+sigma*(e45-ey))-e45*eq*eq/2._rp
3371  ud=sigma*(0.5_rp*(sqrt2p-e45*c1-ey*c2)+sigma*(e45-ey))-e45*(0.5_rp+eq* &
3372  eq/2._rp-eq)
3373  ENDIF
3374  ee=ee/fe
3375  ud=ud/fe
3376  end subroutine prof5
3377  ! ------------------------------------------------------------------------
3378  subroutine tpmix2dd(p,thes,ts,qs)!,i,j)
3379  !
3380  ! Lookup table variables:
3381  ! integer, PARAMETER :: (KFNT=250,KFNP=220)
3382  ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3383  ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3384  ! REAL, SAVE, DIMENSION(1:200) :: ALU
3385  ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3386  ! End of Lookup table variables:
3387  !-----------------------------------------------------------------------
3388  implicit none
3389  !-----------------------------------------------------------------------
3390  real(RP), intent(in) :: P,THES
3391  real(RP), intent(inout) :: TS,QS
3392 
3393  real(RP) :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
3394  integer :: IPTB,ITHTB
3395  !-----------------------------------------------------------------------
3396 
3397  !
3398  !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
3399  ! parameter(kfnt=250,kfnp=220)
3400  !
3401  ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
3402  ! alu(200),rdpr,rdthk,plutop
3403  !***************************************************************
3404  !
3405  !***********************************************************************
3406  ! scaling pressure and tt table index
3407  !***********************************************************************
3408  !
3409  tp=(p-plutop)*rdpr
3410  qq=tp-aint(tp)
3411  iptb=int(tp)+1
3412  !
3413  !***********************************************************************
3414  ! base and scaling factor for the
3415  !***********************************************************************
3416  !
3417  ! scaling the and tt table index
3418  bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
3419  tth=(thes-bth)*rdthk
3420  pp =tth-aint(tth)
3421  ithtb=int(tth)+1
3422  !
3423  t00=ttab(ithtb ,iptb )
3424  t10=ttab(ithtb+1,iptb )
3425  t01=ttab(ithtb ,iptb+1)
3426  t11=ttab(ithtb+1,iptb+1)
3427  !
3428  q00=qstab(ithtb ,iptb )
3429  q10=qstab(ithtb+1,iptb )
3430  q01=qstab(ithtb ,iptb+1)
3431  q11=qstab(ithtb+1,iptb+1)
3432  !
3433  !***********************************************************************
3434  ! parcel temperature and saturation mixing ratio
3435  !***********************************************************************
3436  !
3437  ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
3438  !
3439  qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
3440  !
3441  return
3442  end subroutine tpmix2dd
3443  !!-----------------------------------------------------------------------
3444  subroutine envirtht(P1,T1,Q1,THT1)!,ALIQ,BLIQ,CLIQ,DLIQ)
3445  !
3446  !-----------------------------------------------------------------------
3447  use scale_precision
3448  use scale_const,only : &
3449  p00 => const_pre00
3450  !! C1 ->
3451  implicit none
3452  !-----------------------------------------------------------------------
3453  real(RP), intent(in) :: P1,T1,Q1!,ALIQ,BLIQ,CLIQ,DLIQ module variables
3454  real(RP), intent(out) :: THT1
3455 
3456  real(RP) :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT
3457 ! real(RP) :: T00,P00,C1,C2,C3,C4,C5
3458  real(RP),parameter :: C1=3374.6525_rp
3459  real(RP),parameter :: C2=2.5403_rp
3460  !-----------------------------------------------------------------------
3461  ! DATA T00,P00,C1,C2,C3,C4,C5/273.16_RP,1.E5_RP,3374.6525_RP,2.5403_RP,3114.834_RP, &
3462  ! 0.278296_RP,1.0723E-3_RP/
3463  !
3464  ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
3465  !
3466  ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
3467  ! For example, KF90 Eq. 10 no longer used
3468  !
3469  ee=q1*p1/(0.622_rp+q1)
3470  ! TLOG=ALOG(EE/ALIQ)
3471  ! ...calculate LOG term using lookup table...
3472  !
3473 ! astrt=1.e-3_RP
3474 ! ainc=0.075_RP
3475 ! a1=ee/aliq
3476 ! tp=(a1-astrt)/ainc
3477 ! indlu=int(tp)+1
3478 ! value=(indlu-1)*ainc+astrt
3479 ! aintrp=(a1-value)/ainc
3480  ! tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
3481  ! change nouse lookuptable
3482  tlog = log(ee/aliq)
3483  !! Bolton(1980) Dew point temperature[K]
3484  tdpt=(cliq-dliq*tlog)/(bliq-tlog)
3485  !! Bolton(1980)
3486  tsat=tdpt - (0.212_rp+1.571e-3_rp*(tdpt-tem00)-4.36e-4_rp*(t1-tem00))*(t1-tdpt)
3487  ! TSAT = 2840._RP/(3.5_RP - log(ee) -4.805_RP) +55
3488  !! Bolton(1980) emanuel 132p (4.7.9)
3489  tht=t1*(p00/p1)**(0.2854_rp*(1._rp-0.28_rp*q1))
3490  tht1=tht*exp((c1/tsat-c2)*q1*(1._rp+0.81_rp*q1))
3491  !
3492  return
3493  end subroutine envirtht
3494  !***********************************************************************
3495  subroutine kf_lutab!(SVP1,SVP2,SVP3,SVPT0)
3496  !
3497  ! This subroutine is a lookup table.
3498  ! Given a series of series of saturation equivalent potential
3499  ! temperatures, the temperature is calculated.
3500  !
3501  !--------------------------------------------------------------------
3502  use scale_const, only :&
3503  cp => const_cpdry , &
3504  pre00 => const_pre00, &
3505  grav => const_grav
3506  IMPLICIT NONE
3507  !--------------------------------------------------------------------
3508  ! Lookup table variables
3509  ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
3510  ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3511  ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3512  ! REAL, SAVE, DIMENSION(1:200) :: ALU
3513  ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3514  ! End of Lookup table variables
3515  !! use scale_const, only : &
3516  !! SVPT0 -> CONST_TEM00,&
3517  integer :: KP,IT,ITCNT,I
3518  real(RP) :: DTH=1._rp,tmin=150._rp,toler=0.001_rp
3519  real(RP) :: PBOT,DPR, &
3520  TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
3521  ASTRT,AINC,A1,THTGS
3522  ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
3523  !! to module variables
3524  ! real(RP) :: ALIQ,BLIQ,CLIQ,DLIQ
3525  !! real(RP), INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3526  !
3527  ! equivalent potential temperature increment
3528  ! data dth/1._RP/
3529  ! minimum starting temp
3530  ! data tmin/150._RP/
3531  ! tolerance for accuracy of temperature
3532  ! data toler/0.001_RP/
3533  ! top pressure (pascals)
3534  plutop=5000.0_rp
3535  ! bottom pressure (pascals)
3536  pbot=110000.0_rp
3537  !! to module variable
3538  !! ALIQ = SVP1*1000.
3539  !! BLIQ = SVP2
3540  !! CLIQ = SVP2*SVPT0
3541  !! DLIQ = SVP3
3542 
3543  !
3544  ! compute parameters
3545  !
3546  ! 1._over_(sat. equiv. theta increment)
3547  rdthk=1._rp/dth
3548  ! pressure increment
3549  !
3550  dpr=(pbot-plutop)/REAL(kfnp-1)
3551  ! dpr=(pbot-plutop)/REAL(kfnp-1)
3552  ! 1._over_(pressure increment)
3553  rdpr=1._rp/dpr
3554  ! compute the spread of thes
3555  ! thespd=dth*(kfnt-1)
3556  !
3557  ! calculate the starting sat. equiv. theta
3558  !
3559  temp=tmin
3560  p=plutop-dpr
3561  do kp=1,kfnp
3562  p=p+dpr
3563  es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
3564  qs=0.622_rp*es/(p-es)
3565  pi=(1.e5_rp/p)**(0.2854_rp*(1.-0.28_rp*qs))
3566  the0k(kp)=temp*pi*exp((3374.6525_rp/temp-2.5403_rp)*qs* &
3567  (1._rp+0.81_rp*qs))
3568  enddo
3569  !
3570  ! compute temperatures for each sat. equiv. potential temp.
3571  !
3572  p=plutop-dpr
3573  do kp=1,kfnp
3574  thes=the0k(kp)-dth
3575  p=p+dpr
3576  do it=1,kfnt
3577  ! define sat. equiv. pot. temp.
3578  thes=thes+dth
3579  ! iterate to find temperature
3580  ! find initial guess
3581  if(it.eq.1) then
3582  tgues=tmin
3583  else
3584  tgues=ttab(it-1,kp)
3585  endif
3586  es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
3587  qs=0.622_rp*es/(p-es)
3588  pi=(1.e5_rp/p)**(0.2854_rp*(1._rp-0.28_rp*qs))
3589  thgues=tgues*pi*exp((3374.6525_rp/tgues-2.5403_rp)*qs* &
3590  (1._rp + 0.81_rp*qs))
3591  f0=thgues-thes
3592  t1=tgues-0.5_rp*f0
3593  t0=tgues
3594  itcnt=0
3595  ! iteration loop
3596  do itcnt=1,11
3597  es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
3598  qs=0.622_rp*es/(p-es)
3599  pi=(1.e5_rp/p)**(0.2854_rp*(1._rp-0.28_rp*qs))
3600  thtgs=t1*pi*exp((3374.6525_rp/t1-2.5403_rp)*qs*(1._rp + 0.81_rp*qs))
3601  f1=thtgs-thes
3602  if(abs(f1).lt.toler)then
3603  exit
3604  endif
3605  ! itcnt=itcnt+1
3606  dt=f1*(t1-t0)/(f1-f0)
3607  t0=t1
3608  f0=f1
3609  t1=t1-dt
3610  enddo
3611  ttab(it,kp)=t1
3612  qstab(it,kp)=qs
3613  enddo
3614  enddo
3615  !
3616  ! lookup table for tlog(emix/aliq)
3617  !
3618  ! set up intial variable for lookup tables
3619  !
3620  astrt=1.e-3_rp
3621  ainc=0.075_rp
3622  !
3623  a1=astrt-ainc
3624  do i=1,200
3625  a1=a1+ainc
3626  alu(i)=log(a1)
3627  enddo
3628  !GdCP is g/cp add for SCALE
3629  gdcp = - grav/cp ! inital set
3630  return
3631  end subroutine kf_lutab
3632  !
3633 end module scale_atmos_phy_cp_kf
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_grid::dy
real(rp), public dy
length in the main region [m]: y
Definition: scale_grid_cartesian.F90:43
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_atmos_phy_cp_kf::cp_kf_main
subroutine cp_kf_main(dens, MOMZ, MOMX, MOMY, RHOT, QTRC, w0avg, nca, DENS_t_CP, DT_RHOT, DT_RHOQ, rainrate_cp, cldfrac_sh, cldfrac_dp, timecp, cloudtop, zlcl, I_convflag)
Definition: scale_atmos_phy_cp_kf.F90:513
scale_grid::dx
real(rp), public dx
length in the main region [m]: x
Definition: scale_grid_cartesian.F90:42
scale_tracer::qqe
integer, public qqe
Definition: scale_tracer.F90:61
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
scale_atmos_phy_cp_kf::atmos_phy_cp_kf_setup
subroutine, public atmos_phy_cp_kf_setup(CP_TYPE)
Setup.
Definition: scale_atmos_phy_cp_kf.F90:134
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:95
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_real::real_fz
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
Definition: scale_grid_real.F90:41
scale_grid_real::real_cz
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
Definition: scale_grid_real.F90:40
scale_atmos_phy_cp_kf::kf_dt
real(rp), save, public kf_dt
Definition: scale_atmos_phy_cp_kf.F90:51
scale_atmos_phy_cp_kf::cp_kf_param
subroutine cp_kf_param(RATE_in, TRIGGER_in, FLAG_QS_in, FLAG_QI_in, KF_DT_in, DELCAPE_in, DEEPLIFETIME_in, SHALLOWLIFETIME_in, DEPTH_USL_in, KF_prec_in, KF_threshold_in, WARMRAIN_in, KF_LOG_in, stepkf_in)
Definition: scale_atmos_phy_cp_kf.F90:442
scale_atmos_phy_cp_kf::kf_lutab
subroutine kf_lutab
Definition: scale_atmos_phy_cp_kf.F90:3496
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_time::time_dtsec
real(dp), public time_dtsec
time interval of model [sec]
Definition: scale_time.F90:36
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_tracer::i_qv
integer, public i_qv
Definition: scale_tracer.F90:47
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:93
scale_tracer::i_mp2all
integer, dimension(:), allocatable, public i_mp2all
Definition: scale_tracer.F90:99
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_time
module TIME
Definition: scale_time.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_cp_kf
module ATMOSPHERE / Physics Cumulus Parameterization
Definition: scale_atmos_phy_cp_kf.F90:25
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:40
scale_const::const_emelt
real(rp), parameter, public const_emelt
Definition: scale_const.F90:74
scale_tracer::i_qs
integer, public i_qs
Definition: scale_tracer.F90:51
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:132
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_tracer::i_qi
integer, public i_qi
Definition: scale_tracer.F90:50
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_tracer::qqs
integer, public qqs
Definition: scale_tracer.F90:60
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_atmos_phy_cp_kf::atmos_phy_cp_kf
subroutine, public atmos_phy_cp_kf(DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, DENS_t_CP, MOMZ_t_CP, MOMX_t_CP, MOMY_t_CP, RHOT_t_CP, RHOQ_t_CP, MFLX_cloudbase, SFLX_convrain, cloudtop, cloudbase, cldfrac_dp, cldfrac_sh, nca, w0avg)
Definition: scale_atmos_phy_cp_kf.F90:298
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_atmos_phy_cp_kf::kf_wmean
subroutine kf_wmean(W0_avg, DENS, MOMZ)
running mean vertical wind speed
Definition: scale_atmos_phy_cp_kf.F90:389
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:178
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_tracer::i_qr
integer, public i_qr
Definition: scale_tracer.F90:49
scale_tracer::i_qc
integer, public i_qc
Definition: scale_tracer.F90:48
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13