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