1 !--------------------------------------------------------------------
2 ! Kain-Fritsch + CuP Cumulus Parameterization
5 ! kf_cup_cps* - the top-level driver routine
6 ! kf_cup_para* - the guts of the KF scheme
23 ! * = Subroutine either modified or added for CuP compared to the
24 ! original kfeta scheme.
25 !--------------------------------------------------------------------
27 !--------------------------------------------------------------------
29 ! - Add variable descriptions with units and other code docs
30 ! - Should we vary rBinSize based on t2 to get more sensitivity when cold?
31 ! - Figure out appropriate limiting values for the slopes and sigmas
32 ! that ensure the jfd sums to one and gives at least some
34 ! - Figure out how to make minimum frequency settings dependent upon
35 ! the chosen bin sizes.
36 ! - Tie cloud radius calc. to dx or the shallow trigger.
37 ! - When run with a small dx, deep convection should never be allowed
38 ! to trigger. Right now, it can.
39 ! - Figure out how to do cloud fraction feedback.
40 ! - Figure out how to handle combination of liquid and ice for cloud
41 ! fraction calculation.
42 ! - Clean up cldfratend_cup once we are sure that it will never be
44 ! - When fluxes are negative, wstar goes negative and then the
45 ! time scales go negative for tstar and taucloud. The neg. cancels
46 ! out for the cloud fraction, but it is troublesome none the less.
47 ! - Deep convective clouds don't necessarily develop concurrent
48 ! condensed phase mass. This has impacts for radiation and should
50 !--------------------------------------------------------------------
52 MODULE module_cu_kfcup
58 !--------------------------------------------------------------------
59 ! Lookup table variables:
60 INTEGER, PARAMETER, PRIVATE :: KFNT=250,KFNP=220
61 REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
62 REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
63 REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
64 REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
66 ! Note: KF Lookup table is used by subroutines KF_cup_PARA, TPMIX2,
68 ! End of Lookup table variables:
70 real, parameter, private :: eps=0.622 !used to be epsilon
71 !real, parameter, private :: reallysmall=1e-30 !for div by 0 checks
72 real, parameter, private :: reallysmall=5e-4 !for div by 0 checks
74 ! if ==1, apply barahona and nenes (2007) entrainment adjustment to activation
75 ! at cloud base ; if =/1, do not apply this
76 integer, parameter, private :: qndrop_cldbase_entrain_opt = 1
77 ! if ==1, updraft qndrop above cloud base is reduced by entrainment (dilution) ;
79 integer, parameter, private :: qndrop_incloud_entrain_opt = 1
80 ! minimum vertical velocity (m/s) passed to activate routine
81 real, parameter, private :: w_act_min = 0.2
83 ! for testing -- multiply aerosol number/volume by this before activation calculation
84 ! real, parameter, private :: naero_adjust_factor = 1.0
85 ! for testing -- if ==1, set aerosol size to dcen_sect for activation calcs
86 ! if /=1, do not adjust aerosol size
87 ! integer, parameter, private :: vaero_dsect_adjust_opt = 0
92 SUBROUTINE KF_cup_CPS( grid_id, & ! rce 10-may-2012
93 ids,ide, jds,jde, kds,kde &
94 ,ims,ime, jms,jme, kms,kme &
95 ,its,ite, jts,jte, kts,kte &
98 ,U,V,TH,T,W,dz8w,Pcps,pi &
99 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
100 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
101 ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
103 ,xland & !LD 18-Oct-2011
104 ,psfc,z,z_at_w,ht,tsk,hfx,qfx,mavail & !CuP, wig, 24-Aug-2006
105 ,sf_sfclay_physics & !CuP, wig, 24-Aug-2006
106 ,br,regime,pblh,kpbl,t2,q2 & !CuP, wig, 24-Aug-2006
107 ,slopeSfc,slopeEZ,sigmasfc,sigmaEZ & !CuP, wig, 24-Aug-2006
108 ,cupflag,cldfra_cup,cldfratend_cup & !CuP, wig, 18-Sep-2006
109 ,shall,taucloud,tactive & !CuP, wig, 18-Sep-2006
110 ,activeFrac & !CuP, lkb 5-May-2010
111 ,tstar, lnterms & !CuP, wig 4-Oct-2006
112 ,lnint & !CuP, wig 4-Oct-2006
113 ,numBins, thBinSize, rBinSize & !CuP, lkb 4-Nov-2009
114 ,minDeepFreq, minShallowFreq & !CuP, lkb 4-Nov-2009
115 ,wCloudBase & !CuP, lkb 4-April-2010
116 ,wact_cup & !CuP, rce 10-may-2012
117 ,wulcl_cup & !CuP, rce 10-may-2012
118 ,wup_cup & !CuP, rce 15-mar-2013
119 ,qc_ic_cup & !CuP, rce 10-may-2012
120 ,qndrop_ic_cup & !CuP, rce 10-may-2012
121 ,qc_iu_cup & !CuP, rce 10-may-2012
122 ,fcvt_qc_to_pr_cup & !CuP, rce 10-may-2012
123 ,fcvt_qc_to_qi_cup & !CuP, rce 10-may-2012
124 ,fcvt_qi_to_pr_cup & !CuP, rce 10-may-2012
125 ,mfup_cup & !CuP, rce 10-may-2012
126 ,mfup_ent_cup & !CuP, rce 10-may-2012
127 ,mfdn_cup & !CuP, rce 10-may-2012
128 ,mfdn_ent_cup & !CuP, rce 10-may-2012
129 ,updfra_cup & !CuP, rce 10-may-2012
130 ,tcloud_cup & !CuP, rce 10-may-2012
131 ,shcu_aerosols_opt & !CuP, rce 10-may-2012
133 ,chem_opt & !CuP, rce 10-may-2012
134 ,chem & !CuP, rce 10-may-2012
135 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
136 ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
140 USE module_state_description, only: num_chem
141 #if ( WRF_CHEM == 1 )
142 USE module_state_description, only: cbmz_mosaic_4bin, cbmz_mosaic_4bin_aq, &
143 cbmz_mosaic_8bin, cbmz_mosaic_8bin_aq, &
144 saprc99_mosaic_8bin_vbs2_aq_kpp, &
145 saprc99_mosaic_8bin_vbs2_kpp
146 USE module_data_mosaic_asect, only: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
147 ntype_aer, nsize_aer, ncomp_aer, &
148 ai_phase, msectional, massptr_aer, numptr_aer, &
149 dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer
152 !-------------------------------------------------------------
154 !-------------------------------------------------------------
155 INTEGER, INTENT(IN ) :: grid_id, & !rce 10-may-2012
156 ids,ide, jds,jde, kds,kde, &
157 ims,ime, jms,jme, kms,kme, &
158 its,ite, jts,jte, kts,kte
160 INTEGER, INTENT(IN ) :: STEPCU
161 LOGICAL, INTENT(IN ) :: warm_rain
163 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
164 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
165 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
167 INTEGER, INTENT(IN ) :: KTAU, &
168 sf_sfclay_physics, & !CuP, wig, 24-Aug-2006
169 shcu_aerosols_opt !CuP, rce, 10-may-2012
171 INTEGER, DIMENSION( ims:ime , jms:jme ) , & !CuP, wig, 24-Aug-2006
172 INTENT(IN ) :: & !CuP, wig, 24-Aug-2006
173 kpbl !Note that this is different from kpbl in the main KF scheme below. CuP, wig, 24-Aug-2006
175 REAL, DIMENSION( ims:ime , jms:jme ) , & !CuP, wig, 24-Aug-2006
176 INTENT(IN ) :: & !CuP, wig, 24-Aug-2006
177 psfc, & !CuP, wig, 24-Aug-2006
178 ht, & !CuP, wig, 24-Aug-2006
179 tsk, & !CuP, wig, 24-Aug-2006
180 hfx, & !CuP, wig, 24-Aug-2006
181 qfx, & !CuP, wig, 24-Aug-2006
182 mavail, & !CuP, wig, 24-Aug-2006
183 br, & !CuP, wig, 24-Aug-2006
184 regime, & !CuP, wig, 24-Aug-2006
185 pblh, & !CuP, wig, 24-Aug-2006
186 t2, & !CuP, wig, 24-Aug-2006
187 q2, & !CuP, wig, 24-Aug-2006
188 xland !LD 18-Oct-2011
190 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
202 z, & !CuP, wig, 24-Aug-2006
203 z_at_w !CuP, wig 5-Oct-2006
205 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
208 cldfra_cup, & !CuP, wig, 18-Sep-2006
209 cldfratend_cup !CuP, wig, 18-Sep-2006
211 REAL, INTENT(IN ) :: DT, DX
213 REAL, DIMENSION( ims:ime , jms:jme ), &
214 INTENT(INOUT) :: RAINCV
216 REAL, DIMENSION( ims:ime , jms:jme ), &
217 INTENT(INOUT) :: NCA, &
218 shall !CuP, wig, 18-Sep-2006 This has to be "real" because "integer" would only output zeros to the history file.
220 REAL, DIMENSION( ims:ime , jms:jme ), &
221 INTENT(OUT) :: CUBOT, &
223 slopeSfc, & !CuP, wig, 24-Aug-2006
224 slopeEZ, & !CuP, wig, 24-Aug-2006
225 sigmaSfc, & !CuP, wig, 24-Aug-2006
226 sigmaEZ, & !CuP, wig, 24-Aug-2006
227 taucloud, & !CuP, wig, 1-Oct-2006
228 tactive, & !CuP, wig, 1-Oct-2006
229 tstar, & !CuP, wig, 4-Oct-2006
230 lnint, & !CuP, wig, 4-Oct-2006
231 activeFrac, & !CuP, lkb, 5-May-2010
232 wCloudBase !CuP, lkb, 10-April-2010
234 REAL, DIMENSION( ims:ime , jms:jme ), &
235 INTENT(INOUT) :: wact_cup, & !CuP, rce 10-may-2012
236 wulcl_cup, & !CuP, rce 10-may-2012
237 tcloud_cup !CuP, rce 10-may-2012
239 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
241 wup_cup, & !CuP, rce 15-mar-2013
242 qc_ic_cup, & !CuP, rce 10-may-2012
243 qndrop_ic_cup, & !CuP, rce 10-may-2012
244 qc_iu_cup, & !CuP, rce 10-may-2012
245 fcvt_qc_to_pr_cup, & !CuP, rce 10-may-2012
246 fcvt_qc_to_qi_cup, & !CuP, rce 10-may-2012
247 fcvt_qi_to_pr_cup, & !CuP, rce 10-may-2012
248 mfup_cup, & !CuP, rce 10-may-2012
249 mfup_ent_cup, & !CuP, rce 10-may-2012
250 mfdn_cup, & !CuP, rce 10-may-2012
251 mfdn_ent_cup, & !CuP, rce 10-may-2012
252 updfra_cup !CuP, rce 10-may-2012
254 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
256 lnterms !CuP, wig 4-Oct-2006
258 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
259 INTENT(INOUT) :: CU_ACT_FLAG, &
260 cupflag !CuP, wig 9-Oct-2006
261 INTEGER, INTENT(IN) :: numBins
262 REAL, INTENT(IN) :: thBinSize, rBinSize
263 REAL, INTENT(IN) :: minDeepFreq, minShallowFreq
267 INTEGER, OPTIONAL, INTENT(IN ) :: chem_opt !CuP, rce 10-may-2012
269 REAL, DIMENSION( ims:ime , kms:kme, jms:jme, 1:num_chem ),&
270 OPTIONAL, INTENT(IN) :: &
271 chem !CuP, rce 10-may-2012
273 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
285 ! Flags relating to the optional tendency arrays declared above
286 ! Models that carry the optional tendencies will provdide the
287 ! optional arguments at compile time; these flags all the model
288 ! to determine at run-time whether a particular tracer is in
291 LOGICAL, OPTIONAL :: &
301 LOGICAL :: flag_qr, flag_qi, flag_qs
302 LOGICAL :: flag_chem ! rce 10-may-2012
304 REAL, DIMENSION( kts:kte ) :: &
308 th1d, & !wig, CuP, 24-Aug-2006
309 z1d, & !wig, CuP, 15-Sep-2006
310 z_at_w1d, & !wig, CuP 5-Oct-2006
316 cldfra_cup1d, & !wig, CuP, 20-Sep-2006
317 cldfratend_cup1d, & !wig, CuP, 20-Sep-2006
318 qndrop1d, & !rce, CuP, 11-may-2012
319 qc1d, & !rce, CuP, 11-may-2012
320 qi1d, & !rce, CuP, 11-may-2012
321 fcvt_qc_to_pr, & !rce, CuP, 11-may-2012
322 fcvt_qc_to_qi, & !rce, CuP, 11-may-2012
323 fcvt_qi_to_pr !rce, CuP, 11-may-2012
325 REAL, DIMENSION( kts:kte ):: &
333 REAL, DIMENSION( kts:kte, 1:num_chem ) :: &
334 chem1d !rce, CuP, 11-may-2012
336 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp,RTHCUMAX
339 INTEGER :: i,j,k,NTST,ICLDCK
341 ! Local vars specific to CuP... wig, 24-Aug-2006
342 !~sensitivity test for 41 integer, parameter :: numBins = 21 !Number of perturbations for each variable (theta & qvapor)
343 !! integer, parameter :: numBins = 41 !41!Number of perturbations for each variable (theta & qvapor)
344 !! ! Should be an odd value.
345 !! real, parameter :: thBinSize = 0.1 !0.1 !Size of potential temp. perturbation increment (K)
346 !! real, parameter :: rBinSize = 1.0e-4 !1e-4 !Size of mxing ratio perturbation increment (kg/kg)
347 ! real, parameter :: minFreq = 1e-5 !Minimum frequency required for a perturbation to be used ~should be dependent upon bin sizes
348 !! real, parameter :: minDeepFreq= 50e-2 !Cumulative freq. threshold before deep convection is allowed ~this was 5e-2 before
350 integer :: ipert, ishall, jpert, kcubot, kcutop
351 !!real :: activeFrac, biggestDeepFreq, cumDeepFreq, cumShallFreq, &
352 real :: biggestDeepFreq, cumDeepFreq, cumShallFreq, &
353 cubot_deep, cutop_deep, nca_deep, raincv_deep, &
354 cubot_shall, cutop_shall, nca_shall, raincv_shall, &
356 real, dimension(numBins) :: r_perturb, th_perturb
357 real, dimension(numBins, numBins) :: jfd
358 real, dimension(kts:kte) :: dqdt_deep, dqidt_deep, dqcdt_deep, &
359 dqrdt_deep, dqsdt_deep, dtdt_deep, &
360 dqdt_shall, dqidt_shall, dqcdt_shall, &
361 dqrdt_shall, dqsdt_shall, dtdt_shall, &
362 qlg, qlg_shall, qig, qig_shall
363 character(len=200) :: message
365 ! rce 11-may-2012 mods start -------------------------------------------
366 integer :: idiagee, idiagff
367 integer :: ipert_deepsv, jpert_deepsv
368 integer :: kcubotmin, kcubotmax, kcutopmin, kcutopmax
369 integer :: kupdrbot_deep, kupdrbot_shall
374 real :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg, tmph, tmpi, tmpj
375 real :: tmpr, tmps, tmpx, tmpy, tmpz
377 real :: tmp_nca, tmp_updfra
378 real :: tmpveca(1:999)
379 real :: updfra, updfra_deep, updfra_shall
380 real :: wact, wact_deep, wact_shall
381 real :: wcb_v2, wcb_v2_shall, wcb_v2_deep
382 real :: wulcl, wulcl_deep, wulcl_shall
383 real :: wcloudbase_shall, wcloudbase_deep
385 real, dimension(kts:kte) :: &
386 qlg_deep, qig_deep, &
387 qndrop_ic_deep, qndrop_ic_shall, &
388 qc_ic_deep, qc_ic_shall, &
389 qi_ic_deep, qi_ic_shall, &
390 fcvt_qc_to_pr_deep, fcvt_qc_to_pr_shall, &
391 fcvt_qc_to_qi_deep, fcvt_qc_to_qi_shall, &
392 fcvt_qi_to_pr_deep, fcvt_qi_to_pr_shall, &
394 umfout, uerout, udrout, &
395 umf_deep, uer_deep, udr_deep, &
396 umf_shall, uer_shall, udr_shall, &
397 dmfout, derout, ddrout, &
398 dmf_deep, der_deep, ddr_deep, & ! only deep has downdraft
399 wup, wup_deep, wup_shall
400 ! rce 11-may-2012 mods end ---------------------------------------------
404 !----------------------
410 IF ( PRESENT(F_QR) ) flag_qr = F_QR
411 IF ( PRESENT(F_QI) ) flag_qi = F_QI
412 IF ( PRESENT(F_QS) ) flag_qs = F_QS
414 ! flag_chem is .TRUE. only when chem is present, shcu_aerosols_opt >= 2, and chem_opt is appropriate
416 #if ( WRF_CHEM == 1 )
417 if ( PRESENT( chem ) .and. shcu_aerosols_opt >= 2) then
418 if ( chem_opt == cbmz_mosaic_4bin .or. &
419 chem_opt == cbmz_mosaic_4bin_aq .or. &
420 chem_opt == cbmz_mosaic_8bin .or. &
421 chem_opt == cbmz_mosaic_8bin_aq .or. &
422 chem_opt == saprc99_mosaic_8bin_vbs2_aq_kpp .or. &
423 chem_opt == saprc99_mosaic_8bin_vbs2_kpp ) then !BSINGH (04/08/2014): Added for non-aq vbs
427 CALL wrf_error_fatal( 'kf_cup_cps - bad chem_opt for shcu_aerosols_opt >= 2' )
432 idiagff = 0 ; idiagee = 0 ! rce 11-may-2012 start
433 if ((ide-ids <= 3) .and. (jde-jds <= 3)) then
434 idiagff = 1 ! turn on diagnostics at i=j=1 for single column runs
435 ! idiagff = 0 ! (do this to turn off extra diagnostics)
436 end if ! rce 11-may-2012 end
442 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
443 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
444 ! RHOE=Pcps(I,K,J)/(R*TV)
445 ! W0=-101.9368*SCR1/RHOE
446 W0=0.5*(w(I,K,J)+w(I,K+1,J)) !~this can probably be passed in instead of recalced
447 W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
448 ! CLDFRA_CUP(I,K,j) = 0.0 ! Start with 0 cloud fraction, added by LK Berg 10/29/09 01/11/2012
453 !...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
455 !----------------------
456 ICLDCK=MOD(KTAU,NTST)
458 ! rce 11-may-2012 mods start -------------------------------------------
459 if (idiagff > 0) then
461 write(*,'(a,i5,1p,4e11.3)') 'kfcup_control numbins, ...binsize, min...freq', numbins, thbinsize, rbinsize, mindeepfreq, minshallowfreq
462 write(*,'(a,3i5)') 'kfcup_control -- qndrop_cldbase_entrain_opt, ...incloud', &
463 qndrop_cldbase_entrain_opt, qndrop_incloud_entrain_opt
464 write(*,'(a,1p,2e11.35)') 'kfcup_control -- w_act_min', w_act_min
465 write(*,'(a,2i5/(a,3(i9,i5)))') 'kfcup_control -- grid_id, ktau', grid_id, ktau, &
466 'kfcup_control -- d indices', ids,ide, jds,jde, kds,kde, &
467 'kfcup_control -- m indices', ims,ime, jms,jme, kms,kme, &
468 'kfcup_control -- e indices', its,ite, jts,jte, kts,kte
471 write(*,'(a)') 'kfcup', 'kfcup', 'kfcup--------------------------------------------------------------------------------'
472 write(*,'(a,l5)') 'kfcup -- flag_chem', flag_chem
473 write(*,'(a,3i5,l5,3i5,f10.1,1p,2e10.2)') 'kfcup a00 ktau,ntst,icldck; cupflag,ishall,bot/top; nca,cldfra', &
474 ktau, ntst, icldck, cupflag(its,jts), nint(shall(its,jts)), nint(cubot(its,jts)), nint(cutop(its,jts)), nca(its,jts), &
475 maxval(cldfra_cup(its,kts:kte-2,jts)), maxval(rqvcuten(its,kts:kte-2,jts))
476 write(*,'(a,i5,1p,4e11.3)') 'kfcup numbins, ...binsize, min...freq', numbins, thbinsize, rbinsize, mindeepfreq, minshallowfreq
477 end if ! (idiagff > 0)
478 ! rce 11-may-2012 mods end ---------------------------------------------
480 if ((ide-ids <= 3) .and. (jde-jds <= 3)) then ! rce 11-may-2012
481 ! for single column, skip ktau=1
482 ltmpa = (ICLDCK .EQ. 0) .and. (KTAU .gt. 1)
484 ltmpa = (ICLDCK .EQ. 0) .or. (KTAU .eq. 1)
487 main_test_on_ktau_ntst: & ! rce 11-may-2012
489 ! IF(ICLDCK.EQ.0 .or. KTAU .eq. 1) then
492 !write(message,*)'~trying convection...'
493 !call wrf_message(message)
496 CU_ACT_FLAG(i,j) = .true.
500 main_loop_on_j: & ! rce 11-may-2012
502 main_loop_on_i: & ! rce 11-may-2012
505 idiagee = 0 ! rce 11-may-2012
506 if (idiagff > 0) then
507 ! turn on diagnostics at i=j=1 for single column runs
508 if (i==its .and. j==jts) idiagee = 1
511 ishall = int(shall(i,j)) !CuP, wig 19-Sep-2006
512 !write(message,*)'~i,j,nca,shall=',i,j,nca(i,j),ishall
513 !call wrf_message(message)
515 main_test_on_nca: & ! rce 11-may-2012
516 IF ( NCA(I,J) .ge. 0.5*DT ) then !byang 26 aug 2011
517 ! A previous call to KF triggered a cloud, and now we have to wait for
518 ! the appropriate time scale before triggering another cloud.
519 CU_ACT_FLAG(i,j) = .false.
522 !call wrf_message("~doing convection...")
535 qc_ic_cup(i,:,j) = 0.0 ! rce 11-may-2012 start
536 qndrop_ic_cup(i,:,j) = 0.0
537 qc_iu_cup(i,:,j) = 0.0
538 fcvt_qc_to_pr_cup(i,:,j) = 0.0
539 fcvt_qc_to_qi_cup(i,:,j) = 0.0
540 fcvt_qi_to_pr_cup(i,:,j) = 0.0
544 tcloud_cup(i,j) = 0.0
545 updfra_cup(i,:,j) = 0.0
546 mfup_cup(i,:,j) = 0.0
547 mfup_ent_cup(i,:,j) = 0.0
548 mfdn_cup(i,:,j) = 0.0
549 mfdn_ent_cup(i,:,j) = 0.0 ! rce 11-may-2012 end
551 ! assign vars from 3D to 1D
556 th1d(k) = th(i,k,j) !wig, CuP 24-Aug-2006
560 W0AVG1D(K) =W0AVG(I,K,J)
561 z1d(k) = z(i,k,j) !wig, CuP 15-Sep-2006
562 z_at_w1d(k) = z_at_w(i,k,j) !wig, CuP 15-Sep-2006
564 cldfra_cup1d(k) = cldfra_cup(i,k,j) !wig, CuP 20-Sep-2006
567 if ( flag_chem ) then ! rce 11-may-2012 start
570 chem1d(k,l) = chem(i,k,j,l)
584 jpert_deepsv = -999 ! rce 11-may-2012 end
586 ! CuP, wig: begin, Aug-2006
587 ! Get the slopes and std. dev. for CuP
590 !!$print*,dx, psfc(i,j), p1d, rho1d
591 !!$print*, dz1d, z1d, ht(i,j)
592 !!$print*, t1d, th1d, tsk(i,j), u1d, v1d
593 !!$print*, qv1d, hfx(i,j), qfx(i,j), mavail(i,j)
594 !!$print*, sf_sfclay_physics, br(i,j), regime(i,j), pblh(i,j)
595 !!$print*, kpbl(i,j), t2(i,j), q2(i,j)
596 !!$print*, slopeSfc(i,j), slopeEZ(i,j)
597 !!$print*, sigmaSfc(i,j), sigmaEZ(i,j)
598 !!$print*, wstar, cupflag(i,j)
599 !!$print*, kms, kme, kts, kte
601 !!$print*,'~entering cupSlopeSigma',i,j
603 call cupSlopeSigma(dx, psfc(i,j), p1d, rho1d, &
604 dz1d, z1d, ht(i,j), &
605 t1d, th1d, tsk(i,j), u1d, v1d, &
606 qv1d, hfx(i,j),xland(i,j), qfx(i,j), mavail(i,j), & ! add xland LD 19-Oct-2011
607 sf_sfclay_physics, br(i,j), regime(i,j), pblh(i,j),&
608 kpbl(i,j), t2(i,j), q2(i,j), &
609 slopeSfc(i,j), slopeEZ(i,j), &
610 sigmaSfc(i,j), sigmaEZ(i,j), &
611 wstar, cupflag(i,j), shall(i,j), &
614 if (idiagee>0) then ! rce 11-may-2012
615 write(*,'(a,l5,i5)') 'kfcup cupslopesigma cupflag, ishall', cupflag(i,j), nint(shall(i,j))
616 write(*,'(a,i10,1p,5e10.2)') 'kfcup kpbl, pblh, ht, z1d, dz', kpbl(i,j), pblh(i,j), ht(i,j), z1d(1), dz1d(1)
617 write(*,'(a, 1p,5e10.2)') 'kfcup hfx, qfx, regime // w0', hfx(i,j), qfx(i,j), regime(i,j)
618 write(*,'( 1p,10e10.2)') w0avg1d(kts:kts+19)
621 ! If the CuP scheme is activated, use the CuP perturbations.
622 ! Otherwise, default to the standard KF algorithm.
623 main_test_on_cupflag: & ! rce 11-may-2012
624 if( cupflag(i,j) ) then
626 ! Get the joint frequency distribution and the associated perturbations
628 !~The pert. calcs can be pulled out of the i/j do loops for speed, but
629 !~are left in right now in case we want to vary the pert. values based
630 !~on environmental conditions.
631 call cup_jfd(slopeSfc(i,j), slopeEZ(i,j), &
632 sigmaSfc(i,j), sigmaEZ(i,j), &
633 numBins, thBinSize, rBinSize, &
634 th_perturb, r_perturb, jfd )
636 ! Determine the minimum frequency of occurance that we will allow to
637 ! contribute to the results. This serves two purposes. It prevents large
638 ! excursions from the mean that might creep in from mal-conditioned
639 ! PBL structures. And, it also speeds up overall calculation time by
640 ! limiting which bins to send to the KF scheme for lifting.
642 minFreq = minShallowFreq*jfd(int(numBins/2)+1, int(numBins/2)+1)
643 !!minFreq = 1e-2*jfd(int(numBins/2)+1, int(numBins/2)+1)
644 if (idiagee>0) write(*,'(a,2i5,1p,2e11.3)') 'kfcup minfreq stuff', &
645 int(numBins/2)+1, int(numBins/2)+1, minshallowfreq, minfreq ! rce 11-may-2012
647 ! Setup some vars and then loop through all the perturbation
650 biggestDeepFreq = -999.
666 ! rce 11-may-2012 mods start -------------------------------------------
671 fcvt_qc_to_pr_shall = 0.
672 fcvt_qc_to_qi_shall = 0.
673 fcvt_qi_to_pr_shall = 0.
689 ! rce 11-may-2012 mods end ---------------------------------------------
692 PERTLOOPS: do jpert = 1,numBins
695 ! Only consider the perturbations that exceed a threshold value. Also,
696 ! skip this perturbation if we already know deep convection will be
697 ! output and the current probability is lower than a previous deep
698 ! convective possibility.
700 if( (jfd(ipert,jpert) < minFreq) .or. &
701 !!(jfd(ipert,jpert) > 0.001) .or. & ! lkb, 18-Aug-2008
702 !!(th_perturb(ipert) <= 0) .or. & ! lkb, 18-Aug-2008 : Commented out for tests run on 11/3/09 looking at lower freq. of DC
703 !!(r_perturb(ipert) <= 0) .or. & ! lkb, 18-Aug-2008 : COmmented out for tests run on 11/3/09
704 (cumDeepFreq > minDeepFreq .and. & ! lkb, 18-Aug_2008
705 jfd(ipert,jpert) < biggestDeepFreq) ) cycle
709 ! write(*,*) raincv ,'in if before KF_cup_PARA' !LD, 20-April-2011
710 if (idiagee>0) then ! rce 11-may-2012
711 write(*,'(a,2i5,1p,2e11.3)') 'kfcup calling kf_cup_para'
712 write(98,'(///a,i5,2i5,5x,a,2i5,1pe11.3)') 'kfcup calling kf_cup_para, ktau, i, j', ktau, i, j, &
713 'ijpert, jdf', ipert, jpert, jfd(ipert,jpert)
716 CALL KF_cup_PARA( GRID_ID, KTAU, & ! rce 11-may-2012
718 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
719 W0AVG1D,DT,DX,DXSQ,RHO1D, &
720 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
721 EP2,SVP1,SVP2,SVP3,SVPT0, &
722 pblh(i,j),z_at_w1d,cupflag(i,j), & !wig, 21-Feb-2008
723 th_perturb(ipert),r_perturb(jpert), & !wig, 25-Aug-2006
724 jfd(ipert,jpert), & !lkb, 15-Aug-2008
725 ishall,qlg,qig, & !wig, 20-Sep-2006
726 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
727 RAINCV,NCA,NTST, & !LD add PRATEC 21-April-2011
728 flag_QI,flag_QS,warm_rain, &
730 ids,ide, jds,jde, kds,kde, &
731 ims,ime, jms,jme, kms,kme, &
732 its,ite, jts,jte, kts,kte, &
733 idiagee, updfra, wulcl, wup, &
734 umfout, uerout, udrout, & ! rce 11-may-2012
735 dmfout, derout, ddrout, & ! "
736 shcu_aerosols_opt, & ! "
737 flag_chem, num_chem, & ! "
738 wact, qndrop1d, qc1d, qi1d, & ! "
739 fcvt_qc_to_qi, fcvt_qc_to_pr, & ! "
740 fcvt_qi_to_pr, chem1d, & ! "
741 #if ( WRF_CHEM == 1 )
742 maxd_acomp, maxd_aphase, & ! "
743 maxd_atype, maxd_asize, & ! "
744 ntype_aer, nsize_aer, ncomp_aer, & ! "
745 ai_phase, msectional, & ! "
746 massptr_aer, numptr_aer, & ! "
747 dlo_sect, dhi_sect, & ! "
748 dens_aer, hygro_aer, sigmag_aer ) ! "
751 1, 1 ) ! rce 11-may-2012
754 if (idiagee>0) then ! rce 11-may-2012
755 if (ishall==0 .or. ishall==1) then
756 write(*,'(a,3i5,1p,e11.3,a)') 'kfcup 1 ishall, cubot/top, nca', &
757 ishall, nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j), ' triggered'
759 write(*,'(a,3i5,1p,e11.3,a)') 'kfcup 1 ishall, cubot/top, nca', &
760 ishall, nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j)
764 ! if( raincv(i,j).ne. 0 ) then &!LD, 20-April-2011
765 ! write(*,*) raincv,'after cup_PARA' !LD, 20-April-2011
766 ! end if &!LD, 20-April-2011
768 ! Move these tendency applications to after the averaging for all the
769 ! different CuP perturbations.
770 !!$ IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
772 !!$ RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
773 !!$! RTHCUMAX=max(abs(RTHCUTEN(I,K,J)),RTHCUMAX)
774 !!$ RQVCUTEN(I,K,J)=DQDT(K)
778 ! If deep convection triggered, accumulate the deep convective
779 ! probability. This will also be the case if no convection occurred
780 ! and then the results would be small. Save the results if this deep
781 ! possibility is more probable than previous possibilities...
783 if( ishall == 0 ) then
784 cumDeepFreq = cumDeepFreq + jfd(ipert,jpert)
785 if( jfd(ipert,jpert) > biggestDeepFreq ) then
786 biggestDeepFreq = jfd(ipert,jpert)
787 do k = kts, kte ! Added by lkb
788 dqdt_deep(k) = dqdt(k)
789 dqidt_deep(k) = dqidt(k)
790 dqcdt_deep(k) = dqcdt(k)
791 dqrdt_deep(k) = dqrdt(k)
792 dqsdt_deep(k) = dqsdt(k)
793 dtdt_deep(k) = dtdt(k)
796 raincv_deep = raincv(i,j)
797 cubot_deep = cubot(i,j)
798 cutop_deep = cutop(i,j)
800 ipert_deepsv = ipert ! rce 11-may-2012 start
804 qndrop_ic_deep = qndrop1d
807 fcvt_qc_to_pr_deep = fcvt_qc_to_pr
808 fcvt_qc_to_qi_deep = fcvt_qc_to_qi
809 fcvt_qi_to_pr_deep = fcvt_qi_to_pr
814 wcb_v2_deep = max( wlcl, wulcl )
815 wcloudbase_deep = wlcl
817 kcubot = nint(cubot_deep)
818 kupdrbot_deep = kcubot
819 do k = kcubot-1, kts, -1
820 if ((umfout(k) > 0.0) .or. (uerout(k) > 0.0)) kupdrbot_deep = k
823 umf_deep(k) = max( 0.0, umfout(k) )
824 uer_deep(k) = max( 0.0, uerout(k) )
825 udr_deep(k) = max( 0.0, udrout(k) )
826 dmf_deep(k) = min( 0.0, dmfout(k) )
827 der_deep(k) = max( 0.0, derout(k) )
828 ddr_deep(k) = max( 0.0, ddrout(k) )
829 enddo ! rce 11-may-2012 end
832 ! Or if shallow convection ocurred and we need to accumulate
833 ! frequency weighted running sums of the results...
835 else if( ishall == 1 ) then
836 cumShallFreq = cumShallFreq + jfd(ipert,jpert) ! lkb-9/02/08 changed to just use JFD
837 !!dqdt_shall = dqdt_shall + dqdt*jfd(ipert,jpert)
839 do k = kts, kte ! Added by lkb
840 !!!dqdt_shall = dqdt_shall + dqdt*jfd(ipert,jpert)
841 dqdt_shall(k) = dqdt_shall(k) + dqdt(k)
842 !!!dqidt_shall = dqidt_shall + dqidt*jfd(ipert,jpert)
843 dqidt_shall(k) = dqidt_shall(k) + dqidt(k)
844 !!!dqcdt_shall = dqcdt_shall + dqcdt*jfd(ipert,jpert)
845 dqcdt_shall(k) = dqcdt_shall(k) + dqcdt(k)
846 !!!dqrdt_shall = dqrdt_shall + dqrdt*jfd(ipert,jpert)
847 dqrdt_shall(k) = dqrdt_shall(k) + dqrdt(k)
848 !!!dqsdt_shall = dqsdt_shall + dqsdt*jfd(ipert,jpert)
849 dqsdt_shall(k) = dqsdt_shall(k) + dqsdt(k)
850 !!!dtdt_shall(k) = dtdt_shall(k) + dtdt(k)*jfd(ipert,jpert)
851 dtdt_shall(k) = dtdt_shall(k) + dtdt(k)
852 ! in kf_cup_para, when you have shallow conv,
853 ! ainc (and so updraft area and mass fluxes) get multiplied by jfd(ipert,jpert)
854 ! which is passed in as "freq"
855 ! thus the following variables that are averaged over perts
856 ! should not be weighted by jfd (same for updfra_shall)
857 umf_shall(k) = umf_shall(k) + max( 0.0, umfout(k) ) ! rce 11-may-2012 start
858 uer_shall(k) = uer_shall(k) + max( 0.0, uerout(k) )
859 udr_shall(k) = udr_shall(k) + max( 0.0, udrout(k) ) ! rce 11-may-2012 end
861 nca_shall = nca(i,j)!NINT(TIMEC_SHALL/DT)*DT ! add 01/11/2012 real(ntst)*DT !add dt 01/11/2012 All shallow clouds have a 40 min time scale per KF code.
862 raincv_shall = raincv_shall + raincv(i,j)*jfd(ipert,jpert)
863 !!!raincv_shall = raincv_shall + raincv(i,j)
864 cubot_shall = cubot_shall + z1d(nint(cubot(i,j)))*jfd(ipert,jpert) !Average the heights, then back out index later
865 cutop_shall = cutop_shall + z1d(nint(cutop(i,j)))*jfd(ipert,jpert) !ditto
866 qlg_shall = qlg_shall + qlg*jfd(ipert,jpert)
867 !!!qlg_shall = qlg_shall + qlg
868 qig_shall = qig_shall + qig*jfd(ipert,jpert)
869 !!!qig_shall = qig_shall + qig
870 ! wCloudBase(i,j) = wLCL * jfd(ipert,jpert) + wCloudBase(i,j) ! rce 11-may-2012 start
871 wCloudBase_shall= wLCL * jfd(ipert,jpert) + wCloudBase_shall
872 do k = max( kts, nint(cubot(i,j)) ), min( kte, nint(cutop(i,j)) )
873 ! these are "in cloud" values, so only do them for cubot <= k <= cutop
874 cumshallfreq1d(k) = cumshallfreq1d(k) + jfd(ipert,jpert)
875 qndrop_ic_shall(k) = qndrop_ic_shall(k) + qndrop1d(k)*jfd(ipert,jpert)
876 qc_ic_shall(k) = qc_ic_shall(k) + qc1d(k)*jfd(ipert,jpert)
877 qi_ic_shall(k) = qi_ic_shall(k) + qi1d(k)*jfd(ipert,jpert)
878 ! fcvt_qc_to_pr is fraction of qc converted to precip as air moves through the updraft layer
879 ! compute average as: sum( fcvt_qc_to_pr * qc1d * jfd ) / sum( qc1d * jfd )
880 fcvt_qc_to_pr_shall(k) = fcvt_qc_to_pr_shall(k) + fcvt_qc_to_pr(k)*qc1d(k)*jfd(ipert,jpert)
881 fcvt_qc_to_qi_shall(k) = fcvt_qc_to_qi_shall(k) + fcvt_qc_to_qi(k)*qc1d(k)*jfd(ipert,jpert)
882 fcvt_qi_to_pr_shall(k) = fcvt_qi_to_pr_shall(k) + fcvt_qi_to_pr(k)*qi1d(k)*jfd(ipert,jpert)
884 wup_shall = wup_shall + wup*jfd(ipert,jpert)
885 wact_shall = wact_shall + wact*jfd(ipert,jpert)
886 wulcl_shall = wulcl_shall + wulcl*jfd(ipert,jpert)
887 updfra_shall = updfra_shall + updfra
888 wcb_v2_shall = wcb_v2_shall + jfd(ipert,jpert)*max( wlcl, wulcl )
889 kcubotmin = min( kcubotmin, nint(cubot(i,j)) )
890 kcubotmax = max( kcubotmax, nint(cubot(i,j)) )
891 kcutopmin = min( kcutopmin, nint(cutop(i,j)) )
892 kcutopmax = max( kcutopmax, nint(cutop(i,j)) ) ! rce 11-may-2012 end
897 ! Otherwise, no convection occurred so do nothing.
902 ! Now that we know what kind of convection will occur, copy the
903 ! appropriate type, shallow or deep, into the output arrays that
904 ! KF normally expects. Shallow convection needs to be turned into
905 ! an average from a running sum.
907 ! write(*,*) 'raincv_deep',raincv_deep,ishall,'raincv_deep' !LD, 20-April-2011
909 main_test_on_deep_shall_freq: & ! rce 11-may-2012
910 if( cumDeepFreq > minDeepFreq ) then !Deep convection
913 do k = kts, kte ! Added by lkb
914 dqdt(k) = dqdt_deep(k)
915 dqidt(k) = dqidt_deep(k)
916 dqcdt(k) = dqcdt_deep(k)
917 dqrdt(k) = dqrdt_deep(k)
918 dqsdt(k) = dqsdt_deep(k)
919 dtdt(k) = dtdt_deep(k)
923 raincv(i,j) = raincv_deep
924 cubot(i,j) = cubot_deep
925 cutop(i,j) = cutop_deep
926 ! write(*,*) 'raincv',raincv,ishall,'raincv' !LD, 20-April-2011
928 qc_iu_cup(i,kts:kte,j) = qc_ic_deep(kts:kte) ! rce 11-may-2012 start
929 qc_ic_cup(i,kts:kte,j) = qc_ic_deep(kts:kte)
930 qndrop_ic_cup(i,kts:kte,j) = qndrop_ic_deep(kts:kte)
931 wup_cup(i,kts:kte,j) = wup_deep(kts:kte)
932 wact_cup(i,j) = wact_deep
933 wulcl_cup(i,j) = wulcl_deep
934 wCloudBase(i,j) = wCloudBase_deep
937 kcutop = nint(cutop_deep)
938 fcvt_qc_to_pr_cup(i,kts:kcutop,j) = fcvt_qc_to_pr_deep(kts:kcutop)
939 fcvt_qc_to_qi_cup(i,kts:kcutop,j) = fcvt_qc_to_qi_deep(kts:kcutop)
940 fcvt_qi_to_pr_cup(i,kts:kcutop,j) = fcvt_qi_to_pr_deep(kts:kcutop)
942 call adjust_mfentdet_kfcup( idiagee, grid_id, ktau, &
943 i, j, kts, kte, kcutop, ishall, &
944 umf_deep, uer_deep, udr_deep, dmf_deep, der_deep, ddr_deep )
946 ! mfup_ent_cup(k) is at center of layer k, and is 0 for k > kcutop
947 mfup_ent_cup(i,kts:kcutop,j) = uer_deep(kts:kcutop)
948 ! mfup_cup(k) is at bottom of layer k, and is 0 for k > kcutop
949 ! umf_deep(k) is at top of layer k
950 mfup_cup(i,kts+1:kcutop,j) = umf_deep(kts:kcutop-1)
951 mfdn_ent_cup(i,kts:kcutop,j) = der_deep(kts:kcutop)
952 mfdn_cup(i,kts+1:kcutop,j) = dmf_deep(kts:kcutop-1)
954 updfra_cup(i,kupdrbot_deep:kcutop,j) = updfra_deep
955 tcloud_cup(i,j) = nca_deep ! rce 11-may-2012 end
957 !main_test_on_deep_shall_freq: & ! rce 11-may-2012
958 else if( cumShallFreq > 0. ) then !Shallow convection
960 activeFrac(i,j) = cumShallFreq
962 do k = kts, kte ! Added by lkb
963 !!!dqdt = dqdt_shall / cumShallFreq
964 dqdt(k) = dqdt_shall(k)
965 !!!dqidt = dqidt_shall / cumShallFreq
966 dqidt(k) = dqidt_shall(k)
967 !!!dqcdt = dqcdt_shall / cumShallFreq
968 dqcdt(k) = dqcdt_shall(k)
969 !!!dqrdt = dqrdt_shall / cumShallFreq
970 dqrdt(k) = dqrdt_shall(k)
971 !!!dqsdt = dqsdt_shall / cumShallFreq
972 dqsdt(k) = dqsdt_shall(k)
973 !!!dtdt(k) = dtdt_shall(k) / cumShallFreq
974 dtdt(k) = dtdt_shall(k)
977 nca(i,j) = nca_shall ! shallow convection timescale is locked to convective time scale
978 raincv(i,j) = raincv_shall / cumShallFreq
979 !!!raincv(i,j) = raincv_shall
981 cubot_shall = cubot_shall / cumShallFreq !This gives the average height in AMSL
982 cutop_shall = cutop_shall / cumShallFreq !ditto
983 cubot(i,j) = findIndex(cubot_shall, z_at_w1d)-1 !Now, get the index of the level
984 cutop(i,j) = findIndex(cutop_shall, z_at_w1d)-1 !ditto
985 qlg = qlg_shall / cumShallFreq
987 qig = qig_shall / cumShallFreq
990 ! wCloudBase(i,j) = wCloudBase(i,j) / cumShallFreq ! rce 11-may-2012 start
991 wCloudBase_shall= wCloudBase_shall/ cumShallFreq
992 wCloudBase(i,j) = wCloudBase_shall
995 ! these are "in cloud" values
996 if (cumshallfreq1d(k) > 0.0) then
997 fcvt_qc_to_pr_shall(k) = fcvt_qc_to_pr_shall(k) / max( 1.0e-20, qc_ic_shall(k) )
998 fcvt_qc_to_qi_shall(k) = fcvt_qc_to_qi_shall(k) / max( 1.0e-20, qc_ic_shall(k) )
999 fcvt_qi_to_pr_shall(k) = fcvt_qi_to_pr_shall(k) / max( 1.0e-20, qi_ic_shall(k) )
1000 qndrop_ic_shall(k) = qndrop_ic_shall(k)/cumshallfreq1d(k)
1001 qc_ic_shall(k) = qc_ic_shall(k)/cumshallfreq1d(k)
1002 qi_ic_shall(k) = qi_ic_shall(k)/cumshallfreq1d(k)
1004 cumshallfreq1d(k) = cumshallfreq1d(k)/cumshallfreq
1006 wup_shall = wup_shall/cumshallfreq
1007 wact_shall = wact_shall/cumshallfreq
1008 wulcl_shall = wulcl_shall/cumshallfreq
1009 wcb_v2_shall = wcb_v2_shall / cumshallfreq
1010 wup_cup(i,kts:kte,j) = wup_shall(kts:kte)
1011 wact_cup(i,j) = wact_shall
1012 wulcl_cup(i,j) = wulcl_shall
1013 wcb_v2 = wcb_v2_shall
1015 kcubot = nint(cubot(i,j))
1016 kcutop = nint(cutop(i,j))
1017 ! qc_ic_cup(k) and qndrop_ic_cup(k) are at center of layer k, and are 0 for k > kcutop
1018 qc_ic_cup(i,kts:kcutop,j) = qc_ic_shall(kts:kcutop)
1019 qndrop_ic_cup(i,kts:kcutop,j) = qndrop_ic_shall(kts:kcutop)
1020 ! note: qc_ic_shall = qc1d from subr. kf_cup_para is really for updraft
1021 ! if an empirical "in cumulus" cloud-water is used for radiation,
1022 ! it should be put into qc_ic_cup, and used for cloud-chemistry too
1023 qc_iu_cup(i,kts:kcutop,j) = qc_ic_shall(kts:kcutop)
1024 ! for qc_ic_cup, use the value in module_ra_cam (1.0 g/kg)
1025 ! For shallow convection, use a representative condensate value based on
1026 ! observations from CHAPS (Oklahoma area) and Florida (Blyth et al. 2005)
1027 qc_ic_cup(i,kcubot:kcutop,j) = 1.0e-3
1029 fcvt_qc_to_pr_cup(i,kts:kcutop,j) = fcvt_qc_to_pr_shall(kts:kcutop)
1030 fcvt_qc_to_qi_cup(i,kts:kcutop,j) = fcvt_qc_to_qi_shall(kts:kcutop)
1031 fcvt_qi_to_pr_cup(i,kts:kcutop,j) = fcvt_qi_to_pr_shall(kts:kcutop)
1033 call adjust_mfentdet_kfcup( idiagee, grid_id, ktau, &
1034 i, j, kts, kte, kcutop, ishall, &
1035 umf_shall, uer_shall, udr_shall, dmfout, derout, ddrout )
1037 ! mfup_ent_cup(k) is at center of layer k, and is 0 for k > kcutop
1038 mfup_ent_cup(i,kts:kcutop,j) = uer_shall(kts:kcutop)
1039 ! mfup_cup(k) is at bottom of layer k, and is 0 for k > kcutop
1040 ! umf_shall(k) is at top of layer k
1041 mfup_cup(i,kts+1:kcutop,j) = umf_shall(kts:kcutop-1)
1043 kupdrbot_shall = kcubot
1044 do k = kcubot-1, kts, -1
1045 if ((umf_shall(k) > 0.0) .or. (uer_shall(k) > 0.0)) kupdrbot_shall = k
1047 updfra_cup(i,kupdrbot_shall:kcutop,j) = updfra_shall
1048 tcloud_cup(i,j) = nca_shall ! rce 11-may-2012 end
1050 !main_test_on_deep_shall_freq: & ! rce 11-may-2012
1053 activeFrac(i,j) = 0.
1062 cubot(i,j) = 1.! add 1 replace 0 LD 01/11/2012
1064 end if main_test_on_deep_shall_freq ! rce 11-may-2012
1066 if (idiagee>0) write(*,'(a,3i5,1p,3e11.3)') 'kfcup 2 ishall, cubot/top, nca', &
1067 ishall, nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j) ! rce 11-may-2012
1069 shall(i,j) = real(ishall)
1070 kcubot = int(cubot(i,j))
1071 kcutop = int(cutop(i,j))
1072 call cupCloudFraction(qlg, qig, qv1d, t1d, z1d, p1d, &
1073 kcubot, kcutop, ishall, wStar, wCloudBase(i,j), pblh(i,j), dt, &
1074 activeFrac(i,j), cldfra_cup1d, cldfratend_cup1d, &
1075 taucloud(i,j), tActive(i,j), tstar(i,j), lnterms(i,:,j), &
1077 kts, kte, mfup_cup(i,:,j)) ! add mfup_cup LD 06 29 2012
1080 cldfra_cup(i,k,j) = cldfra_cup1d(k)
1084 if (idiagee > 0) then
1085 call cu_kfcup_diagee01( & ! rce 11-may-2012
1086 ims, ime, jms, jme, kms, kme, kts, kte, &
1088 idiagee, idiagff, ishall, ktau, &
1089 kcubotmin, kcubotmax, kcutopmin, kcutopmax, &
1090 activefrac, cldfra_cup1d, &
1091 cubot, cutop, cumshallfreq1d, &
1092 ddr_deep, der_deep, dmf_deep, dt, dz1d, &
1093 fcvt_qc_to_pr_deep, fcvt_qc_to_qi_deep, fcvt_qi_to_pr_deep, &
1094 fcvt_qc_to_pr_shall, fcvt_qc_to_qi_shall, fcvt_qi_to_pr_shall, &
1095 nca_deep, nca_shall, p1d, pblh, &
1096 qc_ic_deep, qc_ic_shall, qi_ic_deep, qi_ic_shall, qndrop_ic_cup, rho1d, &
1097 tactive, taucloud, tstar, &
1098 udr_deep, udr_shall, uer_deep, uer_shall, umf_deep, umf_shall, &
1099 updfra_deep, updfra_shall, updfra_cup, &
1100 wact_cup, wcloudbase, wcb_v2, wcb_v2_shall, &
1101 wulcl_cup, wstar, z1d, z_at_w1d )
1105 !!$ write(message,'(2i4,a,f10.5,a,f10.5)') i,j," Frequencies: cumDeepFreq=",cumDeepFreq," cumShallFreq=",cumShallFreq
1106 !!$ call wrf_message(message)
1109 !main_test_on_cupflag ! rce 11-may-2012
1112 ! CuP did not trigger due to stable conditions so default to standard
1115 !!CALL KF_cup_PARA(I, J, &
1116 !! U1D,V1D,T1D,QV1D,P1D,DZ1D, &
1117 !! W0AVG1D,DT,DX,DXSQ,RHO1D, &
1118 !! XLV0,XLV1,XLS0,XLS1,CP,R,G, &
1119 !! EP2,SVP1,SVP2,SVP3,SVPT0, &
1120 !! pblh(i,j),z_at_w1d,cupflag(i,j), & !wig, 21-Feb-2008
1121 !! th_perturb(1),r_perturb(1), & !wig, 9-Oct-2006
1122 !! 0.01, & !lkb, 15-Aug-2008, replace mass flux with default
1123 !! ishall,qlg,qig, & !wig, 20-Sep-2006
1124 !! DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
1125 !! RAINCV,NCA,NTST, &
1126 !! flag_QI,flag_QS,warm_rain, &
1128 !! ids,ide, jds,jde, kds,kde, &
1129 !! ims,ime, jms,jme, kms,kme, &
1130 !! its,ite, jts,jte, kts,kte)
1132 CALL KF_cup_PARA( GRID_ID, KTAU, & ! rce 11-may-2012
1134 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
1135 W0AVG1D,DT,DX,DXSQ,RHO1D, &
1136 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
1137 EP2,SVP1,SVP2,SVP3,SVPT0, &
1138 pblh(i,j),z_at_w1d,cupflag(i,j), & !wig, 21-Feb-2008
1139 th_perturb(1),r_perturb(1), & !wig, 9-Oct-2006
1140 0.01, & !lkb, 15-Aug-2008, replace mass flux with default
1141 ishall,qlg,qig, & !wig, 20-Sep-2006
1142 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
1143 RAINCV,NCA,NTST, & !LD, add PRATEC 21-Apr-2011
1144 flag_QI,flag_QS,warm_rain, &
1146 ids,ide, jds,jde, kds,kde, &
1147 ims,ime, jms,jme, kms,kme, &
1148 its,ite, jts,jte, kts,kte, &
1149 idiagee, updfra, wulcl, wup, &
1150 umfout, uerout, udrout, & ! rce 11-may-2012
1151 dmfout, derout, ddrout, & ! "
1152 shcu_aerosols_opt, & ! "
1153 flag_chem, num_chem, & ! "
1154 wact, qndrop1d, qc1d, qi1d, & ! "
1155 fcvt_qc_to_qi, fcvt_qc_to_pr, & ! "
1156 fcvt_qi_to_pr, chem1d, & ! "
1157 #if ( WRF_CHEM == 1 )
1158 maxd_acomp, maxd_aphase, & ! "
1159 maxd_atype, maxd_asize, & ! "
1160 ntype_aer, nsize_aer, ncomp_aer, & ! "
1161 ai_phase, msectional, & ! "
1162 massptr_aer, numptr_aer, & ! "
1163 dlo_sect, dhi_sect, & ! "
1164 dens_aer, hygro_aer, sigmag_aer ) ! "
1167 1, 1 ) ! rce 11-may-2012
1170 !!shall(i,j) = real(ishall)
1172 !! cldfra_cup(i,k,j) = 0.
1175 ! rce 11-may-2012 *** currently, clouds produce by this call to kf_cup_para do not
1176 ! rce 11-may-2012 *** produce any "cup" diagnostics and do not used by chem_cup
1177 ! rce 11-may-2012 *** may want to change that eventually
1179 end if main_test_on_cupflag ! rce 11-may-2012
1182 ! This was moved from earlier in the routine...
1183 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
1185 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
1186 RQVCUTEN(I,K,J)=DQDT(K)
1191 IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
1194 RQRCUTEN(I,K,J)=DQRDT(K)
1195 RQCCUTEN(I,K,J)=DQCDT(K)
1198 ! This is the case for Eta microphysics without 3d rain field
1201 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
1206 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
1208 IF(PRESENT( rqicuten )) THEN
1211 RQICUTEN(I,K,J)=DQIDT(K)
1216 IF(PRESENT( rqscuten )) THEN
1219 RQSCUTEN(I,K,J)=DQSDT(K)
1224 if (idiagee>0) then ! rce 11-may-2012
1225 write(*,'(a,3i5,1p,3e11.3)') 'kfcup 3 ishall, cubot/top, nca', ishall, nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j)
1226 write(*,'(a,5i5,1p,3e11.3)') 'kfcup a08 ishall, i/jpert_deep, cubot/top, nca', ishall, &
1227 ipert_deepsv, jpert_deepsv, nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j)
1230 ENDIF main_test_on_nca ! rce 11-may-2012
1232 ENDDO main_loop_on_i ! rce 11-may-2012
1233 ENDDO main_loop_on_j ! rce 11-may-2012
1234 ENDIF main_test_on_ktau_ntst ! rce 11-may-2012
1236 ! write(*,*) 'end',raincv,ishall,'end' !LD, 20-April-2011
1238 if (idiagff > 0) then ! rce 11-may-2012
1240 write(*,'(a,i5,10x,l5,3i5,f10.1,1p,2e10.2)') 'kfcup a09 ktau; cupflag,ishall,bot/top; nca,cldfra,rqvcuten', &
1241 ktau, cupflag(i,j), nint(shall(i,j)), nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j), &
1242 maxval(cldfra_cup(i,kts:kte-2,j)), maxval(rqvcuten(i,kts:kte-2,j))
1243 write(*,'(a,10i5)') 'kfcup a10 maxlocs for cldfra_cup & rqvcuten', &
1244 maxloc(cldfra_cup(i,kts:kte-2,j)), maxloc(rqvcuten(i,kts:kte-2,j))
1245 write(*,'(a,i7,l5,3i5,2f10.1)') 'kfcup_a20 ktau, cupflag, ishall, bot/top, nca, tcloud', &
1246 ktau, cupflag(i,j), nint(shall(i,j)), nint(cubot(i,j)), nint(cutop(i,j)), nca(i,j), tcloud_cup(i,j)
1249 END SUBROUTINE KF_cup_CPS
1250 ! ****************************************************************************
1251 !-----------------------------------------------------------
1252 SUBROUTINE KF_cup_PARA ( GRID_ID, KTAU, & ! rce 11-may-2012
1254 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
1256 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
1257 EP2,SVP1,SVP2,SVP3,SVPT0, &
1258 pblh,z_at_w1d,cupflag, & !wig, 21-Feb-2008
1259 th_perturb,r_perturb, & !wig, 25-Aug-2006
1260 freq, & !lkb, 15-Aug-2008
1261 ishall,qlg,qig, & !wig, 25-Aug-2006
1262 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
1263 RAINCV,NCA,NTST, & !LD, add PRATEC 21-Apr-2011
1264 F_QI,F_QS,warm_rain, &
1265 CUTOP,CUBOT, wLCL, &
1266 ids,ide, jds,jde, kds,kde, &
1267 ims,ime, jms,jme, kms,kme, &
1268 its,ite, jts,jte, kts,kte, & ! rce 11-may-2012
1269 idiagee, updfra, wulcl, wup, & ! "
1270 umfout, uerout, udrout, & ! "
1271 dmfout, derout, ddrout, & ! "
1272 shcu_aerosols_opt, & ! "
1273 flag_chem, num_chem, & ! "
1274 wact, qndrop1d, qc1d, qi1d, & ! "
1275 fcvt_qc_to_qi, fcvt_qc_to_pr, & ! "
1276 fcvt_qi_to_pr, chem1d, & ! "
1277 maxd_acomp, maxd_aphase, & ! "
1278 maxd_atype, maxd_asize, & ! "
1279 ntype_aer, nsize_aer, ncomp_aer, & ! "
1280 ai_phase, msectional, & ! "
1281 massptr_aer, numptr_aer, & ! "
1282 dlo_sect, dhi_sect, & ! "
1283 dens_aer, hygro_aer, sigmag_aer ) ! rce 11-may-2012
1285 !-----------------------------------------------------------
1286 !***** The KF scheme that is currently used in experimental runs of EMCs
1287 !***** Eta model....jsk 8/00
1290 !-----------------------------------------------------------
1291 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1292 ims,ime, jms,jme, kms,kme, &
1293 its,ite, jts,jte, kts,kte, &
1295 GRID_ID, KTAU ! rce 11-may-2012
1296 ! ,P_QI,P_QS,P_FIRST_SCALAR
1298 LOGICAL, INTENT(IN ) :: F_QI, F_QS
1300 LOGICAL, INTENT(IN ) :: warm_rain, &
1301 cupflag !CuP, wig 9-Oct-2006
1303 REAL, DIMENSION( kts:kte ), &
1304 INTENT(IN ) :: U0, &
1312 z_at_w1d !wig, 21-Feb-2008
1314 REAL, INTENT(IN ) :: DT,DX,DXSQ, &
1315 pblh, & !wig, 21-Feb-2008
1316 th_perturb, r_perturb, & !wig, 25-Aug-2006
1317 freq !lkb, 15-Aug-2008
1320 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
1321 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
1324 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
1332 REAL, DIMENSION( ims:ime , jms:jme ), &
1333 INTENT(INOUT) :: NCA
1335 REAL, DIMENSION( ims:ime , jms:jme ), &
1336 INTENT(INOUT) :: RAINCV !LD, add PRATEC 21-Apr-2011
1338 integer, intent(out) :: ishall !wig, 25-Aug-2006 (was local before)
1339 real, intent(out) :: wLCL !lkb, 29-April-2010
1340 REAL, DIMENSION( kts:kte ), INTENT(OUT) :: &
1341 qlg, & !wig, 20-Sep-2006 (was local before)
1342 qig !wig, 20-Sep-2006 (was local before)
1343 REAL, DIMENSION( ims:ime , jms:jme ), &
1344 INTENT(OUT) :: CUBOT, &
1347 ! rce 11-may-2012 mods start -------------------------------------------
1348 INTEGER, INTENT(IN ) :: idiagee, &
1349 shcu_aerosols_opt, &
1352 LOGICAL, INTENT(IN ) :: flag_chem
1354 REAL, INTENT(OUT ) :: updfra, &
1358 REAL, DIMENSION( kts:kte ), &
1359 INTENT(INOUT) :: umfout, &
1367 REAL, DIMENSION( kts:kte ), &
1368 INTENT(INOUT) :: qndrop1d, &
1375 REAL, DIMENSION( kts:kte, 1:num_chem ), &
1376 INTENT(INOUT) :: chem1d
1378 INTEGER, INTENT(IN ) :: maxd_acomp, &
1383 INTEGER, INTENT(IN ), OPTIONAL :: ntype_aer, &
1384 nsize_aer(maxd_atype), &
1385 ncomp_aer(maxd_atype), &
1388 massptr_aer(maxd_acomp,maxd_asize,maxd_atype,maxd_aphase), &
1389 numptr_aer(maxd_asize,maxd_atype,maxd_aphase)
1391 REAL, DIMENSION( maxd_asize, maxd_atype ), &
1392 INTENT(IN ), OPTIONAL :: dlo_sect, dhi_sect, &
1395 REAL, DIMENSION( maxd_acomp, maxd_atype ), &
1396 INTENT(IN ), OPTIONAL :: dens_aer, hygro_aer
1397 ! rce 11-may-2012 mods end ---------------------------------------------
1400 !...DEFINE LOCAL VARIABLES...
1402 REAL, DIMENSION( kts:kte ) :: &
1403 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
1404 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
1405 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
1406 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
1407 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
1408 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
1409 DETLQ2,DETIC2,RATIO,RATIO2
1412 REAL, DIMENSION( kts:kte ) :: &
1413 DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
1414 QDT,FXM,THTAG,THPA,THFXOUT, &
1415 THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
1416 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
1417 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
1418 QL0,QI0,QR0,QRG,QS0,QSG
1421 REAL, DIMENSION( kts:kte+1 ) :: OMG
1422 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
1423 REAL, DIMENSION( kts:kte ) :: &
1424 CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
1428 REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
1430 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
1432 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
1433 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
1434 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
1435 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT, &
1436 !!ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
1437 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
1438 UPNEW,ABE,WKLCL,TTEMP,FRC1, &
1439 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
1440 DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
1441 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
1442 UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
1443 THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
1444 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
1445 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
1446 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
1447 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
1448 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
1449 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
1450 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
1451 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
1452 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
1453 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
1454 DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
1455 REAL :: TIMEC_SHALL ! Added by lkb, 10/31/10
1456 REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
1457 QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
1459 INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
1460 REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
1461 REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
1466 INTEGER, DIMENSION (kts:kte) :: KCHECK
1468 INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
1469 LC,MXLAYR,LLFC,NLAYRS,NK, &
1470 !KPBL,KLCL,LCL,LET,IFLAG, &
1471 KCLDLAYER,KLCL,LCL,LET,IFLAG, &
1472 NK1,LTOP,NJ,LTOP1, &
1473 LTOPM1,LVF,KSTART,KMIN,LFS, &
1474 ND,NIC,LDB,LDT,ND1,NDK, &
1475 NM,LMAX,NCOUNT,NOITR, &
1476 NSTEP,NTC,NCHM,NSHALL
1478 CHARACTER*1024 message
1480 INTEGER :: ksvaa ! rce 11-may-2012
1481 REAL :: rho_act, tk_act, w_act, w_act_eff ! rce 11-may-2012
1482 REAL :: qndrop_tmp ! rce 11-may-2012
1483 REAL :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg, tmph, tmpi
1484 REAL :: tmp_alphabn, tmp_ebn, tmp_escale, tmp_lv ! rce 11-may-2012
1485 REAL :: tmp_deltarh, tmp_deltatk, tmp_deltatkfact ! rce 11-may-2012
1486 REAL :: qndropbb(kts:kte) ! rce 11-may-2012
1489 DATA P00,T00/1.E5,273.16/
1491 DATA RHIC,RHBC/1.,0.90/
1492 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
1495 !-----------------------------------------------------------
1503 ! rce 11-may-2012 mods start -------------------------------------------
1504 if (idiagee > 0) IPRNT=.TRUE.
1521 ! rce 11-may-2012 mods end ---------------------------------------------
1534 !****************************************************************************
1536 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
1537 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
1538 !...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
1539 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
1540 FBFRC=0.0 ! PPT FB MODS
1541 !...mods to allow shallow convection...
1548 !... Set time constant for shallow convection
1549 TIMEC_SHALL = 1800.0 ! Set to the min value allowed for all convection
1552 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
1553 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
1554 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
1556 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
1557 !...FROM BOTTOM-UP IN THE KF SCHEME...
1560 !SUE tmprpsb=1./PSB(I,J)
1561 !SUE CELL=PTOP*tmprpsb
1565 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
1567 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
1568 QES(K)=0.622*ES/(P0(K)-ES)
1569 Q0(K)=AMIN1(QES(K),QV0(K))
1570 Q0(K)=AMAX1(0.000001,Q0(K))
1575 RH(K) = Q0(K)/QES(K)
1577 TV0(K)=T0(K)*(1.+0.608*Q0(K))
1578 ! RHOE(K)=P0(K)/(R*TV0(K))
1579 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
1580 DP(K)=rhoe(k)*g*DZQ(k)
1581 ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
1582 ! use it for shallow convection...For now, assume it is not available....
1583 ! TKE(K) = Q2(I,J,NK)
1586 ! IF(P0(K).GE.500E2)L5=K
1587 IF(P0(K).GE.0.5*P0(1))L5=K
1588 IF(P0(K).GE.P300)LLFC=K
1589 IF(T0(K).GT.T00)ML=K
1592 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
1596 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
1597 DZA(K-1)=Z0(K)-Z0(K-1)
1602 ! To save time, specify a pressure interval to move up in sequential
1603 ! check of different ~50 mb deep groups of adjacent model layers in
1604 ! the process of identifying updraft source layer (USL). Note that
1605 ! this search is terminated as soon as a buoyant parcel is found and
1606 ! this parcel can produce a cloud greater than specifed minimum depth
1607 ! (CHMIN)...For now, set interval at 15 mb...
1613 IF(P0(K).LT.PM15)THEN
1624 IF(NU.GT.NCHECK)THEN
1630 IF(CLDHGT(NNN).GT.CHMAX)THEN
1640 ! wig, 29-Aug-2006: I think this is where no convecion occurs. So, force
1641 ! ishall to a flag value to indicate this for accounting purposes.
1651 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
1652 !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
1653 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
1654 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
1659 IF ( NK+1 .LT. KTS ) THEN
1660 WRITE(message,*)'WOULD GO OFF BOTTOM: KF_CUP_PARA I,J,NK',I,J,NK
1661 CALL wrf_message (TRIM(message))
1665 IF ( NK .GT. KTE ) THEN
1667 'WOULD GO OFF TOP: KF_CUP_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
1668 CALL wrf_message (TRIM(message))
1671 DPTHMX=DPTHMX+DP(NK)
1673 IF(DPTHMX.GT.DPMIN)THEN
1678 IF(DPTHMX.LT.DPMIN)THEN
1679 ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
1684 KCLDLAYER=LC+NLAYRS-1 ! Added new veriable for top of cloud layer
1685 !!if(ishall .eq. 0) KPBL=LC !lkb, changed to only adjust mixed layer depth for deep convection
1687 !...********************************************************
1688 !...for computational simplicity without much loss in accuracy,
1689 !...mix temperature instead of theta for evaluating convective
1690 !...initiation (triggering) potential...
1697 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
1698 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
1704 TMIX=TMIX+DP(NK)*T0(NK)
1705 QMIX=QMIX+DP(NK)*Q0(NK)
1706 ZMIX=ZMIX+DP(NK)*Z0(NK)
1707 PMIX=PMIX+DP(NK)*P0(NK)
1709 ! THMIX=THMIX/DPTHMX
1714 EMIX=QMIX*PMIX/(0.622+QMIX)
1716 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
1718 ! TLOG=ALOG(EMIX/ALIQ)
1719 ! ...calculate dewpoint using lookup table...
1726 value=(indlu-1)*ainc+astrt
1727 aintrp=(a1-value)/ainc
1728 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
1729 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1730 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
1731 TLCL=AMIN1(TLCL,TMIX)
1732 TVLCL=TLCL*(1.+0.608*QMIX)
1733 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
1738 IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
1743 ! wig, 29-Aug-2006: Indicate no convection occurred.
1748 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
1750 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1752 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
1753 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
1754 TVEN=TENV*(1.+0.608*QENV)
1756 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
1757 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
1758 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
1759 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
1760 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
1761 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
1762 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
1763 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
1764 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
1765 IF(ZLCL.LT.2.E3)THEN
1766 WKLCL=0.02*ZLCL/2.E3
1770 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
1772 ! CuP, wig, 28-Aug-2006, begin:
1774 ! Replace KF perturbation temperatures with CuP perturbations. CuP
1775 ! perturbations are in potential temp. so convert the theta difference
1776 ! to a temperature difference. For the moisture perturbation, convert
1777 ! the CuP mixing ratio (kg/kg) into a virtual temperature adjustment.
1779 ! Standard KF way...
1780 if( .not. cupflag ) then
1781 IF(WKL.LT.0.0001)THEN
1784 DTLCL=4.64*WKL**0.33
1786 DTRH = 0. !CuP, wig: Move this from a few lines below since
1787 ! it is commented out there for CuP.
1790 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1791 dtlcl = th_perturb*(p00/p0(k))**rocp
1792 dtrh = 0.608*r_perturb
1797 !...for ETA model, give parcel an extra temperature perturbation based
1798 !...the threshold RH for condensation (U00)...
1800 !...for now, just assume U00=0.75...
1801 !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
1804 ! QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
1805 ! RHLCL = QENV/QSLCL
1806 ! DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
1807 ! IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
1808 ! DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
1809 ! ELSEIF(RHLCL.GT.0.95)THEN
1810 ! DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
1812 !!$wig, 28-Aug-2006 DTRH = 0.
1815 ! IF(ISHALL.EQ.1)IPRNT=.TRUE.
1817 ! IF(TLCL+DTLCL.GT.TENV)GOTO 45
1820 ! CuP, wig 28-Aug-2006, begin: Change parcel temperature adjustment
1821 ! comparison to use virtual temperature instead of "normal"
1823 !~Check to see if this should be switched back if cupflag==F. Why isn't
1824 ! the virt. temp. used in the standard scheme?
1825 !!$trigger: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
1826 TVLCL=TLCL*(1.+0.608*QMIX)
1827 trigger: if( tvlcl+dtlcl+dtrh < tven ) then
1831 ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
1835 ELSE ! Parcel is buoyant, determine updraft
1837 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
1838 !...EQUIVALENT POTENTIAL TEMPERATURE
1839 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
1841 CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
1843 !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
1845 ! CuP, wig 28-Aug-2006: The original KF algorithm sets the parcel's
1846 ! initial pert. vertical velocity at the LCL based on the pert.
1847 ! temperature, with a minimum W of 3. But, if the pert. temp. is
1848 ! negative, a smaller minimum positive W is set (==1). For CuP,
1849 ! allow the perturbation to set the W without any constraints
1850 ! except that the pert. must be positive.
1852 IF(DTTOT.GT.1.E-4)THEN
1853 GDT=2.*G*DTTOT*500./TVEN
1854 WLCL=1.+0.5*SQRT(GDT)
1855 if( .not. cupflag ) WLCL = AMIN1(WLCL,3.) !wig 9-Oct-2006
1863 !print*,'~ dttot and wlcl=',dttot,wlcl
1865 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1868 TVLCL=TLCL*(1.+0.608*QMIX)
1869 RHOLCL=PLCL/(R*TVLCL)
1873 ! make RAD a function of background vertical velocity...
1876 ELSEIF(WKL.GT.0.1)THEN
1879 RAD = 1000.+1000*WKL/0.1
1882 !*******************************************************************
1884 ! COMPUTE UPDRAFT PROPERTIES *
1886 !*******************************************************************
1890 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
1895 !!UMF(K)=freq*dxsq*WU(K)*RHOLCL ! Added by lkb
1899 ksvaa = k ! rce 11-may-2012
1901 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
1902 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
1903 !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
1924 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
1925 !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
1926 !...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
1927 !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
1928 !...PREVIOUS MODEL LEVEL...
1932 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
1933 !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
1934 !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
1941 qndropbb(:) = 0.0 ! rce 11-may-2012
1943 updraft: DO NK=K,KL-1
1945 RATIO2(NK1)=RATIO2(NK)
1948 THETEU(NK1)=THETEU(NK)
1952 call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
1953 qice(nk1),qnewlq,qnewic,XLV1,XLV0)
1956 !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
1957 !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
1958 !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
1959 !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
1960 !...LIQUID WATER IS FROZEN AT EACH LEVEL...
1962 IF(TU(NK1).LE.TTFRZ)THEN
1963 IF(TU(NK1).GT.TBFRZ)THEN
1964 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
1965 FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
1972 ! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
1973 !...IS BELOW TTFRZ...
1975 ! rce 11-may-2012 - added lines with tmpa/c and fcvt_qc_to_qi
1976 tmpa = max( 0.0, qliq(nk1)+qnewlq ) ! qliq before freezing calc
1977 QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
1978 QNEWIC=QNEWIC+QNEWLQ*FRC1
1979 QNEWLQ=QNEWLQ-QNEWLQ*FRC1
1980 QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
1981 QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
1982 tmpc = max( 0.0, qliq(nk1)+qnewlq ) ! qliq after freezing calc
1983 fcvt_qc_to_qi(nk1) = max( 0.0, tmpa-tmpc ) / max( 1.0e-10, tmpa )
1984 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
1985 QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1987 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
1989 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
1992 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
1993 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
1996 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
1997 BOTERM=2.*DZA(NK)*G*BE/1.5
2000 ENTERM=2.*REI*WTW/UPOLD
2002 ! rce 11-may-2012 - added lines with tmpa/b/c and fcvt_q?_to_pr
2003 tmpa = max( 0.0, qliq(nk1)+qnewlq ) ! qliq before precip calc
2004 tmpb = max( 0.0, qice(nk1)+qnewic ) ! qice before precip calc
2005 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
2006 RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
2007 tmpc = max( 0.0, qliq(nk1)+qnewlq ) ! qliq after precip calc
2008 fcvt_qc_to_pr(nk1) = max( 0.0, tmpa-tmpc ) / max( 1.0e-10, tmpa )
2009 tmpc = max( 0.0, qice(nk1)+qnewic ) ! qice after precip calc
2010 fcvt_qi_to_pr(nk1) = max( 0.0, tmpb-tmpc ) / max( 1.0e-10, tmpb )
2013 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
2014 !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
2016 IF(WTW.LT.1.E-3)THEN
2021 !...Calculate value of THETA-E in environment to entrain into updraft...
2023 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
2025 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
2027 REI=VMFLCL*DP(NK1)*0.03/RAD
2028 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
2030 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
2032 DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
2034 IF(DILBE.GT.0.)ABE=ABE+DILBE*G
2036 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
2037 !...ENTRAINMENT (0.5*REI) IS IMPOSED...
2039 IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
2047 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
2051 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
2052 QTMP=F1*Q0(NK1)+F2*QU(NK1)
2055 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
2056 qnewlq,qnewic,XLV1,XLV0)
2057 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
2058 IF(TU95.GT.TV0(NK1))THEN
2065 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
2066 QTMP=F1*Q0(NK1)+F2*QU(NK1)
2069 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
2070 qnewlq,qnewic,XLV1,XLV0)
2071 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
2072 TVDIFF = ABS(TU10-TVQU(NK1))
2073 IF(TVDIFF.LT.1.e-3)THEN
2078 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
2079 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
2080 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
2081 IF(EQFRC(NK1).EQ.1)THEN
2084 ELSEIF(EQFRC(NK1).EQ.0.)THEN
2089 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
2090 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
2092 CALL PROF5(EQFRC(NK1),EE2,UD2)
2096 ENDIF ! End of Entrain/Detrain IF BLOCK
2099 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
2100 ! VALUES IN THE LAYER...
2102 EE2 = AMAX1(EE2,0.5)
2104 UER(NK1)=0.5*REI*(EE1+EE2)
2105 UDR(NK1)=0.5*REI*(UD1+UD2)
2107 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
2108 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
2110 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
2112 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
2113 ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
2114 ! First, correct ABE calculation if needed...
2120 ! WRITE(98,1015)P0(NK1)/100.
2125 UPOLD=UMF(NK)-UDR(NK1)
2126 UPNEW=UPOLD+UER(NK1)
2128 DILFRC(NK1) = UPNEW/UPOLD
2130 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
2131 !...ICE IN THE DETRAINING UPDRAFT MASS...
2133 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
2134 DETIC(NK1)=QICE(NK1)*UDR(NK1)
2136 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
2137 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
2138 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
2139 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
2141 !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
2142 !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
2143 !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
2144 !...CURRENT MODEL LEVEL...
2146 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
2147 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
2149 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
2150 !!IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
2151 IF(NK1.LE.KCLDLAYER)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
2156 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
2157 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
2158 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BE
2159 ! THE LET AND CLOUD TOP...
2161 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOC
2162 ! FIRST BECOMES NEGATIVE...
2165 CLDHGT(LC)=Z0(LTOP)-ZLCL
2167 !...Instead of using the same minimum cloud height (for deep convection)
2168 !...everywhere, try specifying minimum cloud depth as a function of TLCL...
2172 IF(TLCL.GT.293.)THEN
2174 ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
2175 CHMIN = 2.E3 + 100.*(TLCL-273.)
2176 ELSEIF(TLCL.LT.273.)THEN
2181 !...If cloud top height is less than the specified minimum for deep
2182 !...convection, save value to consider this level as source for
2183 !...shallow convection, go back up to check next level...
2185 !...Try specifying minimum cloud depth as a function of TLCL...
2188 !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
2190 !... 1.) if there is no CAPE, or
2191 !... 2.) cloud top is at model level just above LCL, or
2192 !... 3.) cloud top is within updraft source layer, or
2193 !... 4.) cloud-top detrainment layer begins within
2194 !... updraft source layer.
2196 !!IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
2197 IF(LTOP.LE.KLCL .or. LTOP.LE.KCLDLAYER .or. LET+1.LE.KCLDLAYER)THEN ! No Convection Allowed
2209 ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
2214 !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
2217 EXIT usl ! Shallow Convection from this layer
2219 ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
2234 !!KSTART=MAX0(KPBL,KLCL)
2235 KSTART=MAX0(KCLDLAYER,KLCL)
2236 if (idiagee > 0) write(98,'(a,1p,2i5,2x,2i5)') &
2237 'kfcup let_old, let_new, klcl, ltop', let, kstart, klcl, ltop ! rce 11-may-2012
2241 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
2245 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
2246 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
2247 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
2252 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
2258 DUMFDP=UMF(LET)/DPTT
2260 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
2261 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
2265 !...entrainment is allowed at every level except for LTOP, so disallow
2266 !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
2267 !...so the the dilution factor due to entyrianment is not changed but
2268 !...the actual entrainment rate will change due due forced total
2269 !...detrainment in this layer...
2274 DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
2275 DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
2277 UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
2278 UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
2279 UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
2280 DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
2281 DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
2284 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
2285 PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
2286 PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
2287 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
2292 ! Initialize some arrays below cloud base and above cloud top...
2297 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
2298 UER(NK)=VMFLCL*DP(NK)/DPTHMX
2299 !!ELSEIF(NK.LE.KPBL)THEN
2300 ELSEIF(NK.LE.KCLDLAYER)THEN
2301 UER(NK)=VMFLCL*DP(NK)/DPTHMX
2302 UMF(NK)=UMF(NK-1)+UER(NK)
2307 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
2328 CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
2335 !...DEFINE VARIABLES ABOVE CLOUD TOP...
2350 !IF(NK.GT.LTOP1)THEN
2351 IF(NK.GE.LTOP1)THEN !BSINGH(11/12/2014): So that wu, qu and tu has a value for NK==LTOP1
2370 ! rce 11-may-2012 mods start -------------------------------------------
2371 ! calc droplet number (qndropbb)
2372 if ( flag_chem ) then
2375 if (nk1 == klcl) then
2376 ! calculate aerosol activatation at cloud base
2378 rho_act = p0(nk1)/(r*tu(nk1)*(1.+0.608*qu(nk1)))
2379 ! with cup, wlcl can be 0.0, so use wu(k+1) when wlcl is small
2381 if (wlcl < 0.1) w_act = max( w_act, wu(nk1) )
2383 ! effective w_act accounting for entrainment, from Barahona and Nenes (2007) eqn 14b
2385 ! w_act_effective = w_act * escale
2387 ! escale = 1 + (eBN/alphaBN) * [ (delHv*Mw /(Ru*T*T))*deltaT - deltaRH ]
2388 ! 1 + (eBN/alphaBN) * [ (delHv*ep2/(Ra*T*T))*deltaT - deltaRH ]
2390 ! eBN = entrainment rate = d[ln(updraft_mass_flux)]/dz
2391 ! alphaBN = [g*Mw *delHv/(cp*Ru*T*T)] - [g*Ma/(Ru*T)]
2392 ! = [g*ep2*delHv/(cp*Ra*T*T)] - [g /(Ra*T)]
2394 ! Mw, Ma = molecular weights of water and air ; ep2 = Mw/Ma
2395 ! delHv = latent heat of vaporization
2396 ! Ru = universal gas constant ; Ra = dry-air gas const
2397 ! deltaT = Tupdr - Tenv ; deltaRH = RHupdr - RHenv = 1 - RHenv
2398 tmpa = max( umf(nk), 1.0e-10 )
2399 tmpb = max( uer(nk1), 0.0 )
2400 tmpe = tmpb/(tmpa+0.5*tmpb)
2401 tmp_lv = xlv0 - xlv1*tk_act
2402 tmp_deltatkfact = tmp_lv*ep2/(r*tk_act*tk_act)
2403 tmp_alphabn = tmp_deltatkfact*g/cp - g/(r*tk_act)
2404 tmp_ebn = tmpe/dzq(nk1)
2405 tmp_deltatk = tk_act - t0(nk1)
2406 tmp_deltarh = 1.0 - q0(nk1)/qu(nk1)
2407 tmp_escale = 1.0 + (tmp_ebn/tmp_alphabn) * (tmp_deltatkfact*tmp_deltatk - tmp_deltarh)
2409 if (qndrop_cldbase_entrain_opt == 1) w_act_eff = w_act*tmp_escale
2410 w_act_eff = max( w_act_eff, w_act_min )
2413 if (idiagee > 0) then
2414 write(98,'(//a,8i5)') 'kfcup bb activate_cldbase_kfcup - i, j, nu, kcheck, ksrc1/2', &
2415 i, j, nu, kcheck(nu), lc, kcldlayer
2416 write(98,'( a,3i11 )') 'nk1, klcl, k ', nk1, klcl, k
2417 write(98,'( a,3i11 )') 'cldbase_entopt, incloud_entopt ', qndrop_cldbase_entrain_opt, qndrop_incloud_entrain_opt
2418 write(98,'( a,1p,8e11.3)') 'wlcl, wu(nk1), w_act, _eff, _min ', wlcl, wu(nk1), w_act, w_act_eff, w_act_min
2419 write(98,'( a,1p,8e11.3)') 'r, p, t, q, rho ', r, p0(nk1), tk_act, qu(nk1), rho_act
2420 write(98,'( a,1p,8e11.3)') 'g, r, cp, ep2, xlv0, xlv1, tmp_lv ', g, r, cp, ep2, xlv0, xlv1, tmp_lv
2421 write(98,'( a,1p,8e11.3)') 'tmpa/dx2, tmpb/dx2, tmpe ', tmpa/dxsq, tmpb/dxsq, tmpe
2422 write(98,'( a,1p,8e11.3)') 'ebn, dzq(nk1), dz... ', tmp_ebn, dzq(nk1), z_at_w1d(nk1+1)-z_at_w1d(nk1)
2423 write(98,'( a,1p,8e11.3)') 'deltarh, deltatk, deltatk*factor ', tmp_deltarh, tmp_deltatk, tmp_deltatk*tmp_deltatkfact
2424 write(98,'( a,1p,8e11.3)') 'escale, alphabn, deltatkfact ', tmp_escale, tmp_alphabn, tmp_deltatkfact
2426 call activate_cldbase_kfcup( idiagee, grid_id, ktau, &
2427 i, j, nk1, kts, kte, lc, kcldlayer, &
2428 num_chem, maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2429 ntype_aer, nsize_aer, ncomp_aer, &
2430 ai_phase, msectional, massptr_aer, numptr_aer, &
2431 dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer, &
2432 tk_act, rho_act, dp, w_act_eff, &
2433 chem1d, qndrop_tmp )
2434 qndrop_tmp = qndrop_tmp
2437 ! calculate dilution from entrainment
2438 ! umf(nk) is flux at bottom of layer nk1 ; uer(nk1) is entrainment (delta-umf) in layer nk1
2439 tmpa = max( umf(nk), 1.0e-10 )
2440 tmpb = max( uer(nk1), 0.0 )
2441 if (qndrop_incloud_entrain_opt == 1) then
2442 ! qndrop at center of layer nk1
2443 qndropbb(nk1) = qndrop_tmp*(tmpa/(tmpa+0.5*tmpb))
2444 ! qndrop at top of layer nk1
2445 qndrop_tmp = qndrop_tmp*(tmpa/(tmpa+tmpb))
2447 qndropbb(nk1) = qndrop_tmp
2449 if (idiagee > 0 .and. nk1 <= klcl+4) then
2450 write(98,'( a,i3,1p,8e11.3)') 'nk1, tmpa/dx2, tmpb/dx2, qndrop', nk1, tmpa/dxsq, tmpb/dxsq, qndropbb(nk1)
2453 if (idiagee > 0) write(98,'(a)')
2454 end if ! ( flag_chem ) then
2455 ! rce 11-may-2012 mods end ---------------------------------------------
2458 EMS(NK)=DP(NK)*DXSQ/G
2461 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCH
2463 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
2464 THTAU(NK)=TU(NK)*EXN(NK)
2465 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
2466 THTA0(NK)=T0(NK)*EXN(NK)
2467 DDILFRC(NK) = 1./DILFRC(NK)
2470 ! IF (XTIME.LT.10.)THEN
2471 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
2472 ! * TMIX-T00,PMIX,QMIX,ABE
2473 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
2477 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
2478 !...AND MIDTROPOSPHERE IS USED.
2480 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
2481 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
2482 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
2483 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
2484 !...for ETA model, DX is a function of location...
2485 ! TIMEC=DX(I,J)/VCONV
2488 TIMEC=AMAX1(1800.,TIMEC)
2489 TIMEC=AMIN1(3600.,TIMEC)
2490 !!IF(ISHALL.EQ.1)TIMEC=2400.
2491 IF(ISHALL.EQ.1)TIMEC=TIMEC_SHALL ! Reduced time constant, lkb 3/31/10
2495 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
2497 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
2502 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
2504 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
2505 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
2509 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
2511 CBH=(ZLCL-Z0(1))*3.281E-3
2515 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
2516 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
2518 IF(CBH.GT.25)RCBH=2.4
2520 PEFCBH=AMIN1(PEFCBH,.9)
2522 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
2524 PEFF=.5*(PEF+PEFCBH)
2525 PEFF2 = PEFF ! JSK MODS
2527 WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
2530 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
2531 !*****************************************************************
2533 ! COMPUTE DOWNDRAFT PROPERTIES *
2535 !*****************************************************************
2539 devap:IF(ISHALL.EQ.1)THEN
2541 DMF(1:KX)=0. ! rce 11-may-2012 - zero these out to avoid problems
2542 DER(1:KX)=0. ! with unit=98 diagnostic output
2550 !...start downdraft about 150 mb above cloud base...
2552 ! KSTART=MAX0(KPBL,KLCL)
2553 ! KSTART=KPBL ! Changed 7/23/99
2554 !!KSTART=KPBL+1 ! Changed 7/23/99
2555 KSTART=KCLDLAYER+1 ! Changed 7/23/99
2558 DPPP = P0(KSTART)-P0(NK)
2559 ! IF(DPPP.GT.200.E2)THEN
2560 IF(DPPP.GT.150.E2)THEN
2565 KLFS = MIN0(KLFS,LET-1)
2568 !...if LFS is not at least 50 mb above cloud base (implying that the
2569 !...level of equil temp, LET, is just above cloud base) do not allow a
2572 IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
2573 THETED(LFS) = THETEE(LFS)
2576 !...call tpmix2dd to find wet-bulb temp, qv...
2578 call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
2579 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
2581 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
2583 TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
2584 RDD=P0(LFS)/(R*TVD(LFS))
2589 RHBAR = RH(LFS)*DP(LFS)
2591 DO ND = LFS-1,KSTART,-1
2593 DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
2595 DMF(ND)=DMF(ND1)+DER(ND)
2596 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
2597 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
2599 RHBAR = RHBAR+RH(ND)*DP(ND)
2602 DMFFRC = 2.*(1.-RHBAR)
2604 !...Calculate melting effect
2605 !... first, compute total frozen precipitation generated...
2609 PPTMLT = PPTMLT+PPTICE(NK)
2612 !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
2613 !...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
2615 DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
2619 LDT = MIN0(LFS-1,KSTART-1)
2621 call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
2623 tz(kstart) = tz(kstart)-dtmelt
2624 ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
2625 QSS=0.622*ES/(P0(KSTART)-ES)
2626 THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))* &
2627 EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
2629 LDT = MIN0(LFS-1,KSTART-1) ! Determine the level to start at,
2630 ! KSTART is level of PBL
2633 THETED(ND) = THETED(KSTART)
2636 !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
2638 call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
2641 !...specify RH decrease of 20%/km in downdraft...
2643 RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
2645 !...adjust downdraft TEMP, Q to specified RH:
2648 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
2650 DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
2652 ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ)) ! Teten's equation to find Es
2653 QSRH=0.622*ES/(P0(ND)-ES) ! Find the sat. mixing ratio
2655 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
2656 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
2658 IF(QSRH.LT.QD(ND))THEN
2660 T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
2666 TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
2667 IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
2672 IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
2675 DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
2677 DMF(ND) = DMF(ND1)+DDR(ND)
2678 TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
2680 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
2686 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
2687 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
2689 d_mf: IF(TDER.LT.1.)THEN
2691 !3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
2709 DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
2711 IF(TDER*DDINC.GT.TRPPT)THEN
2716 DMF(NK)=DMF(NK)*DDINC
2717 DER(NK)=DER(NK)*DDINC
2718 DDR(NK)=DDR(NK)*DDINC
2724 write(98,*)'PRECIP EFFICIENCY =',PEFF
2729 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
2730 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
2731 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
2734 ! UMF(NK)=UMF(NK)*UPDINC
2735 ! UDR(NK)=UDR(NK)*UPDINC
2736 ! UER(NK)=UER(NK)*UPDINC
2737 ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
2738 ! PPTICE(NK)=PPTICE(NK)*UPDINC
2739 ! DETLQ(NK)=DETLQ(NK)*UPDINC
2740 ! DETIC(NK)=DETIC(NK)*UPDINC
2743 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
2773 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFL
2774 ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILAB
2775 ! IN THAT LAYER INITIALLY...
2780 !IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
2781 IF((UER(NK)-DER(NK)).GT.1.e-5)THEN
2782 AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
2783 !write(*,*) 'Larry... LMAX ', LMAX, LC, UER(NK), DER(NK)
2784 AINCMX=AMIN1(AINCMX,AINCM1)
2788 IF(AINCMX.LT.AINC)AINC=AINCMX
2790 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
2791 !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
2797 DETLQ2(NK)=DETLQ(NK)
2798 DETIC2(NK)=DETIC(NK)
2811 IF(ISHALL.EQ.1)THEN ! First for shallow convection
2813 ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
2814 ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
2815 ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
2817 !...find the maximum TKE value between LC and KLCL...
2821 ! DO 173 K = LC,KLCL
2823 ! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
2825 ! TKEMAX = AMIN1(TKEMAX,10.)
2826 ! TKEMAX = AMAX1(TKEMAX,5.)
2828 !c...3_24_99...DPMIN was changed for shallow convection so that it is the
2829 !c... the same as for deep convection (5.E3). Since this doubles
2830 !c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
2831 !c... lation of EVAC...
2832 !c EVAC = TKEMAX*0.1
2833 EVAC = 0.5*TKEMAX*0.1
2834 !!EVAC = 0.5*TKEMAX*0.1*freq
2835 ! AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
2836 !!AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
2837 !!AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
2838 AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC) * freq * 2.0 ! Use factor of two becuase only 1/2 of pdf would be expected to rise
2839 !!write(*,*) 'Larry ... old AINC ', AINC
2840 !!AINC = WLCL*freq*DXSQ*RHOLCL/(VMFLCL) ! This version uses mass flux from CuP
2846 UMF(NK)=UMF2(NK)*AINC
2847 DMF(NK)=DMF2(NK)*AINC
2848 DETLQ(NK)=DETLQ2(NK)*AINC
2849 DETIC(NK)=DETIC2(NK)*AINC
2850 UDR(NK)=UDR2(NK)*AINC
2851 UER(NK)=UER2(NK)*AINC
2852 DER(NK)=DER2(NK)*AINC
2853 DDR(NK)=DDR2(NK)*AINC
2855 ENDIF ! Otherwise for deep convection
2856 ! use iterative procedure to find mass fluxes...
2857 iter: DO NCOUNT=1,10
2859 !*****************************************************************
2861 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
2863 !*****************************************************************
2865 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
2866 !...SATISFY MASS CONTINUITY...
2870 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
2872 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
2873 ABSOMG = ABS(OMG(NK))
2874 ABSOMGTC = ABSOMG*TIMEC
2875 FRDP = 0.75*DP(NK-1)
2876 IF(ABSOMGTC.GT.FRDP)THEN
2885 NSTEP=NINT(TIMEC/DTT+1)
2886 DTIME=TIMEC/FLOAT(NSTEP)
2887 FXM(NK)=OMG(NK)*DXSQ/G
2890 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
2894 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
2895 !...SIGN OF OMEGA...
2904 IF(OMG(NK).LE.0.)THEN
2905 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
2906 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
2907 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
2908 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
2910 THFXOUT(NK)=FXM(NK)*THPA(NK)
2911 QFXOUT(NK)=FXM(NK)*QPA(NK)
2912 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
2913 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
2917 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
2920 THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
2921 THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
2923 QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)- &
2924 QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
2932 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORRO
2933 !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
2936 IF(QG(NK).LT.0.)THEN
2937 IF(NK.EQ.1)THEN ! JSK MODS
2938 ! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
2939 ! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
2940 CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
2946 TMA=QG(NK1)*EMS(NK1)
2947 TMB=QG(NK-1)*EMS(NK-1)
2948 TMM=(QG(NK)-1.E-9)*EMS(NK )
2949 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
2950 ACOEFF=BCOEFF*TMA/TMB
2954 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
2955 ! IF(ABS(QVDIFF).GT.1.)THEN
2956 ! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
2958 ! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
2959 ! 'VALUES IN KAIN-FRITSCH'
2963 QG(NK1)=TMA*EMSD(NK1)
2964 QG(NK-1)=TMB*EMSD(NK-1)
2967 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
2968 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
2969 ! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
2970 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
2971 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
2977 !...CONVERT THETA TO T...
2980 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
2981 TG(NK)=THTAG(NK)/EXN(NK)
2982 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
2985 ! write(*,*) 'Larry, exiting iter ',NCOUNT
2986 if (idiagee > 0) write(*,*) 'Larry, exiting iter - ncount,i,j',NCOUNT, I, J ! rce 11-may-2012
2988 ! write(*,*) 'Larry, exited, no more iter'
2991 !*******************************************************************
2993 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
2995 !*******************************************************************
2997 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
3003 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
3004 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
3009 TMIX=TMIX+DP(NK)*TG(NK)
3010 QMIX=QMIX+DP(NK)*QG(NK)
3014 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
3015 QSS=0.622*ES/(PMIX-ES)
3017 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
3021 CPM=CP*(1.+0.887*QMIX)
3022 DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
3023 DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
3029 EMIX=QMIX*PMIX/(0.622+QMIX)
3035 value=(indlu-1)*binc+astrt
3036 aintrp=(a1-value)/binc
3037 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
3038 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
3039 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
3040 TLCL=AMIN1(TLCL,TMIX)
3042 TVLCL=TLCL*(1.+0.608*QMIX)
3043 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
3046 IF(ZLCL.LE.Z0(NK))THEN
3051 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
3053 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
3055 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
3056 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
3057 TVEN=TENV*(1.+0.608*QENV)
3058 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
3059 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
3060 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
3062 !...COMPUTE ADJUSTED ABE(ABEG).
3067 THETEU(NK1) = THETEU(NK)
3069 call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
3071 TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
3074 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
3077 DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
3079 IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
3081 !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
3083 CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
3084 THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
3087 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
3088 !...THE PERIOD TIMEC...
3092 ! write(98,*)'TAU, I, J, =',NTSD,I,J
3093 ! WRITE(98,1060)FABE
3097 DABE=AMAX1(ABE-ABEG,0.1*ABE)
3099 IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
3100 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
3101 ! *GRID POINT; NO CONVECTION ALLOWED!'
3102 ! wig, 29-Aug-2006: Indicate no convection occurred.
3107 IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
3112 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
3121 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
3123 ! write(98,*)'TAU, I, J, =',NTSD,I,J
3124 ! WRITE(98,1055)FABE
3128 ! If there are shallow clouds, relax 90% requiremnt
3129 ! This code is not needed, exit out of shallow cu happens earlier
3130 !!IF(ISHALL .EQ. 1) THEN
3132 IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
3135 IF(NCOUNT.GT.10)THEN
3137 ! write(98,*)'TAU, I, J, =',NTSD,I,J
3138 ! WRITE(98,1060)FABE
3143 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTI
3144 !...MASS FLUX BY THE FACTOR AINC:
3149 IF(DABE.LT.1.e-4)THEN
3154 AINC=AINC*STAB*ABE/DABE
3157 ! AINC=AMIN1(AINCMX,AINC)
3158 AINC=AMIN1(AINCMX,AINC)
3159 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
3160 !...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
3161 IF(AINC.LT.0.05)then
3162 ! wig, 29-Aug-2006: Indicate no convection occurred.
3166 ! AINC=AMAX1(AINC,0.05) ! JSK MODS
3169 ! IF (XTIME.LT.10.)THEN
3170 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
3174 UMF(NK)=UMF2(NK)*AINC
3175 DMF(NK)=DMF2(NK)*AINC
3176 DETLQ(NK)=DETLQ2(NK)*AINC
3177 DETIC(NK)=DETIC2(NK)*AINC
3178 UDR(NK)=UDR2(NK)*AINC
3179 UER(NK)=UER2(NK)*AINC
3180 DER(NK)=DER2(NK)*AINC
3181 DDR(NK)=DDR2(NK)*AINC
3184 !...GO BACK UP FOR ANOTHER ITERATION...
3189 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
3191 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
3192 !...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
3194 ! Redistribute hydormeteors according to the final mass-flux values:
3197 FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
3206 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
3207 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
3211 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAY
3212 !...BASED ON THE SIGN OF OMEGA...
3225 IF(OMG(NK).LE.0.)THEN
3226 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
3227 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
3228 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
3229 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
3230 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
3231 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
3232 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
3233 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
3235 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
3236 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
3237 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
3238 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
3239 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
3240 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
3241 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
3242 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
3246 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
3249 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
3250 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
3251 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
3252 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
3262 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
3265 ! IF (XTIME.LT.10.)THEN
3266 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
3269 WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
3273 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
3277 ! if(I.eq.16 .and. J.eq.41)then
3278 ! IF(ISTOP.EQ.1)THEN
3280 ! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
3281 write(98,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
3282 TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
3283 write(98,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
3285 WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
3286 TMIX-T00,PMIX,QMIX,ABE
3287 WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
3289 WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
3290 write(98,*)'PRECIP EFFICIENCY =',PEFF
3291 WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
3294 WRITE(98,1070)' P ',' DP ',' DT K/D ',' DR K/D ', &
3295 ' OMG ',' DOMGDP ',' UMF ',' UER ', &
3296 ' UDR ',' DMF ',' DER ' ,' DDR ',&
3297 ' EMS ',' W0 ',' DETLQ ',' DETIC '
3298 write(98,*)'just before DO 300...'
3302 DTT=(TG(K)-T0(K))*86400./TIMEC
3304 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
3305 UDFRC=UDR(K)*TIMEC*EMSD(K)
3306 UEFRC=UER(K)*TIMEC*EMSD(K)
3307 DDFRC=DDR(K)*TIMEC*EMSD(K)
3308 DEFRC=-DER(K)*TIMEC*EMSD(K)
3309 WRITE(98,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
3310 UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
3311 W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
3315 ! rce 11-may-2012 mods start -------------------------------------------
3316 if (idiagee > 0) then
3317 write(98,'(/31x,3x,15a11)') 'umf/aeai', 'uer/aeai', 'umf/ae', 'uer/ae'
3318 do k = klcl-2, ltop+2
3321 write(98,'(31x,i3,1p,15e11.3)') k, umf(k)/(dxsq*ainc), uer(k)/(dxsq*ainc), umf(k)/dxsq, uer(k)/dxsq
3324 write(98,'(/a,1p,15i11 )') 'lc, kcldx, klcl, ksvaa, let, ltop', lc, kcldlayer, klcl, ksvaa, let, ltop
3325 write(98,'( a,1p,15e11.3)') 'dt, timec, dx, ae=dxsq, au0, ainc', dt, timec, dx, dxsq, au0, ainc
3326 write(98,'(a,1p,15e11.3 )') 'au0/ae, au0*ainc/ae ', au0/dxsq, au0*ainc/dxsq
3327 write(98,'(a,1p,15e11.3 )') 'vmflcl/ae, vmflcl*ainc/ae ', vmflcl/dxsq, vmflcl*ainc/dxsq
3328 write(98,'(a,1p,15e11.3 )') 'evac, freq, timec, tmp1 / 2 / 3 ', &
3329 evac, freq, timec, (dpthmx/g), (dpthmx/g)*(2.0*evac*freq), (vmflcl*ainc/dxsq)*timec
3330 write(98,'(a,1p,15e11.3 )') 'wlcl, wu(klcl), tpert, rpert ', wlcl, wu(klcl), th_perturb,r_perturb
3331 write(98,'( a,1p,15e11.3)') 'tmpc = (umf/ae)/(wu*rho) tmpd = umf/(wu*rho*au0*ainc)'
3332 write(98,'(3x,15a11)') 'p0', 'dp', 'omg/g', 'umf/ae', 'del-umf', 'uer-udr', 'uer/ae', 'udr/ae', 'wu', 'tmpc', 'tmpd', 'ems'
3333 do k = ltop+2, 1, -1
3335 tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0 ; tmpd = 0.0
3336 if (k > 1 .and. k < ltop) tmpa = (umf(k)-umf(k-1))/dxsq
3337 if (k == ltop) tmpa = (0.0 -umf(k-1))/dxsq
3338 tmpb = (uer(k)-udr(k))/dxsq
3339 if (wu(k) > 1.0e-3) tmpc = umf(k)/(wu(k)*rhoe(k)*dxsq)
3340 if (wu(k) > 1.0e-3) tmpd = umf(k)/(wu(k)*rhoe(k)*au0*ainc)
3341 write(98,'(i3,1p,15e11.3)') k, p0(k), dp(k), omg(k)/g, umf(k)/dxsq, tmpa, tmpb, uer(k)/dxsq, udr(k)/dxsq, wu(k), tmpc, tmpd, ems(k)
3344 write(98,'(/3x,15a11)') 't0', 'p0', 'dp', 'q0', 'qg', 'qu', 'qliq', 'qlg', 'qice', 'qig', 'qndropbb'
3346 write(98,'(i3,f11.2,1p,15e11.3)') k, t0(k)-t00, p0(k), dp(k), q0(k), qg(k), qu(k), qliq(k), qlg(k), qice(k), qig(k), &
3351 ! rce 11-may-2012 mods end ---------------------------------------------
3353 WRITE(98,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0', &
3355 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
3360 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
3362 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
3363 IF(T0(K).LT.T00)THEN
3364 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
3366 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
3368 QGS=ES*0.622/(P0(K)-ES)
3371 WRITE(98,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
3372 TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
3373 1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
3374 QSG(K)*1000.,RH0,RHG
3377 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
3378 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
3380 ! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
3382 ! IF(ISHALL.NE.1)THEN
3383 ! write(98,4421)i,j,iyr,imo,idy,ihr,imn
3384 ! write(98)i,j,iyr,imo,idy,ihr,imn,kl
3388 write(98,'(8a11)') 'p0', 't0', 'q0', 'u0', 'v0', 'w0avg1d', 'dp', 'tke' ! rce 11-may-2012
3391 write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
3392 u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
3393 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
3394 ! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
3395 ! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
3398 CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
3403 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
3404 RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ ! PPT FB MODS
3405 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
3406 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
3408 IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
3410 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
3412 ! EVALUATE MOISTURE BUDGET...
3420 QINIT=QINIT+Q0(NK)*EMS(NK)
3421 QFNL=QFNL+QG(NK)*EMS(NK)
3422 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
3424 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) ! PPT FB MODS
3425 ! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
3426 ERR2=(QFNL-QINIT)*100./QINIT
3427 IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
3428 IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
3429 ! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
3430 ! WRITE(99,1110)QINIT,QFNL,ERR2
3437 ! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
3438 ! u0(k),v0(k),W0AVG1D(K),dp(k)
3439 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
3440 ! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
3441 ! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
3442 WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
3443 U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
3450 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
3452 IF(PPTFLX.GT.0.)THEN
3453 RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
3458 WRITE(98,1120)RELERR
3459 WRITE(98,*)'TDER, CPR, TRPPT =', &
3460 TDER,CPR*AINC,TRPPT*AINC
3463 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
3465 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
3466 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
3468 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
3469 NCA(I,J)=REAL(NIC)*DT !byang
3471 TIMEC = TIMEC_SHALL !! Changed to match other location where TIMEC is set lkb 10/31/10
3473 NCA(I,J) = NINT(TIMEC_SHALL/DT)*DT ! add 01/11/2012
3474 ! NCA(I,J) = NTST*DT !byang
3478 ! IF(IMOIST(INEST).NE.2)THEN
3480 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
3481 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
3482 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
3483 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
3486 ! RLC=XLV0-XLV1*TG(K)
3487 ! RLS=XLS0-XLS1*TG(K)
3488 ! CPM=CP*(1.+0.887*QG(K))
3489 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
3490 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
3497 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
3501 CPM=CP*(1.+0.887*QG(K))
3502 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
3503 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
3505 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
3507 ELSEIF(.NOT. F_QS)THEN
3509 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
3510 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
3512 CPM=CP*(1.+0.887*QG(K))
3514 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
3516 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
3518 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
3520 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
3524 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
3525 !...OF HYDROMETEORS DIRECTLY...
3527 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
3528 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
3529 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
3531 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
3533 DQSDT(K)=DQSDT(K)+(QIG(K)-QI0(K))/TIMEC
3536 ! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
3537 CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS MICROPHYSICS CHOICE IS NOT ALLOWED' )
3539 DTDT(K)=(TG(K)-T0(K))/TIMEC
3540 DQDT(K)=(QG(K)-Q0(K))/TIMEC
3543 ! PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ !LD add PRATEC 21-April-2011
3544 ! RAINCV(I,J)=DT*PRATEC(I,J) !LD add PRATEC 21-April-2011
3546 RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ ! PPT FB MODS
3548 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
3549 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
3551 909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
3552 ! write (98,909)I,J,RNC
3553 ! write (6,909)I,J,RNC
3554 ! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
3557 1000 FORMAT(' ',10A8)
3558 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
3559 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
3560 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
3561 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
3562 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
3563 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
3565 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
3566 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
3568 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
3570 !1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
3571 ! ', NO MORE MASS FLUX IS ALLOWED!')
3572 !1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
3573 ! &DEGREE OF STABILIZATION! FABE= ',F6.4)
3575 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
3576 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
3577 2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
3578 1085 FORMAT (A3,16A7,2A8)
3579 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
3580 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
3581 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
3582 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
3583 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
3584 ' TOTAL WATER CHANGE =',F8.2,'%')
3585 ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
3586 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
3588 !-----------------------------------------------------------------------
3589 !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
3590 !-----------------------------------------------------------------------
3592 IF (ISHALL<2) THEN ! add LKB 12/23/2011 01/11/2012 only define cloud base
3593 CUTOP(I,J)=REAL(LTOP) ! if there are clouds
3594 CUBOT(I,J)=REAL(LCL)
3596 ! rce 11-may-2012 mods start -------------------------------------------
3597 updfra = au0*ainc/dxsq
3602 qndrop1d(:) = qndropbb(:)
3604 ! umf(k) and umfout(k) are at top of layer k
3605 umfout(kts:ltop-1) = max( 0.0, umf(kts:ltop-1)/dxsq )
3606 uerout(kts:ltop) = max( 0.0, uer(kts:ltop)/dxsq )
3607 udrout(kts:ltop) = max( 0.0, udr(kts:ltop)/dxsq )
3608 ! dmf(k) is at bottom of layer k; ! dmfout(k) is at top of layer k [like umf(k) and umfout(k)]
3609 dmfout(kts:ltop-1) = min( 0.0, dmf(kts+1:ltop)/dxsq )
3610 ! der(k) is negative; derout(k) is positive
3611 derout(kts:ltop) = max( 0.0, -der(kts:ltop)/dxsq )
3612 ! ddr(k) is positive so no change needed
3613 ddrout(kts:ltop) = max( 0.0, ddr(kts:ltop)/dxsq )
3615 if ( idiagee > 0 .and. ((ishall == 0) .or. (ishall == 1)) ) then
3616 write(98,'(/a,1p,15i11 )') 'lc, kcldx, klcl, ksvaa, let, ltop', lc, kcldlayer, klcl, ksvaa, let, ltop
3617 write(98,'( a,1p,15e11.3)') 'dt, timec, dx, ae=dxsq, au0, ainc', dt, timec, dx, dxsq, au0, ainc
3618 write(98,'(a,1p,15e11.3 )') 'au0/ae, au0*ainc/ae ', au0/dxsq, au0*ainc/dxsq
3619 write(98,'(a,1p,15e11.3 )') 'wlcl, wu(klcl), tpert, rpert ', wlcl, wu(klcl), th_perturb,r_perturb
3620 tmpa = 0.0 ; tmpb = 0.0
3622 tmpa = tmpa + uerout(k)
3623 if (k >= klcl) tmpb = tmpb + dp(k)
3625 write(98,'(a,1p,15e11.3 )') 'tmpu, ...*tau, tmpv, ...*area/g ', tmpa, tmpa*dt*ntst, tmpb, (tmpb/g)*(au0*ainc/dxsq)
3626 write(98,'(3x,15a11)') 'p0', 'dp', 'omg/g', 'umfout', 'del-umf', 'uer-udr', 'uerout', 'udrout', &
3627 'qc1d', 'qi1d', 'f_qc2qi', 'f_qc2pr', 'f_qi2pr'
3628 do k = ltop+2, 1, -1
3630 tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0 ; tmpd = 0.0
3631 if (k > 1 ) tmpa = umfout(k)-umfout(k-1)
3632 if (k == 1) tmpa = umfout(k)
3633 tmpb = uerout(k)-udrout(k)
3634 write(98,'(i3,1p,15e11.3)') k, p0(k), dp(k), omg(k)/g, umfout(k), tmpa, tmpb, uerout(k), udrout(k), &
3635 qc1d(k), qi1d(k), fcvt_qc_to_qi(k), fcvt_qc_to_pr(k), fcvt_qi_to_pr(k)
3637 write(98,'(3x,15a11)') 'p0', 'dp', ' ', 'dmfout', 'del-dmf', 'der-ddr', 'derout', 'ddrout'
3638 do k = ltop+2, 1, -1
3640 tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0 ; tmpd = 0.0
3641 if (k > 1 ) tmpa = dmfout(k)-dmfout(k-1)
3642 if (k == 1) tmpa = dmfout(k)
3643 tmpb = derout(k)-ddrout(k)
3644 write(98,'(i3,1p,15e11.3)') k, p0(k), dp(k), 0.0, dmfout(k), tmpa, tmpb, derout(k), ddrout(k)
3646 end if ! ( idiagee > 0 .and. ((ishall == 0) .or. (ishall == 1)) ) then
3647 ! rce 11-may-2012 mods end ---------------------------------------------
3650 !-----------------------------------------------------------------------
3651 ! begin: wig, 21-Feb-2008
3652 ! Only allow shallow-Cu to occur if the cloud base is within 500 m of
3653 ! the top of the PBL. This prevents us from getting too many clouds
3654 ! in the mid-troposphere.
3655 if( ishall==1 .and. (z_at_w1d(lcl)-pblh) > 500. ) ishall = 2
3656 ! end: wig, 21-Feb-2008
3658 END SUBROUTINE KF_cup_PARA
3660 !********************************************************************
3661 ! ***********************************************************************
3663 SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
3665 ! Lookup table variables:
3666 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
3667 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3668 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3669 ! REAL, SAVE, DIMENSION(1:200) :: ALU
3670 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3671 ! End of Lookup table variables:
3672 !-----------------------------------------------------------------------
3674 !-----------------------------------------------------------------------
3675 REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
3676 REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
3677 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
3678 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
3679 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
3680 INTEGER :: IPTB,ITHTB
3681 !-----------------------------------------------------------------------
3683 !c******** LOOKUP TABLE VARIABLES... ****************************
3684 ! parameter(kfnt=250,kfnp=220)
3686 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
3687 ! * alu(200),rdpr,rdthk,plutop
3688 !C***************************************************************
3690 !c***********************************************************************
3691 !c scaling pressure and tt table index
3692 !c***********************************************************************
3694 ! plutop = model top pressure
3695 ! p = pressure level
3696 ! rdpr = a pressure (or 1/pressure) increment
3697 ! tp = a number of levels (a pressure difference divided by the increment
3703 !***********************************************************************
3704 ! base and scaling factor for the
3705 !***********************************************************************
3707 ! scaling the and tt table index
3708 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
3709 tth=(thes-bth)*rdthk
3712 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
3713 write(98,*)'**** OUT OF BOUNDS *********'
3717 t00=ttab(ithtb ,iptb )
3718 t10=ttab(ithtb+1,iptb )
3719 t01=ttab(ithtb ,iptb+1)
3720 t11=ttab(ithtb+1,iptb+1)
3722 q00=qstab(ithtb ,iptb )
3723 q10=qstab(ithtb+1,iptb )
3724 q01=qstab(ithtb ,iptb+1)
3725 q11=qstab(ithtb+1,iptb+1)
3727 !***********************************************************************
3728 ! parcel temperature
3729 !***********************************************************************
3731 temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
3733 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
3741 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
3742 ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
3747 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
3748 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
3749 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
3750 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
3751 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
3753 !...subsaturated values only occur in calculations involving various mixtures of
3754 !...updraft and environmental air for estimation of entrainment and detrainment.
3755 !...For these purposes, assume that reasonable estimates can be given using
3756 !...liquid water saturation calculations only - i.e., ignore the effect of the
3757 !...ice phase in this process only...will not affect conservative properties...
3760 qliq=qliq-dq*qliq/(qtot+1.e-10)
3761 qice=qice-dq*qice/(qtot+1.e-10)
3765 CPP=1004.5*(1.+0.89*QU)
3766 IF(QTOT.LT.1.E-10)THEN
3768 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
3769 TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
3772 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
3773 ! THE TEMPERATURE IS GIVEN BY:
3775 TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
3787 END SUBROUTINE TPMIX2
3788 !******************************************************************************
3789 SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
3790 !-----------------------------------------------------------------------
3792 !-----------------------------------------------------------------------
3793 REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
3794 REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
3795 REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
3796 !-----------------------------------------------------------------------
3798 !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
3799 !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
3800 !...TTFRZ TO TBFRZ...
3801 !...FOR COLDER TERMPERATURES, FREEZE ALL LIQUID WATER...
3802 !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
3803 !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
3805 RLC=2.5E6-2369.276*(TU-273.16)
3806 RLS=2833922.-259.532*(TU-273.16)
3808 CPP=1004.5*(1.+0.89*QU)
3810 ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
3811 ! FOR SATURATION VAPOR PRESSURE...
3813 A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
3814 DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
3817 ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
3818 QS = ES*0.622/(P-ES)
3820 !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
3821 !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
3822 !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
3823 !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
3824 !...TEMPERATURE TO THE SATURATION VALUE...
3829 PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
3830 THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
3832 END SUBROUTINE DTFRZNEW
3833 ! --------------------------------------------------------------------------------
3835 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
3836 QNEWIC,QLQOUT,QICOUT,G)
3838 !-----------------------------------------------------------------------
3840 !-----------------------------------------------------------------------
3841 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
3842 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
3843 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
3844 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
3845 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
3847 REAL, INTENT(IN ) :: G
3848 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
3849 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
3850 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
3853 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
3854 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
3855 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
3856 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
3857 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
3861 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
3862 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
3865 QEST=0.5*(QTOT+QNEW)
3866 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
3868 WAVG=0.5*(SQRT(WTW)+SQRT(G1))
3869 CONV=RATE*DZ/max(WAVG,1e-7) !wig, 12-Sep-2006: added div by 0 check
3871 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
3872 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
3873 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
3874 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
3876 RATIO3=QNEWLQ/(QNEW+1.E-8)
3880 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
3881 QTOT=QTOT*EXP(-CONV)
3883 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
3884 ! PARCEL AT THIS LEVEL...
3888 QICOUT=(1.-RATIO4)*DQ
3890 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
3891 ! LATE VERTICAL VELOCITY
3893 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
3894 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
3895 IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
3897 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
3898 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
3900 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
3901 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
3905 END SUBROUTINE CONDLOAD
3907 ! ----------------------------------------------------------------------
3908 SUBROUTINE PROF5(EQ,EE,UD)
3910 !***********************************************************************
3911 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
3912 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3913 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
3914 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
3915 ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
3916 ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
3917 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
3920 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3921 !-----------------------------------------------------------------------
3923 !-----------------------------------------------------------------------
3924 REAL, INTENT(IN ) :: EQ
3925 REAL, INTENT(INOUT) :: EE,UD
3926 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
3928 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
3929 0.9372980,0.33267,0.166666667,0.202765151/
3936 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
3937 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
3939 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
3940 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
3943 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
3944 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
3950 END SUBROUTINE PROF5
3952 ! ------------------------------------------------------------------------
3953 SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
3955 ! Lookup table variables:
3956 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
3957 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3958 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3959 ! REAL, SAVE, DIMENSION(1:200) :: ALU
3960 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3961 ! End of Lookup table variables:
3962 !-----------------------------------------------------------------------
3964 !-----------------------------------------------------------------------
3965 REAL, INTENT(IN ) :: P,THES
3966 REAL, INTENT(INOUT) :: TS,QS
3967 INTEGER, INTENT(IN ) :: i,j ! avail for debugging
3968 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
3969 INTEGER :: IPTB,ITHTB
3970 CHARACTER*256 :: MESS
3971 !-----------------------------------------------------------------------
3974 !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
3975 ! parameter(kfnt=250,kfnp=220)
3977 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
3978 ! alu(200),rdpr,rdthk,plutop
3979 !***************************************************************
3981 !***********************************************************************
3982 ! scaling pressure and tt table index
3983 !***********************************************************************
3989 !***********************************************************************
3990 ! base and scaling factor for the
3991 !***********************************************************************
3993 ! scaling the and tt table index
3994 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
3995 tth=(thes-bth)*rdthk
3999 t00=ttab(ithtb ,iptb )
4000 t10=ttab(ithtb+1,iptb )
4001 t01=ttab(ithtb ,iptb+1)
4002 t11=ttab(ithtb+1,iptb+1)
4004 q00=qstab(ithtb ,iptb )
4005 q10=qstab(ithtb+1,iptb )
4006 q01=qstab(ithtb ,iptb+1)
4007 q11=qstab(ithtb+1,iptb+1)
4009 !***********************************************************************
4010 ! parcel temperature and saturation mixing ratio
4011 !***********************************************************************
4013 ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
4015 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
4017 END SUBROUTINE TPMIX2DD
4019 ! -----------------------------------------------------------------------
4020 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
4022 !-----------------------------------------------------------------------
4024 !-----------------------------------------------------------------------
4025 REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
4026 REAL, INTENT(INOUT) :: THT1
4027 REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
4028 T00,P00,C1,C2,C3,C4,C5
4030 !-----------------------------------------------------------------------
4031 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
4034 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
4036 ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
4039 ! TLOG=ALOG(EE/ALIQ)
4040 ! ...calculate LOG term using lookup table...
4047 value=(indlu-1)*ainc+astrt
4048 aintrp=(a1-value)/ainc
4049 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
4051 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
4052 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
4053 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
4054 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
4056 END SUBROUTINE ENVIRTHT
4057 ! ***********************************************************************
4058 !====================================================================
4059 SUBROUTINE kf_cup_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
4060 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
4061 SVP1,SVP2,SVP3,SVPT0, &
4062 cupflag,cldfra_cup,cldfratend_cup, & !CuP, wig 18-Sep-2006
4063 shall, & !CuP, wig 18-Sep-2006
4064 tcloud_cup, & !CuP, rce 10-may-2012
4065 P_FIRST_SCALAR,restart,allowed_to_read, &
4066 ids, ide, jds, jde, kds, kde, &
4067 ims, ime, jms, jme, kms, kme, &
4068 its, ite, jts, jte, kts, kte )
4069 !--------------------------------------------------------------------
4071 !--------------------------------------------------------------------
4072 LOGICAL , INTENT(IN) :: restart,allowed_to_read
4073 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
4074 ims, ime, jms, jme, kms, kme, &
4075 its, ite, jts, jte, kts, kte
4076 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
4078 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4086 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4088 cldfra_cup, & !CuP, wig 18-Sep-2006
4089 cldfratend_cup !CuP, wig 18-Sep-2006
4091 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA, &
4092 shall, & !CuP, wig 19-Sep-2006
4093 tcloud_cup !CuP, rce 10-may-2012
4095 LOGICAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT):: cupflag !CuP, wig 9-Oct-2006
4097 INTEGER :: i, j, k, itf, jtf, ktf
4098 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
4104 IF(.not.restart)THEN
4113 cldfra_cup(i,k,j) = 0. !CuP, wig 18-Sep-2006
4114 cldfratend_cup(i,k,j) = 0. !CuP, wig 18-Sep-2006
4119 IF (P_QI .ge. P_FIRST_SCALAR) THEN
4129 IF (P_QS .ge. P_FIRST_SCALAR) THEN
4142 shall(i,j) = 2. !Indicate no convection at 1st time step. CuP, wig 18-Sep-2006
4143 cupflag(i,j) = .false. !CuP, wig 9-Oct-2006
4144 tcloud_cup(i,j) = 0.0 !CuP, rce 10-may-2012
4158 CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
4160 END SUBROUTINE kf_cup_init
4162 !-------------------------------------------------------
4164 subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
4166 ! This subroutine is a lookup table.
4167 ! Given a series of series of saturation equivalent potential
4168 ! temperatures, the temperature is calculated.
4170 !--------------------------------------------------------------------
4172 !--------------------------------------------------------------------
4173 ! Lookup table variables
4174 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
4175 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
4176 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
4177 ! REAL, SAVE, DIMENSION(1:200) :: ALU
4178 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
4179 ! End of Lookup table variables
4181 INTEGER :: KP,IT,ITCNT,I
4182 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
4183 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
4185 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
4186 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
4187 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
4189 ! equivalent potential temperature increment
4191 ! minimum starting temp
4193 ! tolerance for accuracy of temperature
4195 ! top pressure (pascals)
4197 ! bottom pressure (pascals)
4206 ! compute parameters
4208 ! 1._over_(sat. equiv. theta increment)
4210 ! pressure increment
4212 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
4213 ! dpr=(pbot-plutop)/REAL(kfnp-1)
4214 ! 1._over_(pressure increment)
4216 ! compute the spread of thes
4217 ! thespd=dth*(kfnt-1)
4219 ! calculate the starting sat. equiv. theta
4225 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
4227 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
4228 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
4232 ! compute temperatures for each sat. equiv. potential temp.
4239 ! define sat. equiv. pot. temp.
4241 ! iterate to find temperature
4242 ! find initial guess
4248 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
4250 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
4251 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
4259 es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
4261 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
4262 thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
4264 if(abs(f1).lt.toler)then
4268 dt=f1*(t1-t0)/(f1-f0)
4278 ! lookup table for tlog(emix/aliq)
4280 ! set up intial values for lookup tables
4291 END SUBROUTINE KF_LUTAB
4294 !--------------------------------------------------------------------
4295 ! Calculates the cloud fraction tendency.
4297 SUBROUTINE cupCloudFraction(qlg, qig, qv1d, t1d, z1d, p1d, &
4298 kcubot, kcutop, ishall, wStar, wParcel, pblh, dt, activeFrac, &
4299 cldfra, cldfraTend, &
4300 taucloud, tActive, tstar, lnterms, lnint, &
4301 kts, kte, mfup_cup) ! add mfup_cup LD 06 29 2012
4304 use module_model_constants, only: r_v, xls0, xls1, xlv0, xlv1
4308 integer, intent(in) :: kts, kte
4309 integer, intent(in) :: ishall ! Flag for cloud type (0=deep, 1=shallow, 2=none)
4311 integer, intent(in) :: kcubot, kcutop ! Indices of cloud top and bottom
4312 real, intent(in) :: wStar, pblh ! Deardorff velocity scale and mixed-layer depth
4313 real, intent(in) :: wParcel ! Vertical velocity of parcel
4314 real, intent(in) :: activeFrac ! Active cloud fraction, determined from cloud model
4315 real, intent(in) :: dt ! Time step used to find cloud fraction
4317 real, dimension(kts:kte), intent(in) :: qlg ! Cloud liquid water
4318 real, dimension(kts:kte), intent(in) :: qig ! Cloud ice
4319 real, dimension(kts:kte), intent(in) :: t1d ! Environment temperature
4320 real, dimension(kts:kte), intent(in) :: qv1d ! Environmental mixing ratio
4321 real, dimension(kts:kte), intent(in) :: z1d ! Height array on cell middles
4322 real, dimension(kts:kte), intent(in) :: p1d ! Pressure array
4324 real, dimension(kts:kte), intent(inout) :: cldfra ! Cloud fraction
4325 real, dimension(kts:kte), intent(in) :: mfup_cup ! LD 06 29 2012
4326 real, dimension(kts:kte), intent(out) :: cldfraTend ! Cloud fraction tendency
4330 integer :: k, kp1 ,kcutop_p1 !BSINGH - Added kcutop_p1
4333 real,intent(out) :: tauCloud ! Cloud time scale ~can make local after testing
4334 real,intent(out) :: tActive ! Cloud time scale ~can make local after testing
4335 !!! real,intent(out) :: wParcel ! Cloud velocity scale ~can make local after testing
4336 real,intent(out) :: tStar ! Boundary-layer time scale ~can make local after testing
4337 !!! real,intent(out) :: activeFrac ! Fraction of PDF that forms clouds
4338 real :: ice_term, liquid_term ! Terms inside of log for gamma
4339 real, dimension(kts:kte),intent(out) :: lnTerms ! Combined log terms to be integrated ~can make local after testing
4340 real,intent(out) :: lnInt ! Integrated log terms for gamma ~can make local after testing
4341 real :: intQC ! Integrated cloud water add 2010/01/17
4342 real, dimension(kts:kte) :: satDef ! Saturation deficit add 2010/01/17
4343 real :: intSatDef ! Integrated saturation deficit aa 2010/01/17
4344 real :: deltaZ ! Height diff. between cell centers
4345 real :: deltaRsInt ! Integrated delta rs
4346 real :: deltaRsTop, deltaRsBot ! Value at deltaRs at the top an bottom of the layer
4347 real :: TEnvTop, TEnvBot ! Env. temperature at top and bottom of the layer
4348 real :: rs, rsi ! Saturation mixing ratios w.r.t liquid and ice
4349 real :: cp , Ls, Lv ! Thermodynamic related "constants"
4351 if( ishall==2 ) then
4352 ! If no convection, then zero out the cloud fraction...
4355 else if( ishall==0 ) then
4356 ! If deep convection formed, then set the cloud fraction to 1.
4359 ! cldfra(kcubot:kcutop) = 1. !!LD
4360 !! print(UMF(?)) unit??
4362 cldfra(k) = max(0.,min(0.1*log(1.+675.*mfup_cup(k)),1.)) !! LD 06 29 2012 :: .1/675 adjustable parameter
4364 ! print*,"mfup_cup(kcubot)=",mfup_cup(kcubot)
4366 tStar = pblh / wStar ! rce 11-may-2012
4368 else if( ishall==1 ) then
4370 ! Shallow convection occurred so we need to be more detailed...
4372 tStar = pblh / wStar ! Find tStar based on mixed-layer depth
4374 ! Integrate the log terms for the cloud time scale over the depth
4375 ! of the cloud and take into account both liquid and ice as
4376 ! separate terms. Do not allow super saturation, and at the same
4377 ! time, preclude divide by zeros by limiting the (rs-r)'s to
4381 !!The determination of the cloud time scales around line 3523.
4382 !!modified the do loop that computes the integrated cloud water and the saturation deficit.
4383 !!code starts at line 3541 and continues through 3560
4385 !!do k=kcubot,kcutop
4386 !!cp = findCp(qv1d(k))
4387 !!rs = findRs(t1d(k), p1d(k))
4388 !!Lv = xlv0 - xlv1*t1d(k)
4389 !!gamma = eps*(Lv**2)*rs / (cp*r_v*t1d(k)**2)
4390 !!liquid_term = (1.+gamma)*qlg(k) / max(rs - qv1d(k),1e-20)
4392 !!Ls = xls0 - xls1*t1d(k)
4393 !!rsi = findRsi(t1d(k), p1d(k))
4394 !!gamma = eps*(Ls**2)*rsi / (cp*r_v*t1d(k)**2)
4395 !!ice_term = (1.+gamma)*qig(k) / max(rsi - qv1d(k),1e-20)
4397 !!lnTerms(k) = 1. + liquid_term !~tmp + ice_term
4400 !lnInt = 0.! add 2011/01/16 start
4405 !BSINGH - Added do-loop to compute 'satDef' before it is being used in the next do-loop
4406 !BSINGH - This loop should go to (kcutop+1) as we are trying to access satDef(k+1) in the next do-loop
4407 kcutop_p1 = min(kcutop + 1,kte)
4408 do k = kcubot, kcutop_p1
4409 rs = findRs(t1d(k), p1d(k))
4410 satDef(k) = max(rs - qv1d(k), 1.0e-20)
4415 kp1 = min(k+1,kte-1)
4416 deltaZ = z1d(kp1) - z1d(k) ! Find the interval
4417 zsum = zsum + deltaz
4418 !!lnInt = lnInt + 0.5*(lnTerms(k) + lnTerms(kp1))*deltaZ
4419 rs = findRs(t1d(k), p1d(k))
4420 satDef(k) = max(rs - qv1d(k), 1e-20)
4421 intQC = 0.5*(qlg(k) + qlg(kp1)) * deltaz + intQC
4422 ! print *, 'Values within cupCloudFraction', intSatDef, satDef(k),satDef(kp1),deltaz,k,kp1
4423 !print *, 'Values within cupCloudFraction',rs,qv1d(k),qlg(k),qlg(kp1),k,kp1
4424 intSatDef = 0.5*(satDef(k) + satDef(kp1)) * deltaz + intSatDef
4426 !!lnInt = lnInt/zsum !Turn the integral into an average
4427 cp = findCp(qv1d(kcubot)) ! Use the thermodynamic properties at cloud base for defining gamma
4428 rs = findRs(t1d(kcubot), p1d(kcubot))
4429 Lv = xlv0 - xlv1*t1d(kcubot)
4430 gamma = (Lv**2)*rs / (cp*r_v*t1d(kcubot)**2)
4431 lnInt = log(1.0 + (1.0 + gamma) * intQC / intSatDef)
4432 lnInt = max(lnInt, 1.0) ! Set the value of lnInt to be 1 or greater
4433 ! add 2011/01/16 end
4434 ! Find the time scale of the cloud lifetime, tauCloud, and the time
4435 ! scale of the cloud formation, tActive...
4436 !!tauCloud = min(tStar*lnInt, 3600.) !Set a max taucld of 60 min.
4437 tauCloud = min(tStar*lnInt, 1800.) !Set a max taucld of 60 min.
4438 if(wParcel .gt. 0) then
4439 tActive = z1d(kcutop)/wParcel
4441 tActive = z1d(kcutop) / wStar
4443 !!!! tActive or tactive matter? Dec-15-2010-LP
4444 !!$ ! Now, find the cloud fraction tendency. Above and below the cloud,
4446 !!$ cldfraTend(kts:max(kcubot-1,kts)) = 0.
4447 !!$ cldfraTend(min(kcutop+1,kte):kte) = 0.
4448 !!$ do k=kcubot,kcutop
4449 !!$ cldfraTend(k) = dt*(activeFrac/tActive - cldfra(k)/tauCloud)
4452 ! Now, get a steady-state cloud fraction and restrict it to the
4456 cldfra(k) = activeFrac*tauCloud/tActive
4457 cldfra(k) = max(cldfra(k), 0.01) ! LKB 9/9/09 Changed from 0 to be 0.1
4458 cldfra(k) = min(cldfra(k), 1.)
4462 !This should never happen!
4463 call wrf_error_fatal("Bad ishall value in kfcup.")
4466 END SUBROUTINE cupCloudFraction
4469 !------------------------------------------------------------------------
4470 SUBROUTINE cup_jfd(slopeSfc, slopeEZ, sigmaSfc, sigmaEZ, &
4471 numBins, thBinSize, rBinSize, th_perturb, r_perturb, jfd )
4473 USE module_model_constants, only: pi2
4477 integer, intent(in) :: numBins
4478 real, intent(in) :: thBinSize, rBinSize
4479 real, intent(inout) :: slopeSfc, slopeEZ, sigmaSfc, sigmaEZ
4480 real, dimension(numBins), intent(out) :: r_perturb, th_perturb
4481 real, dimension(numBins,numBins), intent(out) :: jfd
4485 integer :: centerBin, i, j
4486 real :: bigcheck, c, constants, cterm, dslope, jacobian, jfdCheckSum, m, mterm
4487 character(len=150) :: message
4489 ! Limit the allowable values of the slopes and sigmas ~get the right values for caps
4491 ! slopeSfc = sign( min( abs(slopeSfc), 2e6 ), slopeSfc)
4492 ! slopeEZ = sign( min( abs(slopeEZ), 2e6 ), slopeEZ)
4493 !~ sigmaSfc = max( abs(sigmaSfc), rBinSize ) ! <-- This one is the only one that really limited anything. It was only giving the value rBinSize.
4494 !~ sigmaEZ = max( abs(sigmaEZ), rBinSize )
4496 !!$!~wig begin: testing due to overflow of jfd calc 13-dec-2006
4497 !!$if( abs(slopesfc) < 1e-14 ) print*,"small slopesfc =",slopesfc
4498 !!$if( abs(slopeez) < 1e-14 ) print*,"small slopeez =",slopeez
4499 !!$if( abs(sigmasfc) < 1e-14 ) print*,"small sigmasfc =",sigmasfc
4500 !!$if( abs(sigmaez) < 1e-14 ) print*,"small sigmaez =",sigmaez
4503 slopeSfc = sign(max( abs(slopeSfc), 1e-15 ), slopeSfc)
4504 !!slopeEZ = sign(max( abs(slopeEZ), 1e-10 ), slopeEZ) !1e-15 caused an overflow for the jfd~
4505 if(slopeEZ > 2000) then
4507 else if(slopeEZ < -2000) then
4509 else if(slopeEZ < 10 .and. slopeEZ > 0) then
4511 else if(slopeEZ < 0 .and. slopeEZ > -10.0) then
4514 sigmaSfc = sign(max( abs(sigmaSfc), 1e-15 ), sigmaSfc)
4515 sigmaEZ = sign(max( abs(sigmaEZ), 1e-15 ), sigmaEZ)
4516 !!slopeEZ = 1000.0 ! Larry, set constant value of slopeEZ
4517 !! slopeSfc = sign(min (abs(slopeSfc), 5000.0), slopeSfc) ! lkb Added check on size of slopes
4518 !! slopeSfc = sign(min (abs(slopeEZ), 5000.0), slopeEZ)
4520 ! Calculate all the values that are held constant while looping through
4521 ! the perturbations...
4523 centerBin = numBins / 2 + 1 ! Find the center bin
4524 dslope = sign(max(abs(slopeEZ-slopeSfc),1e-15),slopeEZ-slopeSfc)
4525 jacobian = slopeEZ / dslope ! Compute the jacobian
4526 !wig: 22-Dec-2006 added parentheses that had been inadvertantly dropped...
4527 !wig constants = jacobian*thBinSize*rBinSize / (pi2*sigmaSfc*sigmaEZ)
4528 bigcheck = sqrt(huge(c)) ! 10/30/08 lkb 0.1*huge(c)
4530 ! Loop through all the perturbation possibilities and get the jfd...
4533 do j = 1, numBins ! For each bin of the jfd
4534 r_perturb(j) = rBinSize * (j - centerBin)
4536 th_perturb(i) = thBinSize * (i - centerBin)
4538 ! Convert theta and r to c and m space. This uses eq. 4
4539 ! from Berg and Stull (2004)
4540 c = slopeEZ * (th_perturb(i) - slopeSfc * r_perturb(j)) / dslope
4541 m = (th_perturb(i) - slopeEZ * r_perturb(j)) / dslope
4543 !wig, 22-Dec-2006: Actual desired calc commented since was getting
4544 ! an overflow. So, added code to enforce limits.
4545 ! jfd(i,j) = exp(-0.5 * ( (m/sigmaSfc) * (m/sigmaSfc) + &
4546 ! (c/sigmaEZ) * (c/sigmaEZ) )) * constants
4548 if( abs(cterm) > bigcheck ) then
4550 '("KFCuP setting a bogus cterm for JFD. c=",1e15.6," &
4551 & sigmaEZ=",1e15.6)') &
4553 call wrf_debug(0,trim(message))
4554 cterm = sign(bigcheck,cterm)
4559 !!$ if( abs(mterm) > 0.1*bigcheck ) then
4560 !!$ write(message, &
4561 !!$ '("KFCuP has a big mterm for JFD. m=",1e15.6," sigmaSfc=",1e15.6," dslope=",1e15.6," slopeEZ=",1e15.6," slopeSfc=",1e15.6)') &
4562 !!$ m, sigmaSfc,dslope,slopeEZ,slopeSfc
4563 !!$ call wrf_debug(0,trim(message))
4567 if( abs(mterm) > bigcheck ) then
4569 print*,'bigcheck=',bigcheck
4571 '("KFCuP setting a bogus mterm for JFD. m=",1e15.6, &
4572 & " sigmaSfc=",1e15.6)') &
4574 call wrf_debug(0,trim(message))
4577 mterm = sign(bigcheck,mterm)
4581 !wig: took off constants because they will not affect the outcome after normalizing to one
4582 jfd(i,j) = exp( -0.5*(mterm + cterm) ) !* constants
4583 !wig: end of overflow hack
4585 jfdCheckSum = jfdCheckSum + jfd(i,j)
4590 !!$!~Add check to only output the check sum if it is out of the ordinary...
4591 !!$ write(*,*) "JFD sums to ", jfdCheckSum, " Number of bins is ", numBins
4592 !!$ write(30,*) "~JFD sums to ", jfdCheckSum, " Number of bins is ", numBins
4593 !!$ write(30,'("slope sfc/ez & sigma sfc/ez: ",4g18.8)') slopesfc,slopeez,sigmasfc,sigmaez
4594 !!$ if( count(abs(jfd) > 1e-30) > 1 ) write(30,*) "---Non-spiked JFD---",count(abs(jfd) > 1e-30)
4595 !!$ write(30,'(21g11.4)') 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21
4597 !!$ write(30,'(i3, 17e11.4)') j,jfd(:,j)
4600 ! Force jfd sum to be one...
4601 if( jfdCheckSum /= 0. ) jfd(:,:) = jfd(:,:)/jfdCheckSum !~Re-normalize the jfd to sum to one
4603 !!$ write(30,*) "~adjusted JFD..."
4604 !!$ write(30,'(21g11.4)') 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21
4606 !!$ write(30,'(i3, 21e11.4)') j,jfd(:,j)
4610 END SUBROUTINE cup_jfd
4613 !------------------------------------------------------------------------
4614 SUBROUTINE cupSlopeSigma(dx, psfc, p, rho, dz8w, z, ht, &
4615 t, th, tsk, u, v, qv_curr, hfx,xland, qfxin, mavail, & ! add xland, LD 19-Oct-2011
4616 sf_sfclay_physics, br, regime, pblh, kpbl, t2, q2, &
4617 slopeSfc, slopeEZ, sigmaSfc, sigmaEZ, wStar, cupflag, &
4618 shall, kms, kme, kts, kte )
4620 USE module_model_constants, only: cp, ep_1, ep_2, g, r_d, rcp, &
4621 svp1, svp2, svp3, svpt0, xlv
4623 USE module_state_description, ONLY : KFCUPSCHEME &
4630 ! MPI is needed for the test printouts (to get the rank)...
4631 !#ifdef ( DM_PARALLEL ) && !defined( STUBMPI )
4632 #if ( ! defined(DM_PARALLEL) && ! defined(STUBMPI) )
4633 ! rce_testing turn this off
4640 integer, intent(in) :: kpbl, sf_sfclay_physics, &
4643 real, intent(in) :: &
4644 br, dx, hfx,xland, ht, mavail, pblh, psfc, q2, qfxin, regime, t2, tsk !add xland LD 19-Oct-2011
4646 real, dimension(kms:kme), intent(in) :: &
4647 p, rho, dz8w, z, t, th, qv_curr, u, v
4649 real, intent(out) :: &
4650 slopeSfc, slopeEZ, sigmaSfc, sigmaEZ, wStar
4651 real, intent(inout) :: shall
4653 logical, intent(out) :: cupflag
4657 integer :: docldstep, fout, i, ierr, j, k, kpblmid, numZ
4658 real :: br2, dtcu, e1, dthvdz, flux, govrth, psfccmb, qdiff, qfx, &
4659 qsfc, rhox, thv2, thgb, thv, tskv, tvcon, vconv, vsgd, wspd, za
4660 real, dimension(kts:kte) :: zagl
4661 logical :: UnstableOrNeutral
4662 character(len=50) :: filename
4664 ! Artificially force a latent heat flux that is not close to zero. This
4665 ! prevents sigmaSfc from becoming too small and leading to overflows
4666 ! in the JFD calculation.
4668 if( abs(qfxin) < 1./xlv ) then
4669 qfx = sign(1./xlv,qfxin)
4674 ! Determine if each column is stable or (unstable or neutral). If the regime
4675 ! is already calculated by one of the surface schemes, we can use it. If not,
4676 ! deterimine the stability based on the bulk richardson number. We only
4677 ! care about stable vs. (neutral or unstable).
4679 UnstableOrNeutral = .false.
4680 sfclay_case: SELECT CASE (sf_sfclay_physics)
4682 ! Regime categories:
4683 ! 1 = Stable (nighttime)
4684 ! 2 = Damped mechanical turbulence
4685 ! 3 = Forced convection
4686 ! 4 = Free convection
4687 ! Add condition for positive heat flux because negative heat fluxes
4688 ! were causing the wstar calculation to core dump--can't do a 1/3
4689 ! root of a negative value. wig, 5-Feb-2008
4691 .AND. hfx >= 0. ) UnstableOrNeutral = .true.
4694 if( br <= 0. ) UnstableOrNeutral = .true.
4697 ! The selected sfc scheme does not already provide a stability
4698 ! criteria. So, we will mimic the bulk Richardson calculation from
4699 ! module_sf_sfclay.F.
4701 !!if( pblh <= 0. ) call wrf_error_fatal( &
4702 !! "CuP needs a PBL height from a PBL scheme.")
4704 UnstableOrNeutral = .false. ! Added by LKB 9/8/09
4706 else ! Added by LKB 9/8/09
4709 E1 = SVP1*EXP(SVP2*(TSK-SVPT0)/(TSK-SVP3))
4710 PSFCCMB=PSFC/1000. !converts from Pa to cmb
4711 QSFC = EP_2*E1/(PSFCCMB-E1)
4712 THGB = TSK*(100./PSFCCMB)**RCP
4713 TSKV = THGB*(1.+EP_1*QSFC*MAVAIL)
4714 TVCON = 1.+EP_1*QV_CURR(1)
4720 RHOX = PSFC/(r_d*t(1)*TVCON)
4721 flux = max(hfx/rhox/cp + ep_1*tskv*qfx/rhox,0.)
4722 VCONV = (g/TSK*pblh*flux)**.33
4723 VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
4724 WSPD = SQRT(U(1)*U(1)+V(1)*V(1))
4725 WSPD = SQRT(WSPD*WSPD+VCONV*VCONV+vsgd*vsgd)
4726 WSPD = MAX(WSPD,0.1)
4728 !And finally, the bulk Richardson number...
4729 BR2 = GOVRTH*ZA*DTHVDZ/(WSPD*WSPD)
4731 if( br2 <= 0. ) then
4732 UnstableOrNeutral = .true.
4734 UnstableOrNeutral = .false.
4737 END SELECT sfclay_case
4739 ! If we are in a stable regime, then the assumptions for CuP do not
4740 ! make sense, so default back to the standard KF algorithm. Also, do
4741 ! this if the pbl is at the lowest level since then we cannot
4742 !calculate a proper difference with the surface.
4743 !if( kpbl == 1 .or. (.not. UnstableOrNeutral) .or. hfx < 0 .or. qfx < 0) then !~ lkb 8/25/08 changed to require + heat flux
4744 ! if( kpbl == 1 .or. hfx < 1 ) then !~ lkb 8/25/08 changed to require + heat flux
4745 ! if( kpbl == 1 .or. hfx < 50 ) then !~ lkb 8/25/08 changed to require + heat flux
4746 ! if( kpbl <= 2 .or. hfx < 100 ) then !~ lkb 8/25/08 changed to require + heat flux
4747 if(xland .eq.1 )then !~ LD 18-Oct-2011
4748 if( kpbl <= 2 .or. hfx < 100 ) then !~ lkb 8/25/08 changed to require + heat flux
4754 shall = 2 ! Added by LK Berg on 6/17/09 to stop shallow clouds at night
4755 return ! <---Alternate return point
4760 if( kpbl <= 2 .or. hfx < 1 ) then !~ lkb 8/25/08 changed to require + heat flux
4766 shall = 2 ! Added by LK Berg on 6/17/09 to stop shallow clouds at night
4767 return ! <---Alternate return point
4773 ! Convert height from AMSL to AGL...
4778 !!$ ! Find the index closest to the middle of the PBL...
4781 !!$ if( zagl(k) > pblh(i,j) ) then
4782 !!$ kpblmid = max(1, k/2)
4786 !!$ if( kpblmid == 0 ) &
4787 !!$ call wrf_error("CuP ERROR: PBLH not within the domain.")
4789 if( kpbl == 0 ) call wrf_error_fatal("CuP ERROR: kpbl==0")
4791 ! Calculate the Deardorff velocity, wStar. As a rough
4792 ! approximation of the middle of PBL averaged theta and mixing
4793 ! ratio, use the value at the middle of the PBL.
4794 ! The flux amalgamation formula is from Stull, p.147 and
4795 ! wStar is from Stull, p. 118.
4796 kpblmid = max(kts,kpbl/2)
4797 flux = (1. + EP_1*qv_curr(1))*hfx/rho(1)/cp + &
4798 EP_1*th(1)*qfx/rho(1) !badbad/xlv
4799 tvcon = 1.+EP_1*qv_curr(kpblmid)
4800 thv = th(kpblmid)*tvcon
4801 wStar = (g*pblh*flux/thv)**(1./3.)
4802 !!write(*,*) 'Larry ... wStar', wStar, pblh, flux, thv
4803 ! Calculate the slope (dTemp/dMixRatio) for the surface layer
4804 ! and entrainment zone...
4805 thv = th(kpblmid)*tvcon !Virt. pot. temp. at lowest model level
4806 tvcon = 1.+EP_1*qv_curr(1)
4807 thv2 = th(1)*tvcon !Virt. pot. temp. at lowest model level
4808 qdiff = qv_curr(kpblmid)-qv_curr(1)
4809 if( abs(qdiff) < reallysmall ) qdiff = sign(reallysmall,qdiff)
4810 ! slopeSfc = (thv-thv2) / qdiff
4811 ! Changed slopeSfc to use Bowen ratio
4812 slopeSfc = hfx/(xlv * qfx) * xlv / cp ! Recall that hfx is in W m-2 and LH is also in W m-2
4813 tvcon = 1.+EP_1*qv_curr(min(kpbl+2,kte))
4814 thv = th(min(kpbl+2,kte))*tvcon
4815 tvcon = 1.+EP_1*qv_curr(kpblmid)
4816 thv2 = th(kpblmid)*tvcon
4817 qdiff = qv_curr(min(kpbl+2,kte)) - qv_curr(kpblmid)
4818 if( abs(qdiff) < reallysmall ) then
4819 qdiff = sign(reallysmall,qdiff)
4821 slopeEZ = (thv-thv2) / qdiff
4822 ! Calculate the standard deviations along the theta and
4823 ! mixing ratio axes of the PDF following Berg and Stull (2004)
4824 ! eqs. 17a and 17b. For sigmaSfc, we currently are only using
4825 ! rstar and not rstarNew.
4826 ! For sigmaEZ, reuse the flux var that contains (w'thetav')bar
4827 sigmaEz = flux/wStar* &
4829 (zagl(kpblmid)/pblh)**(-1.8) ) ! Changed by lkb, 1/21/09 to use kpblmid
4830 !!(zagl(min(kpbl+2,kte))/pblh)**(-1.8) )
4832 flux = qfx/rho(1) !badbad /xlv ! (w'qv')bar
4833 !!$ sigmaSfc(i,j) = flux*(1-zagl(1)/pblh(i,j))/wStar * &
4834 !!$ ( 2.3 + 1.1e-2*(zagl(1)/pblh(i,j))**(-1.6) )
4835 sigmaSfc = flux/wStar * &
4836 ( 2.3 + 1.1e-2*(zagl(kpblmid)/pblh)**(-1.6) )
4840 ! Output the inputs to CuP for debugging with offline code...
4842 call wrf_message("Outputting cupin file.")
4844 #ifdef ( DM_PARALLEL ) && !defined( STUBMPI )
4845 CALL MPI_Comm_rank ( MPI_COMM_WORLD, k, ierr ) !this isn't tested with MPI yet
4847 write(filename, '("cupin.",i3.3,".txt")') k
4849 do !Make sure we use an available unit.
4850 inquire(UNIT=fout,OPENED=ierr)
4851 if( ierr==.true. ) exit
4853 if( fout > 100 ) exit
4855 open(UNIT=fout, FILE=trim(filename), FORM="formatted", &
4856 STATUS="unknown", IOSTAT=k)
4857 if( k /= 0 ) call wrf_error_fatal("Could not open cupin file.")
4859 write(UNIT=fout,FMT='(a)') "Inputs to cup_driver..."
4860 !!$ write(UNIT=fout,FMT='("ktau,i,j=",i,2i5)') ktau, i, j
4861 !!$ write(UNIT=fout,FMT='("stepcu, dt=",i,g17.9)') stepcu, dt
4862 !!$ write(UNIT=fout,FMT='("ids,ide, jds, jde, kds, kde=",6i5)') &
4863 !!$ ids,ide, jds, jde, kds, kde
4864 !!$ write(UNIT=fout,FMT='("ims,ime, jms, jme, kms, kme=",6i5)') &
4865 !!$ ims,ime, jms, jme, kms, kme
4866 !!$ write(UNIT=fout,FMT='("its,ite, jts, jte, kts, kte=",6i5)') &
4867 !!$ its,ite, jts, jte, kts, kte
4869 write(UNIT=fout,FMT='("sf_sfclay_physics =",i)') sf_sfclay_physics
4870 write(UNIT=fout,FMT='("dx =",g17.9)') dx
4871 write(UNIT=fout,FMT='("psfc =",g17.9)') psfc
4872 write(UNIT=fout,FMT='("kpbl =",i)') kpbl
4873 write(UNIT=fout,FMT='("pblh =",g17.9)') pblh
4874 write(UNIT=fout,FMT='("ht =",g17.9)') ht
4875 write(UNIT=fout,FMT='("tsk =",g17.9)') tsk
4876 write(UNIT=fout,FMT='("t2 =",g17.9)') t2
4877 write(UNIT=fout,FMT='("q2 =",g17.9)') q2
4878 write(UNIT=fout,FMT='("hfx =",g17.9)') hfx
4879 write(UNIT=fout,FMT='("qfx =",g17.9)') qfx
4880 write(UNIT=fout,FMT='("mavail =",g17.9)') mavail
4881 write(UNIT=fout,FMT='("br =",g17.9)') br
4882 write(UNIT=fout,FMT='("regime =",g17.9)') regime
4884 write(UNIT=fout,FMT='("p,rho, t, th, qv:")')
4886 write(UNIT=fout,FMT='(" ",5g17.9)') &
4887 p(k), rho(k), t(k), th(k), qv_curr(k)
4890 write(UNIT=fout,FMT='("z, dz8w, u, v:")')
4892 write(UNIT=fout,FMT='(" ",4g17.9)') &
4893 z(k), dz8w(k), u(k), v(k)
4896 write(UNIT=fout,FMT='(a)') "Calculated inside cup_driver..."
4897 write(UNIT=fout,FMT='("slopeSfc =",g17.9)') SlopeSfc
4898 write(UNIT=fout,FMT='("slopeEZ =",g17.9)') SlopeEZ
4899 write(UNIT=fout,FMT='("sigmaSfc =",g17.9)') sigmaSfc
4900 write(UNIT=fout,FMT='("sigmaEZ =",g17.9)') sigmaEZ
4901 write(UNIT=fout,FMT='("wStar =",g17.9)') wStar
4902 write(UNIT=fout,FMT='("dtcu =",g17.9)') dtcu
4904 write(UNIT=fout,FMT='("zagl:")')
4906 write(UNIT=fout,FMT='(" ",1g17.9)') &
4912 END SUBROUTINE cupSlopeSigma
4915 !------------------------------------------------------------------------
4916 ! Find Cp for moist air
4921 real, intent(in) :: r ! Mixing ratio
4923 findCp = 1004.67 * (1.0 + 0.84 * r)
4927 !------------------------------------------------------------------------
4928 ! Finds the index when an ordered list becomes bigger than a given value.
4929 ! The list is assumed to be ordered from small to big values.
4930 FUNCTION findIndex(value,list)
4932 integer :: findindex
4933 real, intent(in) :: value
4934 real, intent(in), dimension(:) :: list
4939 do i=1,ubound(list,1)
4940 if( value <= list(i) ) then
4945 END FUNCTION findIndex
4947 !------------------------------------------------------------------------
4948 ! Find the saturation mixing ratio w.r.t. water. This subroutine uses
4950 ! T in K and p in hPa
4951 FUNCTION findRs(t,p)
4953 real, intent(in) :: t, p
4956 es = 610.78 * exp( 17.67 * (t - 273.16) / (t - 29.66))
4957 findRs = eps * es / (p - es)
4961 !------------------------------------------------------------------------
4962 ! Find the saturation mixing ratio w.r.t. ice.
4963 ! T in K and p in hPa
4964 FUNCTION findRsi(t,p)
4966 real, intent(in) :: t, p
4970 ! esi = 10.**(-9.09685*(273.15/t - 1.) - 3.56654*log10(273.15/t) &
4971 ! + 0.87682*(1. - t/273.15) + 0.78614)
4973 ! GoffGratch formula:
4974 esi = 10**(-9.09718*(273.15/t - 1.) - 3.56654*log10(273.15/t) &
4975 + 0.876793*(1. - t/273.15) + log10(6.1071))
4977 findRsi = eps * esi / (p - esi)
4978 END FUNCTION findRsi
4981 !------------------------------------------------------------------------
4982 subroutine activate_cldbase_kfcup( idiagee, grid_id, ktau, &
4983 ii, jj, kk, kts, kte, lc, kcldlayer, &
4984 num_chem, maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
4985 ntype_aer, nsize_aer, ncomp_aer, &
4986 ai_phase, msectional, massptr_aer, numptr_aer, &
4987 dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer, &
4988 tk_act, rho_act, dp, w_act, &
4989 chem1d, qndrop_act )
4991 use module_mixactivate, only: activate
4993 integer, intent(in) :: &
4994 idiagee, grid_id, ktau, &
4995 ii, jj, kk, kts, kte, lc, kcldlayer, &
4996 num_chem, maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
4997 msectional, ntype_aer, ai_phase
4998 integer, intent(in) :: ncomp_aer(maxd_atype), nsize_aer(maxd_atype)
4999 integer, intent(in) :: massptr_aer(maxd_acomp,maxd_asize,maxd_atype,maxd_aphase)
5000 integer, intent(in) :: numptr_aer(maxd_asize,maxd_atype,maxd_aphase)
5002 real, intent(in ) :: chem1d(kts:kte,1:num_chem)
5003 real, intent(in ) :: dens_aer(maxd_acomp,maxd_atype)
5004 real, intent(in ) :: dlo_sect(maxd_asize,maxd_atype), dhi_sect(maxd_asize,maxd_atype)
5005 real, intent(in ) :: dp(kts:kte)
5006 real, intent(in ) :: hygro_aer(maxd_acomp,maxd_atype)
5007 real, intent(inout) :: qndrop_act
5008 real, intent(in ) :: rho_act
5009 real, intent(in ) :: sigmag_aer(maxd_asize,maxd_atype)
5010 real, intent(in ) :: tk_act
5011 real, intent(in ) :: w_act
5013 integer :: icomp, iphase, isize, itype, k, l
5015 real :: flux_fullact
5016 real :: tmpa, tmpdpsum, tmpvol, tmpwght
5017 real, dimension( 1:maxd_asize, 1:maxd_atype ) :: &
5018 fn, fs, fm, fluxn, fluxs, fluxm, &
5019 hygroavg, numbravg, volumavg
5022 ! for each isize and itype, calculate average number, volume, and hygro
5023 ! over the updraft source layers
5025 ! if (idiagee > 0) write(98,'(//a,5i5)') 'kfcup activate_cldbase_kfcup - i, j, ksrc1/2', i, j, lc, kcldlayer
5029 tmpdpsum = sum( dp(lc:kcldlayer) )
5031 do k = lc, kcldlayer
5032 tmpwght = dp(k)/tmpdpsum
5033 do itype = 1, ntype_aer
5034 do isize = 1, nsize_aer(itype)
5035 l = numptr_aer(isize,itype,iphase)
5036 numbravg(isize,itype) = numbravg(isize,itype) + tmpwght*max( 0.0, chem1d(k,l) )
5037 do icomp = 1, ncomp_aer(itype)
5038 l = massptr_aer(icomp,isize,itype,iphase)
5039 tmpvol = max( 0.0, chem1d(k,l) ) / dens_aer(icomp,itype)
5040 volumavg(isize,itype) = volumavg(isize,itype) + tmpwght*tmpvol
5041 hygroavg(isize,itype) = hygroavg(isize,itype) + tmpwght*tmpvol*hygro_aer(icomp,itype)
5047 do itype = 1, ntype_aer
5048 do isize = 1, nsize_aer(itype)
5049 hygroavg(isize,itype) = hygroavg(isize,itype) / max( 1.0e-35, volumavg(isize,itype) )
5050 ! convert numbravg from (#/kg) to (#/m3)
5051 numbravg(isize,itype) = numbravg(isize,itype)*rho_act
5052 ! convert volumavg to (m3/m3) -- need 1e-12 factor because (rho_act*chem1d)/dens_aer = [(ugaero/m3air)/(gaero/cm3aero)]
5053 volumavg(isize,itype) = volumavg(isize,itype)*rho_act*1.0e-12
5055 ! if (vaero_dsect_adjust_opt == 1) then
5056 ! recalc volumavg using particle diameter = dcen_sect
5057 ! tmpvol = sqrt( dlo_sect(isize,itype) * dhi_sect(isize,itype) ) * 1.0e-2 ! particle diameter in (m)
5058 ! tmpvol = (tmpvol**3) * 3.1415926536/6.0 ! particle volume in (m3)
5059 ! volumavg(isize,itype) = numbravg(isize,itype) * tmpvol
5064 ! adjust number and volume for scm sensitivity testing
5065 ! numbravg(:,:) = numbravg(:,:) * max( naero_adjust_factor, 1.0e-2 )
5066 ! volumavg(:,:) = volumavg(:,:) * max( naero_adjust_factor, 1.0e-2 )
5068 call activate( w_act, 0.0, 0.0, 0.0, 1.0, tk_act, rho_act, &
5069 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, &
5070 numbravg, volumavg, dlo_sect, dhi_sect, sigmag_aer, hygroavg, &
5071 fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
5072 grid_id, ktau, ii, jj, kk )
5074 ! subroutine activate( wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
5075 ! msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, &
5076 ! na, volc, dlo_sect, dhi_sect, sigman, hygro, &
5077 ! fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
5078 ! grid_id, ktau, ii, jj, kk )
5082 ! integer,intent(in) :: maxd_atype ! dimension of types
5083 ! integer,intent(in) :: maxd_asize ! dimension of sizes
5084 ! integer,intent(in) :: ntype_aer ! number of types
5085 ! integer,intent(in) :: nsize_aer(maxd_atype) ! number of sizes for type
5086 ! integer,intent(in) :: msectional ! 1 for sectional, 0 for modal
5087 ! integer,intent(in) :: grid_id ! WRF grid%id
5088 ! integer,intent(in) :: ktau ! WRF time step count
5089 ! integer,intent(in) :: ii, jj, kk ! i,j,k of current grid cell
5090 ! real,intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
5091 ! real,intent(in) :: sigw ! subgrid standard deviation of vertical vel (m/s)
5092 ! real,intent(in) :: wdiab ! diabatic vertical velocity (0 if adiabatic)
5093 ! real,intent(in) :: wminf ! minimum updraft velocity for integration (m/s)
5094 ! real,intent(in) :: wmaxf ! maximum updraft velocity for integration (m/s)
5095 ! real,intent(in) :: tair ! air temperature (K)
5096 ! real,intent(in) :: rhoair ! air density (kg/m3)
5097 ! real,intent(in) :: na(maxd_asize,maxd_atype) ! aerosol number concentration (/m3)
5098 ! real,intent(in) :: sigman(maxd_asize,maxd_atype) ! geometric standard deviation of aerosol size distribution
5099 ! real,intent(in) :: hygro(maxd_asize,maxd_atype) ! bulk hygroscopicity of aerosol mode
5100 ! real,intent(in) :: volc(maxd_asize,maxd_atype) ! total aerosol volume concentration (m3/m3)
5101 ! real,intent(in) :: dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm)
5102 ! dhi_sect( maxd_asize, maxd_atype ) ! maximum size of section (cm)
5105 ! real,intent(inout) :: fn(maxd_asize,maxd_atype) ! number fraction of aerosols activated
5106 ! real,intent(inout) :: fs(maxd_asize,maxd_atype) ! surface fraction of aerosols activated
5107 ! real,intent(inout) :: fm(maxd_asize,maxd_atype) ! mass fraction of aerosols activated
5108 ! real,intent(inout) :: fluxn(maxd_asize,maxd_atype) ! flux of activated aerosol number fraction into cloud (m/s)
5109 ! real,intent(inout) :: fluxs(maxd_asize,maxd_atype) ! flux of activated aerosol surface fraction (m/s)
5110 ! real,intent(inout) :: fluxm(maxd_asize,maxd_atype) ! flux of activated aerosol mass fraction into cloud (m/s)
5111 ! real,intent(inout) :: flux_fullact ! flux when activation fraction = 100% (m/s)
5114 do itype = 1, ntype_aer
5115 do isize = 1, nsize_aer(itype)
5116 qndrop_act = qndrop_act + numbravg(isize,itype)*fn(isize,itype)
5117 tmpa = max( numbravg(isize,itype), max(volumavg(isize,itype),1.0)*1.0e-30 )
5118 tmpa = (6.0*volumavg(isize,itype)/(3.1415926536*tmpa))**0.33333333
5119 if (idiagee > 0) write(98,'(a,2i3,1p,9e10.2)') 'bin, numbr, volum, hygro, sg, dlo, dav, dhi, fn, fm', itype, isize, &
5120 numbravg(isize,itype), volumavg(isize,itype), hygroavg(isize,itype), &
5121 sigmag_aer(isize,itype), 0.01*dlo_sect(isize,itype), tmpa, 0.01*dhi_sect(isize,itype), &
5122 fn(isize,itype), fm(isize,itype)
5125 qndrop_act = qndrop_act/rho_act
5126 if (idiagee > 0) write(98,'(a,21x,i6,1p,2e10.2)') 'msectional, w_act, qndrop', msectional, w_act, qndrop_act
5129 end subroutine activate_cldbase_kfcup
5132 !------------------------------------------------------------------------
5133 subroutine adjust_mfentdet_kfcup( idiagee, grid_id, ktau, &
5134 ii, jj, kts, kte, kcutop, ishall, &
5135 umfout, uerout, udrout, dmfout, derout, ddrout )
5137 integer, intent(in) :: &
5138 idiagee, grid_id, ktau, &
5139 ii, jj, kts, kte, kcutop, ishall
5141 real, dimension( kts:kte ), intent(inout) :: &
5142 umfout, uerout, udrout, dmfout, derout, ddrout
5145 real, parameter :: rtol = 1.0e-6
5146 real :: tmpa, tmpb, tmpc, tmpf, tmpg, tmph, tmpold
5149 ! check that delta(dmfout) = derout - ddrout
5150 ! if not, then adjust either derout or ddrout
5151 ! the diagnostic output shows these adjustments to be very small,
5152 ! so this may be unnecessary
5154 if (ishall == 0) then
5156 dmfout(kcutop:kte) = 0.0
5157 if (kcutop < kte) then
5158 derout(kcutop+1:kte) = 0.0
5159 ddrout(kcutop+1:kte) = 0.0
5166 tmpa = dmfout(k) - dmfout(k-1)
5170 tmpb = derout(k) - ddrout(k)
5172 if (tmpc > 0.0) then
5173 if (derout(k) < ddrout(k)*0.05) then
5174 ! der << ddr, so decrease ddr first, then increase der if needed
5176 ddrout(k) = max( 0.0, ddrout(k) - tmpc )
5177 tmpg = tmpg + abs(ddrout(k)-tmpold)
5178 tmpb = derout(k) - ddrout(k)
5180 derout(k) = derout(k) + tmpc
5181 tmpg = tmpg + abs(tmpc)
5184 derout(k) = derout(k) + tmpc
5185 tmpg = tmpg + abs(tmpc)
5188 if (ddrout(k) <= derout(k)*0.05) then
5189 ! ddr << der, so decrease der first, then increase ddr if needed
5191 derout(k) = max( 0.0, derout(k) + tmpc )
5192 tmpg = tmpg + abs(derout(k)-tmpold)
5193 tmpb = derout(k) - ddrout(k)
5195 ddrout(k) = ddrout(k) - tmpc
5196 tmpg = tmpg + abs(tmpc)
5199 ddrout(k) = ddrout(k) - tmpc
5200 tmpg = tmpg + abs(tmpc)
5205 if ( idiagee > 0 ) then
5206 tmpf = sum(derout(kts:kcutop)) + sum(ddrout(kts:kcutop))
5207 tmph = tmpg/max(tmpg,tmpf,1.0e-20)
5208 if (abs(tmph) > rtol) &
5209 write(*,'(a,i9,2i5,1p,4e10.2)') 'kfcupmfadjup', ktau, ii, jj, &
5210 minval(dmfout(kts:kcutop)), tmpf, tmpg, tmph
5215 ! check that delta(umfout) = uerout - udrout
5216 ! if not, then adjust either uerout or udrout
5217 ! the diagnostic output shows these adjustments to mostly be very small,
5218 ! but there is an occasional problem at klcl,
5219 ! this suggests a problem in the code that calculates umf and uer,
5220 ! but i have not been able to locate it, so this bandaid is needed
5222 if ((ishall == 0) .or. (ishall == 1)) then
5224 umfout(kcutop:kte) = 0.0
5225 if (kcutop < kte) then
5226 uerout(kcutop+1:kte) = 0.0
5227 udrout(kcutop+1:kte) = 0.0
5233 tmpa = umfout(k) - umfout(k-1)
5237 tmpb = uerout(k) - udrout(k)
5239 if (tmpc > 0.0) then
5240 if (uerout(k) < udrout(k)*0.05) then
5241 ! uer << udr, so decrease udr first, then increase uer if needed
5243 udrout(k) = max( 0.0, udrout(k) - tmpc )
5244 tmpg = tmpg + abs(udrout(k)-tmpold)
5245 tmpb = uerout(k) - udrout(k)
5247 uerout(k) = uerout(k) + tmpc
5248 tmpg = tmpg + abs(tmpc)
5251 uerout(k) = uerout(k) + tmpc
5252 tmpg = tmpg + abs(tmpc)
5255 if (udrout(k) <= uerout(k)*0.05) then
5256 ! udr << uer, so decrease uer first, then increase udr if needed
5258 uerout(k) = max( 0.0, uerout(k) + tmpc )
5259 tmpg = tmpg + abs(uerout(k)-tmpold)
5260 tmpb = uerout(k) - udrout(k)
5262 udrout(k) = udrout(k) - tmpc
5263 tmpg = tmpg + abs(tmpc)
5266 udrout(k) = udrout(k) - tmpc
5267 tmpg = tmpg + abs(tmpc)
5272 if ( idiagee > 0 ) then
5273 tmpf = sum(uerout(kts:kcutop)) + sum(udrout(kts:kcutop))
5274 tmph = tmpg/max(tmpg,tmpf,1.0e-20)
5275 if (abs(tmph) > rtol) &
5276 write(*,'(a,i9,2i5,1p,4e10.2)') 'kfcupmfadjup', ktau, ii, jj, &
5277 maxval(umfout(kts:kcutop)), tmpf, tmpg, tmph
5284 end subroutine adjust_mfentdet_kfcup
5287 ! rce 11-may-2012 mods start -------------------------------------------
5288 subroutine cu_kfcup_diagee01( &
5289 ims, ime, jms, jme, kms, kme, kts, kte, &
5291 idiagee, idiagff, ishall, ktau, &
5292 kcubotmin, kcubotmax, kcutopmin, kcutopmax, &
5293 activefrac, cldfra_cup1d, &
5294 cubot, cutop, cumshallfreq1d, &
5295 ddr_deep, der_deep, dmf_deep, dt, dz1d, &
5296 fcvt_qc_to_pr_deep, fcvt_qc_to_qi_deep, fcvt_qi_to_pr_deep, &
5297 fcvt_qc_to_pr_shall, fcvt_qc_to_qi_shall, fcvt_qi_to_pr_shall, &
5298 nca_deep, nca_shall, p1d, pblh, &
5299 qc_ic_deep, qc_ic_shall, qi_ic_deep, qi_ic_shall, qndrop_ic_cup, rho1d, &
5300 tactive, taucloud, tstar, &
5301 udr_deep, udr_shall, uer_deep, uer_shall, umf_deep, umf_shall, &
5302 updfra_deep, updfra_shall, updfra_cup, &
5303 wact_cup, wcloudbase, wcb_v2, wcb_v2_shall, &
5304 wulcl_cup, wstar, z1d, z_at_w1d )
5307 integer, intent(in) :: &
5308 ims, ime, jms, jme, kms, kme, kts, kte, &
5310 idiagee, idiagff, ishall, ktau, &
5311 kcubotmin, kcubotmax, kcutopmin, kcutopmax
5313 real, intent(in) :: &
5323 real, dimension( kts:kte ), intent(in) :: &
5330 fcvt_qc_to_pr_deep, &
5331 fcvt_qc_to_pr_shall, &
5332 fcvt_qc_to_qi_deep, &
5333 fcvt_qc_to_qi_shall, &
5334 fcvt_qi_to_pr_deep, &
5335 fcvt_qi_to_pr_shall, &
5351 real, dimension( ims:ime, jms:jme ), intent(in) :: &
5363 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: &
5372 real :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg, tmph, tmpi, tmpj
5373 real :: tmpr, tmps, tmpx, tmpy, tmpz
5375 real :: tmp_nca, tmp_updfra
5376 real :: tmpveca(1:999)
5379 if (idiagee > 0) then
5383 kcubot = nint(cubot(i,j))
5384 kcutop = nint(cutop(i,j))
5385 k = (kcubot+kcutop)/2
5386 updfra = 0.0 ; if (ishall == 0) updfra = updfra_deep ; if (ishall == 1) updfra = updfra_shall
5387 !!! write(*,'(a,1p,5e11.3,3x,3e11.3)') 'activefrac, cldfra(b/m/t), updfra', &
5388 !!! activefrac(i,j), &
5389 !!! cldfra_cup1d(kcubot), cldfra_cup1d(k), cldfra_cup1d(kcutop), updfra
5390 !!! write(*,'(a,1p,4e11.3,3x,3e11.3)') 'wcb, wcb_v2, wulcl, wact ', &
5391 !!! wcloudbase(i,j), wcb_v2_shall, wulcl_cup(i,j), wact_cup(i,j)
5392 !!! write(*,'(a,1p,5e11.3,3x,3e11.3)') 'qndrop(b/m/t/t-b) ', &
5393 !!! qndrop_ic_cup(kcubot), qndrop_ic_cup(k), qndrop_ic_cup(kcutop), &
5394 !!! qndrop_ic_cup(kcutop)-qndrop_ic_cup(kcubot)
5395 !!! write(*,'(a,4i5,f9.5,10(2x,5i5))') 'updfraprofile*1e4', &
5396 !!! ktau, ishall, kcubot, kcutop, updfra, &
5397 !!! ( nint(updfra_cup(i,k,j)*1.0e4), k=kts,min(kte-1,kcutop+3) )
5399 if ((ishall==1 .or. ishall==0) .and. idiagee>0) then
5400 if (ishall == 1) then
5401 tmp_updfra = updfra_shall
5404 tmp_updfra = updfra_deep
5408 tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0 ; tmpd = 0.0 ; tmpe = 0.0 ; tmpf = 0.0
5410 if (ishall == 1) then
5411 tmpa = tmpa + max( 0.0, uer_shall(k) )
5412 tmpx = cumshallfreq1d(k)
5414 tmpa = tmpa + max( 0.0, uer_deep(k) )
5417 tmpcf = cldfra_cup1d(k)*tmpx
5418 tmpc = tmpc + max( 0.0, tmpcf ) * dz1d(k)*rho1d(k)
5419 ! tmpd = tmpd + max( 0.0, cldfra_cup(i,k,j) ) * dz1d(k)*rho1d(k)
5420 tmpd = tmpd + max( 0.0, cldfra_cup1d(k) ) * dz1d(k)*rho1d(k)
5422 tmpe = tmpe + max( 0.0, tmp_updfra*tmpx ) * dz1d(k)*rho1d(k)
5423 if (kcubot <= k .and. k <= kcutop) &
5424 tmpf = tmpf + max( 0.0, tmp_updfra ) * dz1d(k)*rho1d(k)
5427 tmpb = cldfra_cup1d(kcubot)*wcb_v2*rho1d(kcubot)*tmp_nca
5429 ! if (tmpd > 1.0e-10) tmpg = tmpa/tmpd
5431 ! if (idiagee>0) write(*,'(a,1p,6e11.3,0p,f11.3,i8)') 'entrain mass, cloud-vol mass b-e ', &
5432 ! tmpa, tmpf, tmpd, tmpe, tmpb, tmpc, tmpg, ktau
5433 if (idiagee>0) write(*,'(a,1p,2e11.3,0p,2f9.3,2(3x,1p,2e11.3,0p,f9.3),i8,2(2x,3i3))') 'cloudmassaa ', &
5435 tmpa/max(tmpc,1.0e-10), tmpb/max(tmpc,1.0e-10), &
5436 tmpc, tmpd, max(tmpc,1.0e-10)/max(tmpd,1.0e-10), &
5437 tmpe, tmpf, max(tmpe,1.0e-10)/max(tmpf,1.0e-10), &
5438 ktau, kcubot, kcubotmin, kcubotmax, kcutop, kcutopmin, kcutopmax
5440 tmpi = 0.0 ; tmpj = 0.0
5441 do k = kcubot, kcutop
5442 if (ishall == 1) then
5443 tmpi = tmpi + cldfra_cup1d(k)*dz1d(k)*rho1d(k)*qc_ic_shall(k)
5445 tmpi = tmpi + cldfra_cup1d(k)*dz1d(k)*rho1d(k)*qc_ic_deep(k)
5447 tmpj = tmpj + cldfra_cup1d(k)*dz1d(k)*rho1d(k)
5450 tmpveca(1) = tmpa/max(tmpd,1.0e-10)
5451 tmpveca(2) = tmpb/max(tmpd,1.0e-10)
5452 tmpveca(3) = cldfra_cup1d(kcubot)
5453 tmpveca(4) = sum( dz1d(kcubot:kcutop) )
5456 tmpa = tmpa/tmp_nca ! total inflow
5457 tmpg = tmpd * (tmp_updfra/cldfra_cup1d(kcubot)) ! updraft mass
5458 tmpveca(101) = cldfra_cup1d(kcubot)
5459 tmpveca(102) = tmp_updfra
5460 tmpveca(103) = sum( dz1d(kcubot:kcutop) )
5461 tmpveca(104) = wcb_v2 ! w at cloud base
5462 tmpveca(105) = tmpd ! cloud mass
5463 tmpveca(106) = tmpa ! total inflow
5464 if (ishall == 1) then
5465 tmpveca(107) = umf_shall(max(1,kcubot-1)) ! cloud base inflow
5467 tmpveca(107) = umf_deep(max(1,kcubot-1)) ! cloud base inflow
5469 tmpveca(108) = tmpg/tmpa ! time to "fill" updraft
5470 tmpveca(109) = tmpd/tmpa ! time to "fill" cloud
5471 tmpveca(110) = tactive(i,j) ! active cloud time-scale
5472 tmpveca(111) = taucloud(i,j) ! cloud dissipation time-scale
5473 tmpveca(112) = tstar(i,j) ! boundary layer time-scale
5474 tmpveca(113) = wstar ! boundary layer convective velocity scale
5475 tmpveca(114) = pblh(i,j) ! pbl height (m)
5476 tmpveca(115) = z_at_w1d(kcubot ) - z_at_w1d(kts) ! bottom of cloudbase layer (m agl)
5477 tmpveca(116) = z_at_w1d(kcutop+1) - z_at_w1d(kts) ! top of cloudtop layer (m agl)
5478 tmpveca(117) = (tmpi/max(tmpj,1.0e-30))*1.0e3 ! convert kg/kg to g/kg
5480 tmpveca(106:107) = tmpveca(106:107)*60.0 ! convert kg/m2/s to kg/m2/min
5481 tmpveca(108:112) = tmpveca(108:112)/60.0 ! convert s to min
5482 end if ! ((ishall==1 .or. ishall==0) .and. idiagee>0) then
5484 if (idiagee>0 .and. ishall==1) then
5485 write(*,'(a)') 'k, p, z, dz, umf, del(umf), uer-udr, uer, -udr, qc, qi, f_qc2qi, f_qc2pr, f_qi2pr'
5486 do k = min( kcutop+2, kte-1 ), kts, -1
5487 if (k .eq. kts) then
5490 tmpa = umf_shall(k) - umf_shall(k-1)
5492 tmpb = uer_shall(k) - udr_shall(k)
5493 write(*,'(i2,1p,3e11.3,3x,5e11.3,3x,5e11.3)') &
5494 k, p1d(k), z1d(k), dz1d(k), umf_shall(k), tmpa, tmpb, uer_shall(k), -udr_shall(k), &
5495 qc_ic_shall(k), qi_ic_shall(k), fcvt_qc_to_qi_shall(k), fcvt_qc_to_pr_shall(k), fcvt_qi_to_pr_shall(k)
5497 end if ! (idiagee>0 .and. ishall==1) then
5499 if (idiagee>0 .and. ishall==0) then
5500 write(*,'(a)') 'k, p, z, dz, umf, del(umf), uer-udr, uer, -udr, qc, qi, f_qc2qi, f_qc2pr, f_qi2pr'
5501 do k = min( kcutop+2, kte-1 ), kts, -1
5502 if (k .eq. kts) then
5505 tmpa = umf_deep(k) - umf_deep(k-1)
5507 tmpb = uer_deep(k) - udr_deep(k)
5508 write(*,'(i2,1p,3e11.3,3x,5e11.3,3x,5e11.3)') &
5509 k, p1d(k), z1d(k), dz1d(k), umf_deep(k), tmpa, tmpb, uer_deep(k), -udr_deep(k), &
5510 qc_ic_deep(k), qi_ic_deep(k), fcvt_qc_to_qi_deep(k), fcvt_qc_to_pr_deep(k), fcvt_qi_to_pr_deep(k)
5512 write(*,'(a)') 'k, p, z, dz, dmf, del(dmf), der-ddr, der, -ddr, qc'
5513 do k = min( kcutop+2, kte-1 ), kts, -1
5514 if (k .eq. kts) then
5517 tmpa = dmf_deep(k) - dmf_deep(k-1)
5519 tmpb = der_deep(k) - ddr_deep(k)
5520 write(*,'(i2,1p,3e11.3,3x,5e11.3,3x,5e11.3)') &
5521 k, p1d(k), z1d(k), dz1d(k), dmf_deep(k), tmpa, tmpb, der_deep(k), -ddr_deep(k), qc_ic_deep(k)
5523 end if ! (idiagee>0 .and. ishall==0) then
5525 write(*,'(i6,1p,6e11.3,a)') &
5526 ktau, (ktau*dt/3600.0), tmpveca(1:5), &
5527 ' cloudmassbb ktau, t(h), ratio1, ratio2, cldfra, cldhgt, wcb'
5529 write(*,'(i6,i2, f7.2, 2x,2f8.5,f8.2,2f7.3, 2x,f9.4,2f9.5, 2x,5f8.2, 3f9.1,f9.5, 3a)') &
5530 ktau, ishall, (ktau*dt/3600.0), tmpveca(101:104), tmpveca(113), &
5531 min(9999.99,tmpveca(105)), min(99.99,tmpveca(106:107)), &
5532 min(9999.99,tmpveca(108:112)), min(99999.9,tmpveca(114:116)), min(99.9999,tmpveca(117)), &
5533 ' cloudmasscc ktau,ish,t(h), cldfra,updfra,cldhgt,wcb,wstar', &
5534 ', cldmass,uertot,uerbase, tauinupd,tauincld,tactive,taucloud,tstar', &
5535 ', pblh,zbot,ztop, qc_ic_av'
5537 end if ! (idiagee > 0) then
5540 end subroutine cu_kfcup_diagee01
5541 ! rce 11-may-2012 mods end ---------------------------------------------
5544 END MODULE module_cu_kfcup