source: lmdz_wrf/trunk/WRFV3/phys/module_cu_nsas.F @ 1531

Last change on this file since 1531 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 108.4 KB
Line 
1!
2!
3!
4!
5MODULE module_cu_nsas
6CONTAINS
7!
8!-------------------------------------------------------------------------------
9   subroutine cu_nsas(dt,p3di,p3d,pi3d,qc3d,qi3d,rho3d,itimestep,stepcu,       &
10                     hbot,htop,cu_act_flag,cudt,curr_secs,adapt_step_flag,     &
11                     rthcuten,rqvcuten,rqccuten,rqicuten,                      &
12                     rucuten,rvcuten,                                          &
13                     qv3d,t3d,raincv,pratec,xland,dz8w,w,u3d,v3d,              &
14                     hpbl,hfx,qfx,                                             &
15                     mp_physics,                                               &
16                     p_qc,p_qi,p_first_scalar,                                 &
17                     cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,                      &
18                     cice,xls,psat,f_qi,f_qc,                                  &
19                     ids,ide, jds,jde, kds,kde,                                &
20                     ims,ime, jms,jme, kms,kme,                                &
21                     its,ite, jts,jte, kts,kte)
22!-------------------------------------------------------------------------------
23  implicit none
24!-------------------------------------------------------------------------------
25!
26!-- dt          time step (s)
27!-- p3di        3d pressure (pa) at interface level
28!-- p3d         3d pressure (pa)
29!-- pi3d        3d exner function (dimensionless)
30!-- z           height above sea level (m)
31!-- qc3d        cloud water mixing ratio (kg/kg)
32!-- qi3d        cloud ice mixing ratio (kg/kg)
33!-- qv3d        3d water vapor mixing ratio (kg/kg)
34!-- t3d         temperature (k)
35!-- raincv      cumulus scheme precipitation (mm)
36!-- w           vertical velocity (m/s)
37!-- dz8w        dz between full levels (m)
38!-- u3d         3d u-velocity interpolated to theta points (m/s)
39!-- v3d         3d v-velocity interpolated to theta points (m/s)
40!-- ids         start index for i in domain
41!-- ide         end index for i in domain
42!-- jds         start index for j in domain
43!-- jde         end index for j in domain
44!-- kds         start index for k in domain
45!-- kde         end index for k in domain
46!-- ims         start index for i in memory
47!-- ime         end index for i in memory
48!-- jms         start index for j in memory
49!-- jme         end index for j in memory
50!-- kms         start index for k in memory
51!-- kme         end index for k in memory
52!-- its         start index for i in tile
53!-- ite         end index for i in tile
54!-- jts         start index for j in tile
55!-- jte         end index for j in tile
56!-- kts         start index for k in tile
57!-- kte         end index for k in tile
58!-------------------------------------------------------------------------------
59  integer,  intent(in   )   ::       ids,ide, jds,jde, kds,kde,                &
60                                     ims,ime, jms,jme, kms,kme,                &
61                                     its,ite, jts,jte, kts,kte,                &
62                                     itimestep, stepcu,                        &
63                                     p_qc,p_qi,p_first_scalar
64!
65  real,     intent(in   )   ::      cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,       &
66                                    cice,xls,psat
67  real,     intent(in   )   ::      dt
68!
69  real,     dimension( ims:ime, kms:kme, jms:jme ),optional                  , &
70            intent(inout)   ::                                       rthcuten, &
71                                                                      rucuten, &
72                                                                      rvcuten, &
73                                                                     rqccuten, &
74                                                                     rqicuten, &
75                                                                     rqvcuten
76  logical, optional ::                                              F_QC,F_QI
77  real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
78            intent(in   )   ::                                           qv3d, &
79                                                                         qc3d, &
80                                                                         qi3d, &
81                                                                        rho3d, &
82                                                                          p3d, &
83                                                                         pi3d, &
84                                                                          t3d
85  real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
86            intent(in   )   ::                                           p3di
87  real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
88            intent(in   )   ::                                           dz8w, & 
89                                                                            w
90  real,     dimension( ims:ime, jms:jme )                                    , &
91            intent(inout) ::                                           raincv, &
92                                                                       pratec
93  real,     dimension( ims:ime, jms:jme )                                    , &
94            intent(out) ::                                               hbot, &
95                                                                         htop
96  real,     dimension( ims:ime, jms:jme )                                    , &
97            intent(in   ) ::                                            xland
98!
99  real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
100             intent(in   )   ::                                           u3d, &
101                                                                          v3d
102  logical,  dimension( ims:ime, jms:jme )                                    , &
103            intent(inout) ::                                      cu_act_flag
104  real, intent(   in) ::                                                 cudt
105  real, intent(   in) ::                                            curr_secs
106  logical, intent(    in) ::                                  adapt_step_flag
107!
108  real,     dimension( ims:ime, jms:jme )                                    , &
109             intent(in   )   ::                                          hpbl, &
110                                                                          hfx, &
111                                                                          qfx
112  integer,   intent(in   )   ::                                    mp_physics
113  integer :: ncloud
114!
115!local
116!
117  real,  dimension( its:ite, jts:jte )  ::                            raincv1, &
118                                                                      raincv2, &
119                                                                      pratec1, &
120                                                                      pratec2
121  real,   dimension( its:ite, kts:kte )  ::                               del, &
122                                                                        prsll, &
123                                                                          dot, &
124                                                                           u1, &
125                                                                           v1, &
126                                                                           t1, &
127                                                                          q1,  &
128                                                                          qc2, &
129                                                                          qi2
130  real,   dimension( its:ite, kts:kte+1 )  ::                           prsii, &
131                                                                          zii
132  real,   dimension( its:ite, kts:kte )  ::                               zll
133  real,   dimension( its:ite)  ::                                         rain
134  real ::                                                          delt,rdelt
135  integer, dimension (its:ite)  ::                                       kbot, &
136                                                                         ktop, &
137                                                                          kuo
138  logical ::  run_param
139  integer ::  i,j,k,kp
140!
141!-------------------------------------------------------------------------------
142! microphysics scheme --> ncloud
143   if (mp_physics .eq. 1) then
144     ncloud = 0
145   elseif ( mp_physics .eq. 2 .or. mp_physics .eq. 3 ) then
146     ncloud = 2
147   elseif ( mp_physics .eq. 5 ) then
148     ncloud = 3
149   elseif ( mp_physics .eq. 4 .or. mp_physics .eq. 14 ) then
150     ncloud = 4
151   elseif ( mp_physics .eq. 9) then
152     ncloud = 6
153   else
154     ncloud = 5
155   endif 
156!
157!-------------------------------------------------------------------------------
158!
159!*** check to see if this is a convection timestep
160!
161      if (adapt_step_flag) then
162        if ( (itimestep .eq. 0) .or. (cudt .eq. 0) .or.                        &
163          (curr_secs + dt >= (int(curr_secs/(cudt*60))+1)*cudt*60)) then
164          run_param = .true.
165        else
166          run_param = .false.
167        endif
168      else
169        if (MOD(itimestep,stepcu) .EQ. 0 .or. itimestep .eq. 0) then
170          run_param = .true.
171        else
172          run_param = .false.
173        endif
174      endif
175!-------------------------------------------------------------------------------
176   if(run_param) then
177      do j=jts,jte
178        do i=its,ite
179          cu_act_flag(i,j)=.TRUE.
180        enddo
181      enddo
182      delt=dt*stepcu
183      rdelt=1./delt
184!
185   do j = jts,jte  !outer most J_loop
186      do k = kts,kte
187        kp = k+1
188        do i = its,ite
189          dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
190          prsll(i,k)=p3d(i,k,j)*0.001
191          prsii(i,k)=p3di(i,k,j)*0.001
192        enddo
193      enddo
194      do i = its,ite
195        prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
196      enddo
197!
198      do i=its,ite
199        zii(i,1)=0.0
200      enddo     
201!
202      do k=kts,kte                                           
203        do i=its,ite
204          zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
205        enddo
206      enddo
207!
208      do k=kts,kte               
209        do i=its,ite                                                 
210          zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
211        enddo                                                         
212      enddo
213!
214      do k=kts,kte
215        do i=its,ite
216          del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
217          u1(i,k)=u3d(i,k,j)
218          v1(i,k)=v3d(i,k,j)
219          t1(i,k)=t3d(i,k,j)
220          q1(i,k)=qv3d(i,k,j)
221          qi2(i,k) = qi3d(i,k,j)
222          qc2(i,k) = qc3d(i,k,j)
223        enddo
224      enddo
225!
226! NCEP SAS
227      call nsas2d(delt=dt,del=del(its,kts),prsl=prsll(its,kts),                &
228              prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts),       &
229              zi=zii(its,kts),ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
230              q1=q1(its,kts),t1=t1(its,kts),rain=rain(its),                    &
231              kbot=kbot(its),ktop=ktop(its),                                   &
232              kuo=kuo(its),                                                    &
233              lat=j,slimsk=xland(ims,j),dot=dot(its,kts),                      &
234              u1=u1(its,kts), v1=v1(its,kts),                                  &
235              cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv,                      &
236              rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2,                               &
237              cice=cice,xls=xls,psat=psat,                                     &
238              ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,               &
239              ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,               &
240              its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
241!
242      do i=its,ite
243        pratec1(i,j)=rain(i)*1000./(stepcu*dt)
244        raincv1(i,j)=rain(i)*1000./(stepcu)
245      enddo
246!
247! NCEP SCV
248      call nscv2d(delt=dt,del=del(its,kts),prsl=prsll(its,kts),                &
249              prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts),       &
250              zi=zii(its,kts),ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
251              q1=q1(its,kts),t1=t1(its,kts),rain=rain(its),                    &
252              kbot=kbot(its),ktop=ktop(its),                                   &
253              kuo=kuo(its),                                                    &
254              slimsk=xland(ims,j),dot=dot(its,kts),                            &
255              u1=u1(its,kts), v1=v1(its,kts),                                  &
256              cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv,                      &
257              rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2,                               &
258              cice=cice,xls=xls,psat=psat,                                     &
259              hpbl=hpbl(ims,j),hfx=hfx(ims,j),qfx=qfx(ims,j),                  &
260              ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,               &
261              ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,               &
262              its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
263!
264   do i=its,ite
265     pratec2(i,j)=rain(i)*1000./(stepcu*dt)
266     raincv2(i,j)=rain(i)*1000./(stepcu)
267   enddo
268!
269   do i=its,ite
270     raincv(i,j) = raincv1(i,j) + raincv2(i,j)
271     pratec(i,j) = pratec1(i,j) + pratec2(i,j)
272     hbot(i,j) = kbot(i)
273     htop(i,j) = ktop(i)
274   enddo
275!
276      IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
277      do k = kts,kte
278        do i= its,ite
279          rthcuten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
280          rqvcuten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
281        enddo
282      enddo
283      ENDIF
284!
285      IF(PRESENT(rucuten).AND.PRESENT(rvcuten)) THEN
286      do k = kts,kte
287        do i= its,ite
288          rucuten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
289          rvcuten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
290        enddo
291      enddo
292      ENDIF
293!
294      if(PRESENT( rqicuten )) THEN
295        IF ( F_QI ) THEN
296        do k=kts,kte
297          do i=its,ite
298            rqicuten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
299          enddo
300        enddo
301        endif
302      endif
303      if(PRESENT( rqccuten )) THEN
304        IF ( F_QC ) THEN
305        do k=kts,kte
306          do i=its,ite
307            rqccuten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
308          enddo
309        enddo
310        endif
311      endif
312!
313   enddo ! outer most J_loop
314  endif
315!
316   end subroutine cu_nsas
317!
318!==============================================================================
319! NCEP SAS (Deep Convection Scheme)
320!==============================================================================
321   subroutine nsas2d(delt,del,prsl,prsi,prslk,zl,zi,                           &
322            ncloud,                                                            &
323            qc2,qi2,                                                           &
324            q1,t1,rain,kbot,ktop,                                              &
325            kuo,                                                               &
326            lat,slimsk,dot,u1,v1,cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2,     &
327            cice,xls,psat,                                                     &
328            ids,ide, jds,jde, kds,kde,                                         &
329            ims,ime, jms,jme, kms,kme,                                         &
330            its,ite, jts,jte, kts,kte)
331!
332!------------------------------------------------------------------------------
333!
334! subprogram:    phys_cps_sas      computes convective heating and moistening
335!                                                      and momentum transport
336!
337! abstract: computes convective heating and moistening using a one
338!   cloud type arakawa-schubert convection scheme originally developed
339!   by georg grell. the scheme has been revised at ncep since initial
340!   implementation in 1993. it includes updraft and downdraft effects.
341!   the closure is the cloud work function. both updraft and downdraft
342!   are assumed to be saturated and the heating and moistening are
343!   accomplished by the compensating environment. convective momentum transport
344!   is taken into account. the name comes from "simplified arakawa-schubert
345!   convection parameterization".
346!
347! developed by hua-lu pan, wan-shu wu, songyou hong, and jongil han
348!   implemented into wrf by kyosun sunny lim and songyou hong
349!   module with cpp-based options is available in grims
350!
351! program history log:
352!   92-03-01  hua-lu pan       operational development
353!   96-03-01  song-you hong    revised closure, and trigger
354!   99-03-01  hua-lu pan       multiple clouds
355!   06-03-01  young-hwa byun   closure based on moisture convergence (optional)
356!   09-10-01  jung-eun kim     f90 format with standard physics modules
357!   10-07-01  jong-il han      revised cloud model,trigger, as in gfs july 2010
358!   10-12-01  kyosun sunny lim wrf compatible version
359!
360!
361! usage:    call phys_cps_sas(delt,delx,del,prsl,prsi,prslk,prsik,zl,zi,       &
362!                             q2,q1,t1,u1,v1,rcs,slimsk,dot,cldwrk,rain,       &
363!                             jcap,ncloud,lat,kbot,ktop,kuo,                   &
364!                             ids,ide, jds,jde, kds,kde,                       &
365!                             ims,ime, jms,jme, kms,kme,                       &
366!                             its,ite, jts,jte, kts,kte)
367!
368!   delt     - real model integration time step
369!   delx     - real model grid interval
370!   del      - real (kms:kme) sigma layer thickness
371!   prsl     - real (ims:ime,kms:kme) pressure values
372!   prsi     - real (ims:ime,kms:kme) pressure values at interface level
373!   prslk    - real (ims:ime,kms:kme) pressure values to the kappa
374!   prsik    - real (ims:ime,kms:kme) pressure values to the kappa at interface lev.
375!   zl       - real (ims:ime,kms:kme) height above sea level
376!   zi       - real (ims:ime,kms:kme) height above sea level at interface level
377!   rcs      - real
378!   slimsk   - real (ims:ime) land(1),sea(0), ice(2) flag
379!   dot      - real (ims:ime,kms:kme) vertical velocity
380!   jcap     - integer spectral truncation
381!   ncloud   - integer number of hydrometeors
382!   lat      - integer  current latitude index
383!
384! output argument list:
385!   q2       - real (ims:ime,kms:kme) detrained hydrometeors in kg/kg
386!            - in case of the  --> qc2(cloud), qi2(ice)
387!   q1       - real (ims:ime,kms:kme) adjusted specific humidity in kg/kg
388!   t1       - real (ims:ime,kms:kme) adjusted temperature in kelvin
389!   u1       - real (ims:ime,kms:kme) adjusted zonal wind in m/s
390!   v1       - real (ims:ime,kms:kme) adjusted meridional wind in m/s
391!   cldwrk   - real (ims:ime) cloud work function
392!   rain     - real (ims:ime) convective rain in meters
393!   kbot     - integer (ims:ime) cloud bottom level
394!   ktop     - integer (ims:ime) cloud top level
395!   kuo      - integer (ims:ime) bit flag indicating deep convection
396!
397! subprograms called:
398!   fpvs     - function to compute saturation vapor pressure
399!
400! remarks: function fpvs is inlined by fpp.
401!          nonstandard automatic arrays are used.
402!
403! references :
404!   pan and wu    (1995, ncep office note)
405!   hong and pan  (1998, mon wea rev)
406!   park and hong (2007,jmsj)
407!   byun and hong (2007, mon wea rev)
408!   han and pan   (2011, wea. forecasting)
409!
410!------------------------------------------------------------------------------
411!------------------------------------------------------------------------------
412   implicit none
413!------------------------------------------------------------------------------
414!
415! model tunable parameters
416!
417   real,parameter  ::  alphal = 0.5,    alphas = 0.5
418   real,parameter  ::  betal  = 0.05,   betas  = 0.05
419   real,parameter  ::  pdpdwn = 0.0,    pdetrn = 200.0
420   real,parameter  ::  c0     = 0.002,  c1     = 0.002
421   real,parameter  ::  pgcon = 0.55
422   real,parameter  ::  xlamdd = 1.0e-4, xlamde = 1.0e-4
423   real,parameter  ::  clam   = 0.1,    cxlamu = 1.0e-4
424   real,parameter  ::  aafac  = 0.1
425   real,parameter  ::  dthk=25.
426   real,parameter  ::  cincrmax = 180.,cincrmin = 120.
427   real,parameter  ::  W1l = -8.E-3
428   real,parameter  ::  W2l = -4.E-2
429   real,parameter  ::  W3l = -5.E-3
430   real,parameter  ::  W4l = -5.E-4
431   real,parameter  ::  W1s = -2.E-4
432   real,parameter  ::  W2s = -2.E-3
433   real,parameter  ::  W3s = -1.E-3
434   real,parameter  ::  W4s = -2.E-5
435   real,parameter  ::  mbdt = 10., edtmaxl = 0.3, edtmaxs = 0.3
436   real,parameter  ::  evfacts = 0.3, evfactl = 0.3
437!
438   real,parameter  ::  tf=233.16,tcr=263.16,tcrf=1.0/(tcr-tf)
439   real,parameter  ::  xk1=2.e-5,xlhor=3.e4,xhver=5000.,theimax=1.
440   real,parameter  ::  xc1=1.e-7,xc2=1.e4,xc3=3.e3,ecesscr=3.0,edtk1=3.e4
441!
442!  passing variables
443!
444   real            ::  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
445   real            ::  pi_,qmin_,t0c_,cice,xlv0,xls,psat
446   integer         ::  lat,                                                    &
447                       ncloud,                                                 &
448                       ids,ide, jds,jde, kds,kde,                              &
449                       ims,ime, jms,jme, kms,kme,                              &
450                       its,ite, jts,jte, kts,kte
451!
452   real            ::  delt,rcs
453   real            ::  del(its:ite,kts:kte),                                   &
454                       prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme),           &
455                       prsi(its:ite,kts:kte+1),                                &
456                       zl(its:ite,kts:kte),zi(its:ite,kts:kte+1),              &
457                       q1(its:ite,kts:kte),t1(its:ite,kts:kte),                &
458                       u1(its:ite,kts:kte),v1(its:ite,kts:kte),                &
459                       dot(its:ite,kts:kte)
460   real            ::  qi2(its:ite,kts:kte)
461   real            ::  qc2(its:ite,kts:kte)
462!
463   real            ::  rain(its:ite)
464   integer         ::  kbot(its:ite),ktop(its:ite),kuo(its:ite)
465   real            ::  slimsk(ims:ime)
466!
467!
468!  local variables and arrays
469!
470   integer         ::  i,k,kmax,kbmax,kbm,jmn,indx,indp,kts1,kte1,kmax1,kk
471   real            ::  p(its:ite,kts:kte),pdot(its:ite),acrtfct(its:ite)
472   real            ::  uo(its:ite,kts:kte),vo(its:ite,kts:kte)
473   real            ::  to(its:ite,kts:kte),qo(its:ite,kts:kte)
474   real            ::  hcko(its:ite,kts:kte)
475   real            ::  qcko(its:ite,kts:kte),eta(its:ite,kts:kte)
476   real            ::  etad(its:ite,kts:kte)
477   real            ::  qrcdo(its:ite,kts:kte)
478   real            ::  pwo(its:ite,kts:kte),pwdo(its:ite,kts:kte)
479   real            ::  dtconv(its:ite)
480   real            ::  deltv(its:ite),acrt(its:ite)
481   real            ::  qeso(its:ite,kts:kte)
482   real            ::  tvo(its:ite,kts:kte),dbyo(its:ite,kts:kte)
483   real            ::  heo(its:ite,kts:kte),heso(its:ite,kts:kte)
484   real            ::  qrcd(its:ite,kts:kte)
485   real            ::  dellah(its:ite,kts:kte),dellaq(its:ite,kts:kte)
486!
487   integer         ::  kb(its:ite),kbcon(its:ite)
488   integer         ::  kbcon1(its:ite)
489   real            ::  hmax(its:ite),delq(its:ite)
490   real            ::  hkbo(its:ite),qkbo(its:ite),pbcdif(its:ite)
491   integer         ::  kbds(its:ite),lmin(its:ite),jmin(its:ite)
492   integer         ::  ktcon(its:ite)
493   integer         ::  ktcon1(its:ite)
494   integer         ::  kbdtr(its:ite),kpbl(its:ite)
495   integer         ::  klcl(its:ite),ktdown(its:ite)
496   real            ::  vmax(its:ite)
497   real            ::  hmin(its:ite),pwavo(its:ite)
498   real            ::  aa1(its:ite),vshear(its:ite)
499   real            ::  qevap(its:ite)
500   real            ::  edt(its:ite)
501   real            ::  edto(its:ite),pwevo(its:ite)
502   real            ::  qcond(its:ite)
503   real            ::  hcdo(its:ite,kts:kte)
504   real            ::  ddp(its:ite),pp2(its:ite)
505   real            ::  qcdo(its:ite,kts:kte)
506   real            ::  adet(its:ite),aatmp(its:ite)
507   real            ::  xhkb(its:ite),xqkb(its:ite)
508   real            ::  xpwav(its:ite),xpwev(its:ite),xhcd(its:ite,kts:kte)
509   real            ::  xaa0(its:ite),f(its:ite),xk(its:ite)
510   real            ::  xmb(its:ite)
511   real            ::  edtx(its:ite),xqcd(its:ite,kts:kte)
512   real            ::  hsbar(its:ite),xmbmax(its:ite)
513   real            ::  xlamb(its:ite,kts:kte),xlamd(its:ite)
514   real            ::  excess(its:ite)
515   real            ::  plcl(its:ite)
516   real            ::  delhbar(its:ite),delqbar(its:ite),deltbar(its:ite)
517   real,save       ::  pcrit(15), acritt(15)
518   real            ::  acrit(15)
519   real            ::  qcirs(its:ite,kts:kte),qrski(its:ite)
520   real            ::  dellal(its:ite,kts:kte)
521   real            ::  rntot(its:ite),delqev(its:ite),delq2(its:ite)
522!
523   real            ::  fent1(its:ite,kts:kte),fent2(its:ite,kts:kte)
524   real            ::  frh(its:ite,kts:kte)
525   real            ::  xlamud(its:ite),sumx(its:ite)
526   real            ::  aa2(its:ite)
527   real            ::  ucko(its:ite,kts:kte),vcko(its:ite,kts:kte)
528   real            ::  ucdo(its:ite,kts:kte),vcdo(its:ite,kts:kte)
529   real            ::  dellau(its:ite,kts:kte),dellav(its:ite,kts:kte)
530   real            ::  delubar(its:ite),delvbar(its:ite)
531   real            ::  qlko_ktcon(its:ite)
532!
533   real            ::  alpha,beta,                                             &
534                       dt2,dtmin,dtmax,dtmaxl,dtmaxs,                          &
535                       el2orc,eps,fact1,fact2,                                 &
536                       tem,tem1,cincr
537   real            ::  dz,dp,es,pprime,qs,                                     &
538                       dqsdp,desdt,dqsdt,gamma,                                &
539                       dt,dq,po,thei,delza,dzfac,                              &
540                       thec,theb,thekb,thekh,theavg,thedif,                    &
541                       omgkb,omgkbp1,omgdif,omgfac,heom,rh,thermal,chi,        &
542                       factor,onemf,dz1,qrch,etah,qlk,qc,rfact,shear,          &
543                       e1,dh,deta,detad,theom,edtmax,dhh,dg,aup,adw,           &
544                       dv1,dv2,dv3,dv1q,dv2q,dv3q,dvq1,                        &
545                       dv1u,dv2u,dv3u,dv1v,dv2v,dv3v,                          &
546                       dellat,xdby,xqrch,xqc,xpw,xpwd,                         &
547                       w1,w2,w3,w4,qrsk,evef,ptem,ptem1
548!
549   logical         ::  totflg, cnvflg(its:ite),flg(its:ite)
550!
551!  climatological critical cloud work functions for closure
552!
553   data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,               &
554              350.,300.,250.,200.,150./
555   data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,                 &
556              .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
557!
558!-----------------------------------------------------------------------
559!
560! define miscellaneous values
561!
562   pi_   = 3.14159
563   qmin_ = 1.0e-30
564   t0c_ = 273.15
565   xlv0 = hvap_
566   rcs  = 1.
567   el2orc = hvap_*hvap_/(rv_*cp_)
568   eps    = rd_/rv_
569   fact1  = (cvap_-cliq_)/rv_
570   fact2  = hvap_/rv_-fact1*t0c_
571   kts1 = kts + 1
572   kte1 = kte - 1
573   dt2    = delt
574   dtmin  = max(dt2,1200.)
575   dtmax  = max(dt2,3600.)
576!
577!
578!  initialize arrays
579!
580   do i = its,ite
581     rain(i)     = 0.0
582     kbot(i)   = kte+1
583     ktop(i)   = 0
584     kuo(i)    = 0
585     cnvflg(i) = .true.
586     dtconv(i) = 3600.
587     pdot(i)   = 0.0
588     edto(i)   = 0.0
589     edtx(i)   = 0.0
590     xmbmax(i) = 0.3
591     excess(i) = 0.0
592     plcl(i)   = 0.0
593     kpbl(i)   = 1
594     aa2(i) = 0.0
595     qlko_ktcon(i) = 0.0
596     pbcdif(i)= 0.0
597     lmin(i) = 1
598     jmin(i) = 1
599     edt(i) = 0.0
600   enddo
601!
602   do k = 1,15
603     acrit(k) = acritt(k) * (975. - pcrit(k))
604   enddo
605!
606! Define top layer for search of the downdraft originating layer
607! and the maximum thetae for updraft
608!
609   kbmax = kte
610   kbm   = kte
611   kmax  = kte
612   do k = kts,kte
613     do i = its,ite
614       if(prsl(i,k).gt.prsi(i,1)*0.45) kbmax = k + 1
615       if(prsl(i,k).gt.prsi(i,1)*0.70) kbm   = k + 1
616       if(prsl(i,k).gt.prsi(i,1)*0.04) kmax  = k + 1
617     enddo
618   enddo
619   kmax = min(kmax,kte)
620   kmax1 = kmax - 1
621   kbm = min(kbm,kte)
622!
623! convert surface pressure to mb from cb
624!
625   do k = kts,kte
626     do i = its,ite
627       pwo(i,k)  = 0.0
628       pwdo(i,k) = 0.0
629       dellal(i,k) = 0.0
630       hcko(i,k) = 0.0
631       qcko(i,k) = 0.0
632       hcdo(i,k) = 0.0
633       qcdo(i,k) = 0.0
634     enddo
635   enddo
636!
637   do k = kts,kmax
638     do i = its,ite
639       p(i,k) = prsl(i,k) * 10.
640       pwo(i,k) = 0.0
641       pwdo(i,k) = 0.0
642       to(i,k) = t1(i,k)
643       qo(i,k) = q1(i,k)
644       dbyo(i,k) = 0.0
645       fent1(i,k) = 1.0
646       fent2(i,k) = 1.0
647       frh(i,k) = 0.0
648       ucko(i,k) = 0.0
649       vcko(i,k) = 0.0
650       ucdo(i,k) = 0.0
651       vcdo(i,k) = 0.0
652       uo(i,k) = u1(i,k) * rcs
653       vo(i,k) = v1(i,k) * rcs
654     enddo
655   enddo
656!
657! column variables
658!
659!  p is pressure of the layer (mb)
660!  t is temperature at t-dt (k)..tn
661!  q is mixing ratio at t-dt (kg/kg)..qn
662!  to is temperature at t+dt (k)... this is after advection and turbulan
663!  qo is mixing ratio at t+dt (kg/kg)..q1
664!
665   do k = kts,kmax
666     do i = its,ite
667       qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
668       qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
669       qeso(i,k) = max(qeso(i,k),qmin_)
670       qo(i,k)   = max(qo(i,k), 1.e-10 )
671       tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
672     enddo
673   enddo
674!
675! compute moist static energy
676!
677   do k = kts,kmax
678     do i = its,ite
679       heo(i,k)  = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qo(i,k)
680       heso(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qeso(i,k)
681     enddo
682   enddo
683!
684! Determine level with largest moist static energy
685! This is the level where updraft starts
686!
687   do i = its,ite
688     hmax(i) = heo(i,1)
689     kb(i) = 1
690   enddo
691!
692   do k = kts1,kbm
693     do i = its,ite
694       if(heo(i,k).gt.hmax(i)) then
695         kb(i) = k
696         hmax(i) = heo(i,k)
697       endif
698     enddo
699   enddo
700!
701   do k = kts,kmax1
702     do i = its,ite
703       if(cnvflg(i)) then
704         dz = .5 * (zl(i,k+1) - zl(i,k))
705         dp = .5 * (p(i,k+1) - p(i,k))
706         es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
707         pprime = p(i,k+1) + (eps-1.) * es
708         qs = eps * es / pprime
709         dqsdp = - qs / pprime
710         desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
711         dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
712         gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
713         dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
714         dq = dqsdt * dt + dqsdp * dp
715         to(i,k) = to(i,k+1) + dt
716         qo(i,k) = qo(i,k+1) + dq
717         po = .5 * (p(i,k) + p(i,k+1))
718         qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
719         qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
720         qeso(i,k) = max(qeso(i,k),qmin_)
721         qo(i,k)   = max(qo(i,k), 1.e-10)
722         frh(i,k)  = 1. - min(qo(i,k)/qeso(i,k), 1.)
723         heo(i,k)  = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
724                cp_ * to(i,k) + hvap_ * qo(i,k)
725         heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
726                cp_ * to(i,k) + hvap_ * qeso(i,k)
727         uo(i,k)   = .5 * (uo(i,k) + uo(i,k+1))
728         vo(i,k)   = .5 * (vo(i,k) + vo(i,k+1))
729       endif
730     enddo
731   enddo
732!
733! look for convective cloud base as the level of free convection
734!
735   do i = its,ite
736     if(cnvflg(i)) then
737       indx = kb(i)
738       hkbo(i) = heo(i,indx)
739       qkbo(i) = qo(i,indx)
740     endif
741   enddo
742!
743   do i = its,ite
744     flg(i) = cnvflg(i)
745     kbcon(i) = kmax
746   enddo
747!
748     do k = kts,kbmax
749       do i = its,ite
750         if(flg(i).and.k.gt.kb(i)) then
751           hsbar(i) = heso(i,k)
752           if(hkbo(i).gt.hsbar(i)) then
753             flg(i) = .false.
754             kbcon(i) = k
755           endif
756         endif
757       enddo
758     enddo
759     do i = its,ite
760       if(kbcon(i).eq.kmax) cnvflg(i) = .false.
761     enddo
762!
763     totflg = .true.
764     do i = its,ite
765       totflg = totflg .and. (.not. cnvflg(i))
766     enddo
767     if(totflg) return
768!
769     do i = its,ite
770       if(cnvflg(i)) then
771!
772!  determine critical convective inhibition
773!  as a function of vertical velocity at cloud base.
774!
775          pdot(i)  = 10.* dot(i,kbcon(i))
776          if(slimsk(i).eq.1.) then
777            w1 = w1l
778            w2 = w2l
779            w3 = w3l
780            w4 = w4l
781          else
782            w1 = w1s
783            w2 = w2s
784            w3 = w3s
785            w4 = w4s
786          endif
787          if(pdot(i).le.w4) then
788            tem = (pdot(i) - w4) / (w3 - w4)
789          elseif(pdot(i).ge.-w4) then
790            tem = - (pdot(i) + w4) / (w4 - w3)
791          else
792            tem = 0.
793          endif
794          tem = max(tem,-1.)
795          tem = min(tem,1.)
796          tem = 1. - tem
797          tem1= .5*(cincrmax-cincrmin)
798          cincr = cincrmax - tem * tem1
799         pbcdif(i) = -p(i,kbcon(i)) + p(i,kb(i))
800         if(pbcdif(i).gt.cincr) cnvflg(i) = .false.
801       endif
802     enddo
803!
804!
805   totflg = .true.
806   do i = its,ite
807     totflg = totflg .and. (.not. cnvflg(i))
808   enddo
809   if(totflg) return
810!
811   do k = kts,kte1
812     do i = its,ite
813       xlamb(i,k) = clam / zi(i,k+1)
814     enddo
815   enddo
816!
817!  assume that updraft entrainment rate above cloud base is
818!  same as that at cloud base
819!
820   do k = kts1,kmax1
821     do i = its,ite
822       if(cnvflg(i).and.(k.gt.kbcon(i))) then
823         xlamb(i,k) = xlamb(i,kbcon(i))
824       endif
825     enddo
826   enddo
827!
828!   assume the detrainment rate for the updrafts to be same as
829!   the entrainment rate at cloud base
830!
831   do i = its,ite
832     if(cnvflg(i)) then
833       xlamud(i) = xlamb(i,kbcon(i))
834     endif
835   enddo
836!
837!  functions rapidly decreasing with height, mimicking a cloud ensemble
838!    (Bechtold et al., 2008)
839!
840   do k = kts1,kmax1
841     do i = its,ite
842       if(cnvflg(i).and.(k.gt.kbcon(i))) then
843         tem = qeso(i,k)/qeso(i,kbcon(i))
844         fent1(i,k) = tem**2
845         fent2(i,k) = tem**3
846       endif
847     enddo
848   enddo
849!
850!  final entrainment rate as the sum of turbulent part and organized entrainment
851!    depending on the environmental relative humidity
852!    (Bechtold et al., 2008)
853!
854   do k = kts1,kmax1
855     do i = its,ite
856       if(cnvflg(i).and.(k.ge.kbcon(i))) then
857          tem = cxlamu * frh(i,k) * fent2(i,k)
858          xlamb(i,k) = xlamb(i,k)*fent1(i,k) + tem
859       endif
860     enddo
861   enddo
862!
863! determine updraft mass flux
864!
865   do k = kts,kte
866     do i = its,ite
867      if(cnvflg(i)) then
868         eta(i,k) = 1.
869       endif
870     enddo
871   enddo
872!
873   do k = kbmax,kts1,-1
874     do i = its,ite
875       if(cnvflg(i).and.k.lt.kbcon(i).and.k.ge.kb(i)) then
876         dz = zi(i,k+2) - zi(i,k+1)
877         ptem     = 0.5*(xlamb(i,k)+xlamb(i,k+1))-xlamud(i)
878         eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
879       endif
880     enddo
881   enddo
882  do k = kts1,kmax1
883     do i = its,ite
884       if(cnvflg(i).and.k.gt.kbcon(i)) then
885         dz  = zi(i,k+1) - zi(i,k)
886         ptem = 0.5*(xlamb(i,k)+xlamb(i,k-1))-xlamud(i)
887         eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
888       endif
889     enddo
890   enddo
891   do i = its,ite
892     if(cnvflg(i)) then
893       dz = zi(i,3) - zi(i,2)
894       ptem     = 0.5*(xlamb(i,1)+xlamb(i,2))-xlamud(i)
895       eta(i,1) = eta(i,2) / (1. + ptem * dz)
896     endif
897   enddo
898!
899! work up updraft cloud properties
900!
901   do i = its,ite
902     if(cnvflg(i)) then
903       indx = kb(i)
904       hcko(i,indx) = hkbo(i)
905       qcko(i,indx) = qkbo(i)
906       ucko(i,indx) = uo(i,indx)
907       vcko(i,indx) = vo(i,indx)
908       pwavo(i) = 0.
909     endif
910   enddo
911!
912! cloud property below cloud base is modified by the entrainment proces
913!
914   do k = kts1,kmax1
915     do i = its,ite
916       if(cnvflg(i).and.k.gt.kb(i)) then
917         dz   = zi(i,k+1) - zi(i,k)
918         tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
919         tem1 = 0.5 * xlamud(i) * dz
920         factor = 1. + tem - tem1
921         ptem = 0.5 * tem + pgcon
922         ptem1= 0.5 * tem - pgcon
923         hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                          &
924                     (heo(i,k)+heo(i,k-1)))/factor
925         ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)                      &
926                     +ptem1*uo(i,k-1))/factor
927         vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)                      &
928                     +ptem1*vo(i,k-1))/factor
929         dbyo(i,k) = hcko(i,k) - heso(i,k)
930       endif
931     enddo
932   enddo
933!
934!   taking account into convection inhibition due to existence of
935!    dry layers below cloud base
936!
937   do i = its,ite
938     flg(i) = cnvflg(i)
939     kbcon1(i) = kmax
940   enddo
941!
942   do k = kts1,kmax
943     do i = its,ite
944       if(flg(i).and.k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
945         kbcon1(i) = k
946         flg(i) = .false.
947       endif
948     enddo
949   enddo
950!
951   do i = its,ite
952     if(cnvflg(i)) then
953       if(kbcon1(i).eq.kmax) cnvflg(i) = .false.
954     endif
955   enddo
956!
957   do i =its,ite
958     if(cnvflg(i)) then
959       tem = p(i,kbcon(i)) - p(i,kbcon1(i))
960       if(tem.gt.dthk) then
961          cnvflg(i) = .false.
962       endif
963     endif
964   enddo
965!
966   totflg = .true.
967   do i = its,ite
968     totflg = totflg .and. (.not. cnvflg(i))
969   enddo
970   if(totflg) return
971!
972!
973!  determine cloud top
974!
975!
976   do i = its,ite
977     flg(i) = cnvflg(i)
978     ktcon(i) = 1
979   enddo
980!
981!   check inversion
982!
983   do k = kts1,kmax1
984     do i = its,ite
985       if(dbyo(i,k).lt.0..and.flg(i).and.k.gt. kbcon1(i)) then
986         ktcon(i) = k
987         flg(i)   = .false.
988       endif
989     enddo
990   enddo
991!
992!
993! check cloud depth
994!
995   do i = its,ite
996     if(cnvflg(i).and.(p(i,kbcon(i)) - p(i,ktcon(i))).lt.150.)                 &
997            cnvflg(i) = .false.
998   enddo
999!
1000   totflg = .true.
1001   do i = its,ite
1002     totflg = totflg .and. (.not. cnvflg(i))
1003   enddo
1004   if(totflg) return
1005!
1006!
1007!  search for downdraft originating level above theta-e minimum
1008!
1009   do i = its,ite
1010     if(cnvflg(i)) then
1011       hmin(i) = heo(i,kbcon1(i))
1012       lmin(i) = kbmax
1013       jmin(i) = kbmax
1014    endif
1015   enddo
1016!
1017   do k = kts1,kbmax
1018     do i = its,ite
1019       if(cnvflg(i).and.k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
1020         lmin(i) = k + 1
1021         hmin(i) = heo(i,k)
1022       endif
1023     enddo
1024   enddo
1025!
1026! make sure that jmin is within the cloud
1027!
1028   do i = its,ite
1029     if(cnvflg(i)) then
1030       jmin(i) = min(lmin(i),ktcon(i)-1)
1031       jmin(i) = max(jmin(i),kbcon1(i)+1)
1032       if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
1033     endif
1034   enddo
1035!
1036!  specify upper limit of mass flux at cloud base
1037!
1038   do i = its,ite
1039     if(cnvflg(i)) then
1040       k = kbcon(i)
1041       dp = 1000. * del(i,k)
1042       xmbmax(i) = dp / (g_ * dt2)
1043     endif
1044   enddo
1045!
1046!
1047! compute cloud moisture property and precipitation
1048!
1049   do i = its,ite
1050     aa1(i) = 0.
1051   enddo
1052!
1053   do k = kts1,kmax
1054     do i = its,ite
1055       if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1056         dz = .5 * (zl(i,k+1) - zl(i,k-1))
1057         dz1 = (zi(i,k+1) - zi(i,k))
1058         gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1059         qrch = qeso(i,k)                                                      &
1060              + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1061         tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz1
1062         tem1 = 0.5 * xlamud(i) * dz1
1063         factor = 1. + tem - tem1
1064         qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                           &
1065                    (qo(i,k)+qo(i,k-1)))/factor
1066         qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1067!
1068! check if there is excess moisture to release latent heat
1069!
1070         if(qcirs(i,k).gt.0. .and. k.ge.kbcon(i)) then
1071           etah = .5 * (eta(i,k) + eta(i,k-1))
1072           if(ncloud.gt.0..and.k.gt.jmin(i)) then
1073             dp = 1000. * del(i,k)
1074             qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz1)
1075             dellal(i,k) = etah * c1 * dz1 * qlk * g_ / dp
1076           else
1077             qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz1)
1078           endif
1079           aa1(i) = aa1(i) - dz1 * g_ * qlk
1080           qc = qlk + qrch
1081           pwo(i,k) = etah * c0 * dz1 * qlk
1082           qcko(i,k) = qc
1083           pwavo(i) = pwavo(i) + pwo(i,k)
1084         endif
1085       endif
1086     enddo
1087   enddo
1088!
1089! calculate cloud work function at t+dt
1090!
1091   do k = kts1,kmax
1092     do i = its,ite
1093       if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1094         dz1 = zl(i,k+1) - zl(i,k)
1095         gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1096         rfact =  1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1097         aa1(i) = aa1(i) +dz1 * (g_ / (cp_ * to(i,k)))                         &
1098                  * dbyo(i,k) / (1. + gamma)* rfact
1099         aa1(i) = aa1(i)+dz1 * g_ * fv_ *                                      &
1100                  max(0.,(qeso(i,k) - qo(i,k)))
1101       endif
1102     enddo
1103   enddo
1104!
1105   do i = its,ite
1106     if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1107   enddo
1108!
1109      totflg = .true.
1110      do i=its,ite
1111        totflg = totflg .and. (.not. cnvflg(i))
1112      enddo
1113      if(totflg) return
1114!
1115!    estimate the convective overshooting as the level
1116!    where the [aafac * cloud work function] becomes zero,
1117!    which is the final cloud top
1118!
1119      do i = its,ite
1120        if (cnvflg(i)) then
1121          aa2(i) = aafac * aa1(i)
1122        endif
1123      enddo
1124!
1125      do i = its,ite
1126        flg(i) = cnvflg(i)
1127        ktcon1(i) = kmax1
1128      enddo
1129!
1130      do k = kts1, kmax
1131        do i = its, ite
1132          if (flg(i)) then
1133            if(k.ge.ktcon(i).and.k.lt.kmax) then
1134              dz1 = zl(i,k+1) - zl(i,k)
1135              gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1136              rfact =  1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1137              aa2(i) = aa2(i) +dz1 * (g_ / (cp_ * to(i,k)))                    &
1138                       * dbyo(i,k) / (1. + gamma)* rfact
1139              if(aa2(i).lt.0.) then
1140                ktcon1(i) = k
1141                flg(i) = .false.
1142              endif
1143            endif
1144          endif
1145        enddo
1146      enddo
1147!
1148!  compute cloud moisture property, detraining cloud water
1149!  and precipitation in overshooting layers
1150!
1151   do k = kts1,kmax
1152     do i = its,ite
1153       if (cnvflg(i)) then
1154         if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1155           dz = (zi(i,k+1) - zi(i,k))
1156           gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1157           qrch = qeso(i,k)+ gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1158           tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1159           tem1 = 0.5 * xlamud(i) * dz
1160           factor = 1. + tem - tem1
1161           qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                         &
1162                      (qo(i,k)+qo(i,k-1)))/factor
1163           qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1164!
1165!  check if there is excess moisture to release latent heat
1166!
1167           if(qcirs(i,k).gt.0.) then
1168             etah = .5 * (eta(i,k) + eta(i,k-1))
1169             if(ncloud.gt.0.) then
1170               dp = 1000. * del(i,k)
1171               qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz)
1172               dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
1173             else
1174               qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz)
1175             endif
1176             qc = qlk + qrch
1177             pwo(i,k) = etah * c0 * dz * qlk
1178             qcko(i,k) = qc
1179             pwavo(i) = pwavo(i) + pwo(i,k)
1180           endif
1181         endif
1182       endif
1183     enddo
1184   enddo
1185!
1186! exchange ktcon with ktcon1
1187!
1188   do i = its,ite
1189     if(cnvflg(i)) then
1190       kk = ktcon(i)
1191       ktcon(i) = ktcon1(i)
1192       ktcon1(i) = kk
1193     endif
1194   enddo
1195!
1196! this section is ready for cloud water
1197!
1198   if (ncloud.gt.0) then
1199!
1200!  compute liquid and vapor separation at cloud top
1201!
1202   do i = its,ite
1203     if(cnvflg(i)) then
1204     k = ktcon(i)-1
1205       gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1206       qrch = qeso(i,k)                                                     &
1207              + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1208       dq = qcko(i,k) - qrch
1209!
1210!  check if there is excess moisture to release latent heat
1211!
1212       if(dq.gt.0.) then
1213         qlko_ktcon(i) = dq
1214         qcko(i,k) = qrch
1215       endif
1216     endif
1217   enddo
1218   endif
1219!
1220! ..... downdraft calculations .....
1221!
1222! determine downdraft strength in terms of wind shear
1223!
1224   do i = its,ite
1225     if(cnvflg(i)) then
1226       vshear(i) = 0.
1227     endif
1228   enddo
1229!
1230   do k = kts1,kmax
1231     do i = its,ite
1232       if(k.gt.kb(i).and.k.le.ktcon(i).and.cnvflg(i)) then
1233         shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2                                  &
1234                       + (vo(i,k)-vo(i,k-1)) ** 2)
1235         vshear(i) = vshear(i) + shear
1236       endif
1237     enddo
1238   enddo
1239!
1240   do i = its,ite
1241     if(cnvflg(i)) then
1242       vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1243       e1 = 1.591-.639*vshear(i)                                               &
1244           +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1245       edt(i)  = 1.-e1
1246       edt(i)  = min(edt(i),.9)
1247       edt(i)  = max(edt(i),.0)
1248       edto(i) = edt(i)
1249       edtx(i) = edt(i)
1250     endif
1251   enddo
1252!
1253! determine detrainment rate between 1 and kbdtr
1254!
1255   do i = its,ite
1256     if(cnvflg(i)) then
1257       sumx(i) = 0.
1258     endif
1259   enddo
1260!
1261   do k = kts,kmax1
1262     do i = its,ite
1263       if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1264         dz = zi(i,k+2) - zi(i,k+1)
1265         sumx(i) = sumx(i) + dz
1266       endif
1267     enddo
1268   enddo
1269!
1270   do i = its,ite
1271     kbdtr(i) = kbcon(i)
1272     beta = betas
1273     if(slimsk(i).eq.1.) beta = betal
1274     if(cnvflg(i)) then
1275       kbdtr(i) = kbcon(i)
1276       kbdtr(i) = max(kbdtr(i),1)
1277       dz =(sumx(i)+zi(i,2))/float(kbcon(i))
1278       tem = 1./float(kbcon(i))
1279       xlamd(i) = (1.-beta**tem)/dz
1280     endif
1281   enddo
1282!
1283! determine downdraft mass flux
1284!
1285   do k = kts,kmax
1286     do i = its,ite
1287       if(cnvflg(i)) then
1288         etad(i,k) = 1.
1289       endif
1290       qrcdo(i,k) = 0.
1291       qrcd(i,k) = 0.
1292     enddo
1293   enddo
1294!
1295   do k = kmax1,kts,-1
1296     do i = its,ite
1297       if(cnvflg(i)) then
1298         if(k.lt.jmin(i).and.k.ge.kbcon(i))then
1299           dz = (zi(i,k+2) - zi(i,k+1))
1300           ptem = xlamdd-xlamde
1301           etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1302         elseif(k.lt.kbcon(i))then
1303           dz = (zi(i,k+2) - zi(i,k+1))
1304           ptem = xlamd(i)+xlamdd-xlamde
1305           etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1306         endif
1307       endif
1308     enddo
1309   enddo
1310!
1311!
1312! downdraft moisture properties
1313!
1314   do i = its,ite
1315     if(cnvflg(i)) then
1316      pwevo(i) = 0.
1317     endif
1318   enddo
1319!
1320   do i = its,ite
1321     if(cnvflg(i))  then
1322       jmn = jmin(i)
1323       hcdo(i,jmn) = heo(i,jmn)
1324       qcdo(i,jmn) = qo(i,jmn)
1325       qrcdo(i,jmn) = qeso(i,jmn)
1326       ucdo(i,jmn) = uo(i,jmn)
1327       vcdo(i,jmn) = vo(i,jmn)
1328     endif
1329   enddo
1330!
1331   do k = kmax1,kts,-1
1332     do i = its,ite
1333       if (cnvflg(i) .and. k.lt.jmin(i)) then
1334         dz = zi(i,k+2) - zi(i,k+1)
1335         if(k.ge.kbcon(i)) then
1336           tem  = xlamde * dz
1337           tem1 = 0.5 * xlamdd * dz
1338         else
1339           tem  = xlamde * dz
1340           tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1341          endif
1342          factor = 1. + tem - tem1
1343          ptem = 0.5 * tem - pgcon
1344          ptem1= 0.5 * tem + pgcon
1345          hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*                     &
1346                      (heo(i,k)+heo(i,k+1)))/factor
1347          ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)               &
1348                     +ptem1*uo(i,k))/factor
1349          vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)               &
1350                     +ptem1*vo(i,k))/factor
1351          dbyo(i,k) = hcdo(i,k) - heso(i,k)
1352       endif
1353     enddo
1354   enddo
1355!
1356   do k = kmax1,kts,-1
1357     do i = its,ite
1358       if(cnvflg(i).and.k.lt.jmin(i)) then
1359         dq = qeso(i,k)
1360         dt = to(i,k)
1361         gamma = el2orc * dq / dt**2
1362         qrcdo(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dbyo(i,k)
1363         detad = etad(i,k+1) - etad(i,k)
1364         dz = zi(i,k+2) - zi(i,k+1)
1365         if(k.ge.kbcon(i)) then
1366            tem  = xlamde * dz
1367            tem1 = 0.5 * xlamdd * dz
1368         else
1369            tem  = xlamde * dz
1370            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1371         endif
1372         factor = 1. + tem - tem1
1373         qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5*                      &
1374                     (qo(i,k)+qo(i,k+1)))/factor
1375         pwdo(i,k) = etad(i,k+1) * qcdo(i,k) -etad(i,k+1) * qrcdo(i,k)
1376         qcdo(i,k) = qrcdo(i,k)
1377         pwevo(i) = pwevo(i) + pwdo(i,k)
1378       endif
1379     enddo
1380   enddo
1381!
1382! final downdraft strength dependent on precip
1383! efficiency (edt), normalized condensate (pwav), and
1384! evaporate (pwev)
1385!
1386   do i = its,ite
1387     edtmax = edtmaxl
1388     if(slimsk(i).eq.2.) edtmax = edtmaxs
1389     if(cnvflg(i)) then
1390       if(pwevo(i).lt.0.) then
1391         edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1392         edto(i) = min(edto(i),edtmax)
1393       else
1394         edto(i) = 0.
1395       endif
1396     endif
1397   enddo
1398!
1399! downdraft cloudwork functions
1400!
1401   do k = kmax1,kts,-1
1402     do i = its,ite
1403       if(cnvflg(i).and.k.lt.jmin(i)) then
1404         gamma = el2orc * qeso(i,k) / to(i,k)**2
1405         dhh=hcdo(i,k)
1406         dt=to(i,k)
1407         dg=gamma
1408         dh=heso(i,k)
1409         dz=-1.*(zl(i,k+1)-zl(i,k))
1410         aa1(i)=aa1(i)+edto(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))             &
1411                *(1.+fv_*cp_*dg*dt/hvap_)
1412         aa1(i)=aa1(i)+edto(i)*dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1413       endif
1414     enddo
1415   enddo
1416!
1417   do i = its,ite
1418     if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1419   enddo
1420!
1421   totflg = .true.
1422   do i=its,ite
1423     totflg = totflg .and. (.not. cnvflg(i))
1424   enddo
1425   if(totflg) return
1426!
1427! what would the change be, that a cloud with unit mass
1428! will do to the environment?
1429!
1430   do k = kts,kmax
1431     do i = its,ite
1432       if(cnvflg(i)) then
1433         dellah(i,k) = 0.
1434         dellaq(i,k) = 0.
1435         dellau(i,k) = 0.
1436         dellav(i,k) = 0.
1437       endif
1438     enddo
1439   enddo
1440!
1441   do i = its,ite
1442     if(cnvflg(i)) then
1443       dp = 1000. * del(i,1)
1444       dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)                          &
1445                   - heo(i,1)) * g_ / dp
1446       dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1)                          &
1447                   - qo(i,1)) * g_ / dp
1448       dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)                          &
1449                   - uo(i,1)) * g_ / dp
1450       dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)                          &
1451                   - vo(i,1)) * g_ / dp
1452     endif
1453   enddo
1454!
1455! changed due to subsidence and entrainment
1456!
1457   do k = kts1,kmax1
1458     do i = its,ite
1459       if(cnvflg(i).and.k.lt.ktcon(i)) then
1460         aup = 1.
1461         if(k.le.kb(i)) aup = 0.
1462         adw = 1.
1463         if(k.gt.jmin(i)) adw = 0.
1464         dv1= heo(i,k)
1465         dv2 = .5 * (heo(i,k) + heo(i,k-1))
1466         dv3= heo(i,k-1)
1467         dv1q= qo(i,k)
1468         dv2q = .5 * (qo(i,k) + qo(i,k-1))
1469         dv3q= qo(i,k-1)
1470         dv1u = uo(i,k)
1471         dv2u = .5 * (uo(i,k) + uo(i,k-1))
1472         dv3u = uo(i,k-1)
1473         dv1v = vo(i,k)
1474         dv2v = .5 * (vo(i,k) + vo(i,k-1))
1475         dv3v = vo(i,k-1)
1476         dp = 1000. * del(i,k)
1477         dz = zi(i,k+1) - zi(i,k)
1478         tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1))
1479         tem1 = xlamud(i)
1480         if(k.le.kbcon(i)) then
1481           ptem  = xlamde
1482           ptem1 = xlamd(i)+xlamdd
1483         else
1484           ptem  = xlamde
1485           ptem1 = xlamdd
1486         endif
1487         deta = eta(i,k) - eta(i,k-1)
1488         detad = etad(i,k) - etad(i,k-1)
1489         dellah(i,k) = dellah(i,k) +                                           &
1490             ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1               &
1491         - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3               &
1492         - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2*dz              &
1493         +  aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz                  &
1494         +  adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz) *g_/dp
1495         dellaq(i,k) = dellaq(i,k) +                                           &
1496             ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1q              &
1497         - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3q              &
1498         - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz             &
1499         +  aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz                  &
1500         +  adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz) *g_/dp
1501         dellau(i,k) = dellau(i,k) +                                           &
1502             ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1u              &
1503         - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3u              &
1504         - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz             &
1505         +  aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz                  &
1506         +  adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz          &
1507         -  pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u))*g_/dp
1508!
1509         dellav(i,k) = dellav(i,k) +                                           &
1510             ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1v              &
1511         - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3v              &
1512         - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz             &
1513         +  aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz                  &
1514         +  adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz          &
1515         -  pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v))*g_/dp
1516       endif
1517     enddo
1518   enddo
1519!
1520! cloud top
1521!
1522   do i = its,ite
1523     if(cnvflg(i)) then
1524       indx = ktcon(i)
1525       dp = 1000. * del(i,indx)
1526       dv1 = heo(i,indx-1)
1527       dellah(i,indx) = eta(i,indx-1) *                                        &
1528                        (hcko(i,indx-1) - dv1) * g_ / dp
1529       dvq1 = qo(i,indx-1)
1530       dellaq(i,indx) = eta(i,indx-1) *                                        &
1531                        (qcko(i,indx-1) - dvq1) * g_ / dp
1532       dv1u = uo(i,indx-1)
1533       dellau(i,indx) = eta(i,indx-1) *                                        &
1534                        (ucko(i,indx-1) - dv1u) * g_ / dp
1535       dv1v = vo(i,indx-1)
1536       dellav(i,indx) = eta(i,indx-1) *                                        &
1537                        (vcko(i,indx-1) - dv1v) * g_ / dp
1538!
1539!  cloud water
1540!
1541       dellal(i,indx) = eta(i,indx-1) * qlko_ktcon(i) * g_ / dp
1542     endif
1543   enddo
1544!
1545! final changed variable per unit mass flux
1546!
1547   do k = kts,kmax
1548     do i = its,ite
1549       if(cnvflg(i).and.k.gt.ktcon(i)) then
1550         qo(i,k) = q1(i,k)
1551         to(i,k) = t1(i,k)
1552       endif
1553       if(cnvflg(i).and.k.le.ktcon(i)) then
1554         qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1555         dellat  = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1556         to(i,k) = dellat * mbdt + t1(i,k)
1557         qo(i,k) = max(qo(i,k),1.0e-10)
1558       endif
1559     enddo
1560   enddo
1561!
1562!------------------------------------------------------------------------------
1563!
1564! the above changed environment is now used to calulate the
1565! effect the arbitrary cloud (with unit mass flux)
1566! which then is used to calculate the real mass flux,
1567! necessary to keep this change in balance with the large-scale
1568! destabilization.
1569!
1570! environmental conditions again
1571!
1572!------------------------------------------------------------------------------
1573!
1574   do k = kts,kmax
1575     do i = its,ite
1576       if(cnvflg(i)) then
1577         qeso(i,k)=0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1578         qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1579         qeso(i,k) = max(qeso(i,k),qmin_)
1580         tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
1581       endif
1582     enddo
1583   enddo
1584!
1585   do i = its,ite
1586     if(cnvflg(i)) then
1587       xaa0(i) = 0.
1588       xpwav(i) = 0.
1589     endif
1590   enddo
1591!
1592! moist static energy
1593!
1594   do k = kts,kmax1
1595     do i = its,ite
1596       if(cnvflg(i)) then
1597         dz = .5 * (zl(i,k+1) - zl(i,k))
1598         dp = .5 * (p(i,k+1) - p(i,k))
1599         es =0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1600         pprime = p(i,k+1) + (eps-1.) * es
1601         qs = eps * es / pprime
1602         dqsdp = - qs / pprime
1603         desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1604         dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
1605         gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1606         dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
1607         dq = dqsdt * dt + dqsdp * dp
1608         to(i,k) = to(i,k+1) + dt
1609         qo(i,k) = qo(i,k+1) + dq
1610         po = .5 * (p(i,k) + p(i,k+1))
1611         qeso(i,k) =0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1612         qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
1613         qeso(i,k) = max(qeso(i,k),qmin_)
1614         qo(i,k)   = max(qo(i,k), 1.0e-10)
1615         heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
1616                     cp_ * to(i,k) + hvap_ * qo(i,k)
1617         heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                         &
1618                     cp_ * to(i,k) + hvap_ * qeso(i,k)
1619       endif
1620     enddo
1621   enddo
1622!
1623   k = kmax
1624   do i = its,ite
1625     if(cnvflg(i)) then
1626       heo(i,k)  = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qo(i,k)
1627       heso(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qeso(i,k)
1628     endif
1629   enddo
1630!
1631   do i = its,ite
1632     if(cnvflg(i)) then
1633       xaa0(i) = 0.
1634       xpwav(i) = 0.
1635       indx = kb(i)
1636       xhkb(i) = heo(i,indx)
1637       xqkb(i) = qo(i,indx)
1638       hcko(i,indx) = xhkb(i)
1639       qcko(i,indx) = xqkb(i)
1640     endif
1641   enddo
1642!
1643! ..... static control .....
1644!
1645! moisture and cloud work functions
1646!
1647   do k = kts1,kmax1
1648     do i = its,ite
1649       if(cnvflg(i).and.k.gt.kb(i).and.k.le.ktcon(i)) then
1650         dz = zi(i,k+1) - zi(i,k)
1651         tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1652         tem1 = 0.5 * xlamud(i) * dz
1653         factor = 1. + tem - tem1
1654         hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                         &
1655                    (heo(i,k)+heo(i,k-1)))/factor
1656       endif
1657     enddo
1658   enddo
1659!
1660   do k = kts1,kmax1
1661     do i = its,ite
1662       if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1663         dz = zi(i,k+1) - zi(i,k)
1664         gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1665         xdby = hcko(i,k) - heso(i,k)
1666         xqrch = qeso(i,k)                                                     &
1667              + gamma * xdby / (hvap_ * (1. + gamma))
1668         tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1669         tem1 = 0.5 * xlamud(i) * dz
1670         factor = 1. + tem - tem1
1671         qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*(qo(i,k)+qo(i,k-1)))/factor
1672         dq = eta(i,k) * qcko(i,k) - eta(i,k) * xqrch
1673         if(k.ge.kbcon(i).and.dq.gt.0.) then
1674           etah = .5 * (eta(i,k) + eta(i,k-1))
1675           if(ncloud.gt.0..and.k.gt.jmin(i)) then
1676             qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1677           else
1678             qlk = dq / (eta(i,k) + etah * c0 * dz)
1679           endif
1680           if(k.lt.ktcon1(i)) then
1681             xaa0(i) = xaa0(i) - dz * g_ * qlk
1682           endif
1683           qcko(i,k) = qlk + xqrch
1684           xpw = etah * c0 * dz * qlk
1685           xpwav(i) = xpwav(i) + xpw
1686         endif
1687       endif
1688       if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1689         dz1 = zl(i,k+1) - zl(i,k)
1690         gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1691         rfact =  1. + fv_ * cp_ * gamma                                       &
1692                  * to(i,k) / hvap_
1693         xdby = hcko(i,k) - heso(i,k)
1694         xaa0(i) = xaa0(i)                                                     &
1695                 + dz1 * (g_ / (cp_ * to(i,k)))                                &
1696                 * xdby / (1. + gamma)                                         &
1697                 * rfact
1698         xaa0(i)=xaa0(i)+                                                      &
1699                  dz1 * g_ * fv_ *                                             &
1700                  max(0.,(qeso(i,k) - qo(i,k)))
1701       endif
1702     enddo
1703   enddo
1704!
1705! ..... downdraft calculations .....
1706!
1707!
1708! downdraft moisture properties
1709!
1710   do i = its,ite
1711     xpwev(i) = 0.
1712   enddo
1713   do i = its,ite
1714     if(cnvflg(i)) then
1715       jmn = jmin(i)
1716       xhcd(i,jmn) = heo(i,jmn)
1717       xqcd(i,jmn) = qo(i,jmn)
1718       qrcd(i,jmn) = qeso(i,jmn)
1719     endif
1720   enddo
1721!
1722   do k = kmax1,kts, -1
1723     do i = its,ite
1724       if(cnvflg(i).and.k.lt.jmin(i)) then
1725         dz = zi(i,k+2) - zi(i,k+1)
1726         if(k.ge.kbcon(i)) then
1727            tem  = xlamde * dz
1728            tem1 = 0.5 * xlamdd * dz
1729         else
1730            tem  = xlamde * dz
1731            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1732         endif
1733         factor = 1. + tem - tem1
1734         xhcd(i,k) = ((1.-tem1)*xhcd(i,k+1)+tem*0.5*                        &
1735                    (heo(i,k)+heo(i,k+1)))/factor
1736       endif
1737     enddo
1738   enddo
1739!
1740   do k = kmax1,kts, -1
1741     do i = its,ite
1742       if(cnvflg(i).and.k.lt.jmin(i)) then
1743         dq = qeso(i,k)
1744         dt = to(i,k)
1745         gamma = el2orc * dq / dt**2
1746         dh = xhcd(i,k) - heso(i,k)
1747         qrcd(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dh
1748         dz = zi(i,k+2) - zi(i,k+1)
1749         if(k.ge.kbcon(i)) then
1750           tem  = xlamde * dz
1751           tem1 = 0.5 * xlamdd * dz
1752         else
1753           tem  = xlamde * dz
1754           tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1755         endif
1756         factor = 1. + tem - tem1
1757         xqcd(i,k) = ((1.-tem1)*xqcd(i,k+1)+tem*0.5*                           &
1758                   (qo(i,k)+qo(i,k+1)))/factor
1759         xpwd     = etad(i,k+1) * (xqcd(i,k) - qrcd(i,k))
1760         xqcd(i,k)= qrcd(i,k)
1761         xpwev(i) = xpwev(i) + xpwd
1762       endif
1763     enddo
1764   enddo
1765!
1766   do i = its,ite
1767     edtmax = edtmaxl
1768     if(slimsk(i).eq.2.) edtmax = edtmaxs
1769     if(cnvflg(i)) then
1770       if(xpwev(i).ge.0.) then
1771         edtx(i) = 0.
1772       else
1773         edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1774         edtx(i) = min(edtx(i),edtmax)
1775       endif
1776     endif
1777   enddo
1778!
1779! downdraft cloudwork functions
1780!
1781   do k = kmax1,kts, -1
1782     do i = its,ite
1783       if(cnvflg(i).and.k.lt.jmin(i)) then
1784         gamma = el2orc * qeso(i,k) / to(i,k)**2
1785         dhh=xhcd(i,k)
1786         dt= to(i,k)
1787         dg= gamma
1788         dh= heso(i,k)
1789         dz=-1.*(zl(i,k+1)-zl(i,k))
1790         xaa0(i)=xaa0(i)+edtx(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))           &
1791                 *(1.+fv_*cp_*dg*dt/hvap_)
1792         xaa0(i)=xaa0(i)+edtx(i)*                                              &
1793                  dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1794       endif
1795     enddo
1796   enddo
1797!
1798! calculate critical cloud work function
1799!
1800   do i = its,ite
1801     acrt(i) = 0.
1802     if(cnvflg(i)) then
1803       if(p(i,ktcon(i)).lt.pcrit(15))then
1804         acrt(i)=acrit(15)*(975.-p(i,ktcon(i)))/(975.-pcrit(15))
1805       else if(p(i,ktcon(i)).gt.pcrit(1))then
1806         acrt(i)=acrit(1)
1807       else
1808         k = int((850. - p(i,ktcon(i)))/50.) + 2
1809         k = min(k,15)
1810         k = max(k,2)
1811         acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*                               &
1812              (p(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1813        endif
1814      endif
1815    enddo
1816!
1817   do i = its,ite
1818     acrtfct(i) = 1.
1819     w1 = w1s
1820     w2 = w2s
1821     w3 = w3s
1822     w4 = w4s
1823     if(slimsk(i).eq.1.) then
1824       w1 = w1l
1825       w2 = w2l
1826       w3 = w3l
1827       w4 = w4l
1828     endif
1829     if(cnvflg(i)) then
1830       if(pdot(i).le.w4) then
1831         acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1832       elseif(pdot(i).ge.-w4) then
1833       acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1834       else
1835         acrtfct(i) = 0.
1836       endif
1837       acrtfct(i) = max(acrtfct(i),-1.)
1838       acrtfct(i) = min(acrtfct(i),1.)
1839       acrtfct(i) = 1. - acrtfct(i)
1840       dtconv(i) = dt2 + max((1800. - dt2),0.) * (pdot(i) - w2) / (w1 - w2)   
1841       dtconv(i) = max(dtconv(i),dtmin)
1842       dtconv(i) = min(dtconv(i),dtmax)
1843!
1844     endif
1845   enddo
1846!
1847! large scale forcing
1848!
1849   do i= its,ite
1850     if(cnvflg(i)) then
1851       f(i) = (aa1(i) - acrt(i) * acrtfct(i)) / dtconv(i)
1852       if(f(i).le.0.) cnvflg(i) = .false.
1853     endif
1854     if(cnvflg(i)) then
1855       xk(i) = (xaa0(i) - aa1(i)) / mbdt
1856       if(xk(i).ge.0.) cnvflg(i) = .false.
1857     endif
1858!
1859! kernel, cloud base mass flux
1860!
1861     if(cnvflg(i)) then
1862       xmb(i) = -f(i) / xk(i)
1863       xmb(i) = min(xmb(i),xmbmax(i))
1864     endif
1865!
1866     if(cnvflg(i)) then
1867     endif
1868!
1869   enddo
1870   totflg = .true.
1871   do i = its,ite
1872     totflg = totflg .and. (.not. cnvflg(i))
1873   enddo
1874   if(totflg) return
1875!
1876! restore t0 and qo to t1 and q1 in case convection stops
1877!
1878   do k = kts,kmax
1879     do i = its,ite
1880       if (cnvflg(i)) then
1881       to(i,k) = t1(i,k)
1882       qo(i,k) = q1(i,k)
1883       uo(i,k) = u1(i,k)
1884       vo(i,k) = v1(i,k)
1885       qeso(i,k) = 0.01*fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1886       qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1887       qeso(i,k) = max(qeso(i,k),qmin_)
1888       endif
1889     enddo
1890   enddo
1891!
1892! feedback: simply the changes from the cloud with unit mass flux
1893!           multiplied by  the mass flux necessary to keep the
1894!           equilibrium with the larger-scale.
1895!
1896   do i = its,ite
1897     delhbar(i) = 0.
1898     delqbar(i) = 0.
1899     deltbar(i) = 0.
1900     qcond(i) = 0.
1901     qrski(i) = 0.
1902     delubar(i) = 0.
1903     delvbar(i) = 0.
1904   enddo
1905!
1906   do k = kts,kmax
1907     do i = its,ite
1908       if(cnvflg(i).and.k.le.ktcon(i)) then
1909         aup = 1.
1910         if(k.le.kb(i)) aup = 0.
1911         adw = 1.
1912         if(k.gt.jmin(i)) adw = 0.
1913         dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1914         t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1915         q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1916         tem=1./rcs
1917         u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1918         v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1919         dp = 1000. * del(i,k)
1920         delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
1921         delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
1922         deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
1923         delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
1924         delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
1925       endif
1926     enddo
1927   enddo
1928!
1929   do i = its,ite
1930     if(cnvflg(i)) then
1931     endif
1932   enddo
1933!
1934   do k = kts,kmax
1935     do i = its,ite
1936       if (cnvflg(i) .and. k.le.ktcon(i)) then
1937         qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1938         qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
1939         qeso(i,k) = max(qeso(i,k), qmin_ )
1940       endif
1941     enddo
1942   enddo
1943!
1944   do i = its,ite
1945     rntot(i) = 0.
1946     delqev(i) = 0.
1947     delq2(i) = 0.
1948     flg(i) = cnvflg(i)
1949   enddo
1950!
1951!  comptute rainfall 
1952!
1953   do k = kmax,kts,-1
1954     do i = its,ite
1955       if(cnvflg(i).and.k.lt.ktcon(i)) then
1956         aup = 1.
1957         if(k.le.kb(i)) aup = 0.
1958         adw = 1.
1959         if(k.ge.jmin(i)) adw = 0.
1960         rntot(i) = rntot(i)                                                   &
1961               + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
1962               * xmb(i) * .001 * dt2
1963       endif
1964     enddo
1965   enddo
1966!
1967!  conversion rainfall (m) and compute the evaporation of falling raindrops
1968!
1969   do k = kmax,kts,-1
1970     do i = its,ite
1971       delq(i) = 0.0
1972       deltv(i) = 0.0
1973       qevap(i) = 0.0
1974       if(cnvflg(i).and.k.lt.ktcon(i)) then
1975         aup = 1.
1976         if(k.le.kb(i)) aup = 0.
1977         adw = 1.
1978         if(k.ge.jmin(i)) adw = 0.
1979         rain(i) = rain(i)                                                     &
1980               + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
1981               * xmb(i) * .001 * dt2
1982       endif
1983       if(flg(i).and.k.lt.ktcon(i)) then
1984         evef = edt(i) * evfacts
1985         if(slimsk(i).eq.1.) evef = edt(i) * evfactl
1986         qcond(i) = evef * (q1(i,k) - qeso(i,k)) / (1. + el2orc *              &
1987                  qeso(i,k) / t1(i,k)**2)
1988         dp = 1000. * del(i,k)
1989         if(rain(i).gt.0..and.qcond(i).lt.0.) then
1990           qevap(i) = -qcond(i) * (1. - exp(-.32 * sqrt(dt2 * rain(i))))
1991           qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
1992           delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
1993           if (delq2(i).gt.rntot(i)) then
1994             qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
1995             flg(i) = .false.
1996           endif
1997         endif
1998         if(rain(i).gt.0..and.qevap(i).gt.0.) then
1999           q1(i,k) = q1(i,k) + qevap(i)
2000           t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
2001           rain(i) = rain(i) - .001 * qevap(i) * dp / g_
2002           delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
2003           deltv(i) =  - (hvap_/cp_)*qevap(i)/dt2
2004           delq(i) =  + qevap(i)/dt2
2005         endif
2006         dellaq(i,k) = dellaq(i,k) + delq(i)/xmb(i)
2007         delqbar(i)  = delqbar(i) + delq(i)*dp/g_
2008         deltbar(i)  = deltbar(i) + deltv(i)*dp/g_
2009       endif
2010     enddo
2011   enddo
2012!
2013!
2014! consider the negative rain in the event of rain evaporation and downdrafts
2015!
2016   do i = its,ite
2017     if(cnvflg(i)) then
2018       if(rain(i).lt.0..and..not.flg(i)) rain(i) = 0.
2019       if(rain(i).le.0.) then
2020         rain(i) = 0.
2021       else
2022         ktop(i) = ktcon(i)
2023         kbot(i) = kbcon(i)
2024         kuo(i) = 1
2025       endif
2026     endif
2027   enddo
2028!
2029   do k = kts,kmax
2030     do i = its,ite
2031       if(cnvflg(i).and.rain(i).le.0.) then
2032          t1(i,k) = to(i,k)
2033          q1(i,k) = qo(i,k)
2034          u1(i,k) = uo(i,k)
2035          v1(i,k) = vo(i,k)
2036       endif
2037     enddo
2038   enddo
2039!
2040!  detrainment of cloud water and ice
2041!
2042   if (ncloud.gt.0) then
2043     do k = kts,kmax
2044       do i = its,ite
2045         if (cnvflg(i) .and. rain(i).gt.0.) then
2046           if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
2047             tem  = dellal(i,k) * xmb(i) * dt2
2048             tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2049             if (ncloud.ge.4) then
2050               qi2(i,k) = qi2(i,k) + tem * tem1            ! ice
2051               qc2(i,k) = qc2(i,k) + tem *(1.0-tem1)       ! water
2052             else
2053               qc2(i,k) = qc2(i,k) + tem
2054             endif
2055           endif
2056         endif
2057       enddo
2058     enddo
2059   endif
2060!
2061   end subroutine nsas2d
2062!===============================================================================
2063      REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2064!-------------------------------------------------------------------------------
2065      IMPLICIT NONE
2066!-------------------------------------------------------------------------------
2067      REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,      &
2068           xai,xbi,ttp,tr
2069      INTEGER :: ice
2070! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2071      ttp=t0c+0.01
2072      dldt=cvap-cliq
2073      xa=-dldt/rv
2074      xb=xa+hvap/(rv*ttp)
2075      dldti=cvap-cice
2076      xai=-dldti/rv
2077      xbi=xai+hsub/(rv*ttp)
2078      tr=ttp/t
2079      if(t.lt.ttp.and.ice.eq.1) then
2080        fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2081      else
2082        fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2083      endif
2084!
2085      if (t.lt.180.) then
2086        tr=ttp/180.
2087        if(t.lt.ttp.and.ice.eq.1) then
2088          fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2089        else
2090          fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2091        endif
2092      endif
2093!
2094      if (t.ge.330.) then
2095        tr=ttp/330
2096        if(t.lt.ttp.and.ice.eq.1) then
2097          fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2098        else
2099          fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2100        endif
2101      endif
2102!
2103! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2104      END FUNCTION fpvs
2105!===============================================================================
2106   subroutine nsasinit(rthcuten,rqvcuten,rqccuten,rqicuten,                    &
2107                      rucuten,rvcuten,                                         & 
2108                      restart,p_qc,p_qi,p_first_scalar,                        &
2109                      allowed_to_read,                                         &
2110                      ids, ide, jds, jde, kds, kde,                            &
2111                      ims, ime, jms, jme, kms, kme,                            &
2112                      its, ite, jts, jte, kts, kte                  )
2113!-------------------------------------------------------------------------------
2114   implicit none
2115!-------------------------------------------------------------------------------
2116   logical , intent(in)           ::  allowed_to_read,restart
2117   integer , intent(in)           ::  ids, ide, jds, jde, kds, kde,            &
2118                                      ims, ime, jms, jme, kms, kme,            &
2119                                      its, ite, jts, jte, kts, kte
2120   integer , intent(in)           ::  p_first_scalar, p_qi, p_qc
2121   real,     dimension( ims:ime , kms:kme , jms:jme ) , intent(out) ::         &
2122                                                              rthcuten,        &
2123                                                              rqvcuten,        &
2124                                                               rucuten,        &
2125                                                               rvcuten,        &
2126                                                              rqccuten,        &
2127                                                              rqicuten
2128   integer :: i, j, k, itf, jtf, ktf
2129   jtf=min0(jte,jde-1)
2130   ktf=min0(kte,kde-1)
2131   itf=min0(ite,ide-1)
2132   if(.not.restart)then
2133     do j=jts,jtf
2134     do k=kts,ktf
2135     do i=its,itf
2136       rthcuten(i,k,j)=0.
2137       rqvcuten(i,k,j)=0.
2138       rucuten(i,k,j)=0.   
2139       rvcuten(i,k,j)=0.   
2140     enddo
2141     enddo
2142     enddo
2143     if (p_qc .ge. p_first_scalar) then
2144        do j=jts,jtf
2145        do k=kts,ktf
2146        do i=its,itf
2147           rqccuten(i,k,j)=0.
2148        enddo
2149        enddo
2150        enddo
2151     endif
2152     if (p_qi .ge. p_first_scalar) then
2153        do j=jts,jtf
2154        do k=kts,ktf
2155        do i=its,itf
2156           rqicuten(i,k,j)=0.
2157        enddo
2158        enddo
2159        enddo
2160     endif
2161   endif
2162      end subroutine nsasinit
2163!
2164!==============================================================================
2165! NCEP SCV (Shallow Convection Scheme)
2166!==============================================================================
2167   subroutine nscv2d(delt,del,prsl,prsi,prslk,zl,zi,                           &
2168                 ncloud,qc2,qi2,q1,t1,rain,kbot,ktop,                          &
2169                 kuo,                                                          &
2170                 slimsk,dot,u1,v1,                                             &
2171                 cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2,                     &
2172                 cice,xls,psat,                                                &
2173                 hpbl,hfx,qfx,                                                 &
2174                 ids,ide, jds,jde, kds,kde,                                    &
2175                 ims,ime, jms,jme, kms,kme,                                    &
2176                 its,ite, jts,jte, kts,kte)
2177!
2178!-------------------------------------------------------------------------------
2179!
2180! subprogram:    nscv2d           computes shallow-convective heating and moisng
2181!
2182! abstract: computes non-precipitating convective heating and moistening
2183!   using a one cloud type arakawa-schubert convection scheme as in the ncep
2184!   sas scheme. the scheme has been operational at ncep gfs model since july 2010
2185!   the scheme includes updraft and downdraft effects, but the cloud depth is
2186!   limited less than 150 hpa.
2187!
2188! developed by jong-il han and hua-lu pan
2189!   implemented into wrf by jiheon jang and songyou hong
2190!   module with cpp-based options is available in grims
2191!
2192! program history log:
2193!   10-07-01 jong-il han  initial operational implementation at ncep gfs
2194!   10-12-01 jihyeon jang implemented into wrf
2195!
2196! subprograms called:
2197!   fpvs     - function to compute saturation vapor pressure
2198!
2199! references:
2200!   han and pan (2010, wea. forecasting)
2201!   
2202!-------------------------------------------------------------------------------
2203   implicit none
2204!-------------------------------------------------------------------------------
2205!  in/out variables
2206!
2207   integer         ::  ids,ide, jds,jde, kds,kde,                              &
2208                       ims,ime, jms,jme, kms,kme,                              &
2209                       its,ite, jts,jte, kts,kte
2210   real            ::  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
2211   real            ::  pi_,qmin_,t0c_
2212   real            ::  cice,xlv0,xls,psat
2213!
2214   real            ::  delt
2215   real            ::  del(its:ite,kts:kte),                                   &
2216                       prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme),           &
2217                       prsi(its:ite,kts:kte+1),zl(its:ite,kts:kte)
2218   integer         ::  ncloud
2219   real            ::  slimsk(ims:ime)
2220   real            ::  dot(its:ite,kts:kte)
2221   real            ::  hpbl(ims:ime)
2222   real            ::  rcs
2223   real            ::  hfx(ims:ime),qfx(ims:ime)
2224!
2225   real            ::  qi2(its:ite,kts:kte),qc2(its:ite,kts:kte)
2226   real            ::  q1(its:ite,kts:kte),                                    &
2227                       t1(its:ite,kts:kte),                                    &
2228                       u1(its:ite,kts:kte),                                    &
2229                       v1(its:ite,kts:kte)
2230   integer         ::  kuo(its:ite)
2231!
2232   real            ::  rain(its:ite)
2233   integer         ::  kbot(its:ite),ktop(its:ite)
2234!
2235!  local variables and arrays
2236!
2237   integer         ::  i,j,indx, jmn, k, kk, km1
2238   integer         ::  kpbl(its:ite)
2239!
2240   real            ::  dellat,                                                 &
2241                       desdt,   deta,    detad,   dg,                          &
2242                       dh,      dhh,     dlnsig,  dp,                          &
2243                       dq,      dqsdp,   dqsdt,   dt,                          &
2244                       dt2,     dtmax,   dtmin,                                &
2245                       dv1h,    dv2h,    dv3h,                                 &
2246                       dv1q,    dv2q,    dv3q,                                 &
2247                       dv1u,    dv2u,    dv3u,                                 &
2248                       dv1v,    dv2v,    dv3v,                                 &
2249                       dz,      dz1,     e1,      clam,                        &
2250                       aafac,                                                  &
2251                       es,      etah,                                          &
2252                       evef,    evfact,  evfactl,                              &
2253                       factor,  fjcap,                                         &
2254                       gamma,   pprime,  betaw,                                &
2255                       qlk,     qrch,    qs,                                   &
2256                       rfact,   shear,   tem1,                                 &
2257                       tem2,    val,     val1,                                 &
2258                       val2,    w1,      w1l,     w1s,                         &
2259                       w2,      w2l,     w2s,     w3,                          &
2260                       w3l,     w3s,     w4,      w4l,                         &
2261                       w4s,     tem,     ptem,    ptem1,                       &
2262                       pgcon
2263!
2264   integer         ::  kb(its:ite), kbcon(its:ite), kbcon1(its:ite),           &
2265                       ktcon(its:ite), ktcon1(its:ite),                        &
2266                       kbm(its:ite), kmax(its:ite)
2267!
2268   real            ::  aa1(its:ite),                                           &
2269                       delhbar(its:ite), delq(its:ite),                        &
2270                       delq2(its:ite),   delqev(its:ite), rntot(its:ite),      &
2271                       delqbar(its:ite), deltbar(its:ite),                     &
2272                       deltv(its:ite),   edt(its:ite),                         &
2273                       wstar(its:ite),   sflx(its:ite),                        &
2274                       pdot(its:ite),    po(its:ite,kts:kte),                  &
2275                       qcond(its:ite),   qevap(its:ite),  hmax(its:ite),       &
2276                       vshear(its:ite),                                        &
2277                       xlamud(its:ite),  xmb(its:ite),    xmbmax(its:ite)
2278   real            ::  delubar(its:ite), delvbar(its:ite)
2279!
2280   real            ::  cincr
2281!
2282   real            ::  thx(its:ite, kts:kte)
2283   real            ::  rhox(its:ite)
2284   real            ::  tvcon
2285!
2286   real            ::  p(its:ite,kts:kte),       to(its:ite,kts:kte),          &
2287                       qo(its:ite,kts:kte),      qeso(its:ite,kts:kte),        &
2288                       uo(its:ite,kts:kte),      vo(its:ite,kts:kte)
2289!
2290!  cloud water
2291!
2292   real            ::  qlko_ktcon(its:ite),     dellal(its:ite,kts:kte),       &
2293                       dbyo(its:ite,kts:kte),                                  &
2294                       xlamue(its:ite,kts:kte),                                &
2295                       heo(its:ite,kts:kte),    heso(its:ite,kts:kte),         &
2296                       dellah(its:ite,kts:kte), dellaq(its:ite,kts:kte),       &
2297                       dellau(its:ite,kts:kte), dellav(its:ite,kts:kte),       &
2298                       ucko(its:ite,kts:kte),   vcko(its:ite,kts:kte),         &
2299                       hcko(its:ite,kts:kte),   qcko(its:ite,kts:kte),         &
2300                       eta(its:ite,kts:kte),    zi(its:ite,kts:kte+1),         &
2301                       pwo(its:ite,kts:kte)
2302!
2303   logical         ::  totflg, cnvflg(its:ite), flg(its:ite)
2304!
2305!  physical parameters
2306!
2307   real,parameter  ::  c0=.002,c1=5.e-4
2308   real,parameter  ::  cincrmax=180.,cincrmin=120.,dthk=25.
2309   real            ::  el2orc,fact1,fact2,eps
2310   real,parameter  ::  h1=0.33333333
2311   real,parameter  ::  tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)
2312!
2313!-------------------------------------------------------------------------------
2314!
2315   pi_ = 3.14159
2316   qmin_ = 1.0e-30
2317   t0c_ = 273.15
2318   xlv0 = hvap_
2319      km1 = kte - 1
2320!
2321!  compute surface buoyancy flux
2322!
2323      do k = kts,kte
2324        do i = its,ite
2325          thx(i,k) = t1(i,k)/prslk(i,k)
2326        enddo
2327      enddo
2328!
2329      do i=its,ite
2330         tvcon = (1.+fv_*q1(i,1))
2331         rhox(i) = prsl(i,1)*1.e3/(rd_*t1(i,1)*tvcon)
2332      enddo
2333!
2334      do i=its,ite
2335!        sflx(i) = heat(i)+fv_*t1(i,1)*evap(i)
2336         sflx(i) = hfx(i)/rhox(i)/cp_ + qfx(i)/rhox(i)*fv_*thx(i,1)
2337      enddo
2338!
2339!  initialize arrays
2340!
2341      do i=its,ite
2342        cnvflg(i) = .true.
2343        if(kuo(i).eq.1) cnvflg(i) = .false.
2344        if(sflx(i).le.0.) cnvflg(i) = .false.
2345        if(cnvflg(i)) then
2346          kbot(i)=kte+1
2347          ktop(i)=0
2348        endif
2349        rain(i)=0.
2350        kbcon(i)=kte
2351        ktcon(i)=1
2352        kb(i)=kte
2353        pdot(i) = 0.
2354        qlko_ktcon(i) = 0.
2355        edt(i)  = 0.
2356        aa1(i)  = 0.
2357        vshear(i) = 0.
2358      enddo
2359!
2360      totflg = .true.
2361      do i=its,ite
2362        totflg = totflg .and. (.not. cnvflg(i))
2363      enddo
2364      if(totflg) return
2365!
2366      dt2   =  delt
2367      val   =         1200.
2368      dtmin = max(dt2, val )
2369      val   =         3600.
2370      dtmax = max(dt2, val )
2371!  model tunable parameters are all here
2372      clam    = .3
2373      aafac   = .1
2374      betaw   = .03
2375      evfact  = 0.3
2376      evfactl = 0.3
2377      pgcon   = 0.55    ! Zhang & Wu (2003,JAS)
2378      val     =           1.
2379!
2380! define miscellaneous values
2381!
2382     el2orc = hvap_*hvap_/(rv_*cp_)
2383     eps    = rd_/rv_
2384     fact1  = (cvap_-cliq_)/rv_
2385     fact2  = hvap_/rv_-fact1*t0c_
2386!
2387      w1l     = -8.e-3
2388      w2l     = -4.e-2
2389      w3l     = -5.e-3
2390      w4l     = -5.e-4
2391      w1s     = -2.e-4
2392      w2s     = -2.e-3
2393      w3s     = -1.e-3
2394      w4s     = -2.e-5
2395!
2396!  define top layer for search of the downdraft originating layer
2397!  and the maximum thetae for updraft
2398!
2399      do i=its,ite
2400        kbm(i)   = kte
2401        kmax(i)  = kte
2402      enddo
2403!     
2404      do k = kts, kte
2405        do i=its,ite
2406          if (prsl(i,k).gt.prsi(i,1)*0.70) kbm(i) = k + 1
2407          if (prsl(i,k).gt.prsi(i,1)*0.60) kmax(i) = k + 1
2408        enddo
2409      enddo
2410      do i=its,ite
2411        kbm(i)   = min(kbm(i),kmax(i))
2412      enddo
2413!
2414!  hydrostatic height assume zero terr and compute
2415!  updraft entrainment rate as an inverse function of height
2416!
2417      do k = kts, km1
2418        do i=its,ite
2419          xlamue(i,k) = clam / zi(i,k)
2420        enddo
2421      enddo
2422      do i=its,ite
2423        xlamue(i,kte) = xlamue(i,km1)
2424      enddo
2425!
2426!  pbl height
2427!
2428      do i=its,ite
2429        flg(i) = cnvflg(i)
2430        kpbl(i)= 1
2431      enddo
2432!
2433      do k = kts+1, km1
2434        do i=its,ite
2435          if (flg(i).and.zl(i,k).le.hpbl(i)) then
2436            kpbl(i) = k
2437          else
2438            flg(i) = .false.
2439          endif
2440        enddo
2441      enddo
2442!
2443      do i=its,ite
2444        kpbl(i)= min(kpbl(i),kbm(i))
2445      enddo
2446!
2447!   convert surface pressure to mb from cb
2448!
2449      rcs = 1.
2450      do k = kts, kte
2451        do i =its,ite
2452          if (cnvflg(i) .and. k .le. kmax(i)) then
2453            p(i,k) = prsl(i,k) * 10.0
2454            eta(i,k)  = 1.
2455            hcko(i,k) = 0.
2456            qcko(i,k) = 0.
2457            ucko(i,k) = 0.
2458            vcko(i,k) = 0.
2459            dbyo(i,k) = 0.
2460            pwo(i,k)  = 0.
2461            dellal(i,k) = 0.
2462            to(i,k)   = t1(i,k)
2463            qo(i,k)   = q1(i,k)
2464            uo(i,k)   = u1(i,k) * rcs
2465            vo(i,k)   = v1(i,k) * rcs
2466          endif
2467        enddo
2468      enddo
2469!
2470!
2471!  column variables
2472!  p is pressure of the layer (mb)
2473!  t is temperature at t-dt (k)..tn
2474!  q is mixing ratio at t-dt (kg/kg)..qn
2475!  to is temperature at t+dt (k)... this is after advection and turbulan
2476!  qo is mixing ratio at t+dt (kg/kg)..q1
2477!
2478      do k = kts, kte
2479        do i=its,ite
2480          if (cnvflg(i) .and. k .le. kmax(i)) then
2481            qeso(i,k) = 0.01 * fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2482            qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
2483            val1      =             1.e-8
2484            qeso(i,k) = max(qeso(i,k), val1)
2485            val2      =           1.e-10
2486            qo(i,k)   = max(qo(i,k), val2 )
2487          endif
2488        enddo
2489      enddo
2490!
2491!  compute moist static energy
2492!
2493      do k = kts,kte
2494        do i=its,ite
2495          if (cnvflg(i) .and. k .le. kmax(i)) then
2496            tem       = g_ * zl(i,k) + cp_ * to(i,k)
2497            heo(i,k)  = tem  + hvap_ * qo(i,k)
2498            heso(i,k) = tem  + hvap_ * qeso(i,k)
2499          endif
2500        enddo
2501      enddo
2502!
2503!  determine level with largest moist static energy within pbl
2504!  this is the level where updraft starts
2505!
2506      do i=its,ite
2507         if (cnvflg(i)) then
2508            hmax(i) = heo(i,1)
2509            kb(i) = 1
2510         endif
2511      enddo
2512!
2513      do k = kts+1, kte
2514        do i=its,ite
2515          if (cnvflg(i).and.k.le.kpbl(i)) then
2516            if(heo(i,k).gt.hmax(i)) then
2517              kb(i)   = k
2518              hmax(i) = heo(i,k)
2519            endif
2520          endif
2521        enddo
2522      enddo
2523!
2524      do k = kts, km1
2525        do i=its,ite
2526          if (cnvflg(i) .and. k .le. kmax(i)-1) then
2527            dz      = .5 * (zl(i,k+1) - zl(i,k))
2528            dp      = .5 * (p(i,k+1) - p(i,k))
2529            es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2530            pprime  = p(i,k+1) + (eps-1.) * es
2531            qs      = eps * es / pprime
2532            dqsdp   = - qs / pprime
2533            desdt   = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2534            dqsdt   = qs * p(i,k+1) * desdt / (es * pprime)
2535            gamma   = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2536            dt      = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
2537            dq      = dqsdt * dt + dqsdp * dp
2538            to(i,k) = to(i,k+1) + dt
2539            qo(i,k) = qo(i,k+1) + dq
2540            po(i,k) = .5 * (p(i,k) + p(i,k+1))
2541          endif
2542        enddo
2543      enddo
2544!
2545      do k = kts, km1
2546        do i=its,ite
2547          if (cnvflg(i) .and. k .le. kmax(i)-1) then
2548            qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2549            qeso(i,k) = eps * qeso(i,k) / (po(i,k) + (eps-1.) * qeso(i,k))
2550            val1      =             1.e-8
2551            qeso(i,k) = max(qeso(i,k), val1)
2552            val2      =           1.e-10
2553            qo(i,k)   = max(qo(i,k), val2 )
2554            heo(i,k)  = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                      &
2555                        cp_ * to(i,k) + hvap_ * qo(i,k)
2556            heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                      &
2557                        cp_ * to(i,k) + hvap_ * qeso(i,k)
2558            uo(i,k)   = .5 * (uo(i,k) + uo(i,k+1))
2559            vo(i,k)   = .5 * (vo(i,k) + vo(i,k+1))
2560          endif
2561        enddo
2562      enddo
2563!
2564!  look for the level of free convection as cloud base
2565!
2566      do i=its,ite
2567        flg(i)   = cnvflg(i)
2568        if(flg(i)) kbcon(i) = kmax(i)
2569      enddo
2570!
2571      do k = kts+1, km1
2572        do i=its,ite
2573          if (flg(i).and.k.lt.kbm(i)) then
2574            if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
2575              kbcon(i) = k
2576              flg(i)   = .false.
2577            endif
2578          endif
2579        enddo
2580      enddo
2581!
2582      do i=its,ite
2583        if(cnvflg(i)) then
2584          if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
2585        endif
2586      enddo
2587!
2588      totflg = .true.
2589      do i=its,ite
2590        totflg = totflg .and. (.not. cnvflg(i))
2591      enddo
2592      if(totflg) return
2593!
2594!  determine critical convective inhibition
2595!  as a function of vertical velocity at cloud base.
2596!
2597      do i=its,ite
2598        if(cnvflg(i)) then
2599          pdot(i)  = 10.* dot(i,kbcon(i))
2600        endif
2601      enddo
2602!
2603      do i=its,ite
2604        if(cnvflg(i)) then
2605          if(slimsk(i).eq.1.) then
2606            w1 = w1l
2607            w2 = w2l
2608            w3 = w3l
2609            w4 = w4l
2610          else
2611            w1 = w1s
2612            w2 = w2s
2613            w3 = w3s
2614            w4 = w4s
2615          endif
2616          if(pdot(i).le.w4) then
2617            ptem = (pdot(i) - w4) / (w3 - w4)
2618          elseif(pdot(i).ge.-w4) then
2619            ptem = - (pdot(i) + w4) / (w4 - w3)
2620          else
2621            ptem = 0.
2622          endif
2623          val1    =             -1.
2624          ptem = max(ptem,val1)
2625          val2    =             1.
2626          ptem = min(ptem,val2)
2627          ptem = 1. - ptem
2628          ptem1= .5*(cincrmax-cincrmin)
2629          cincr = cincrmax - ptem * ptem1
2630          tem1 = p(i,kb(i)) - p(i,kbcon(i))
2631          if(tem1.gt.cincr) then
2632             cnvflg(i) = .false.
2633          endif
2634        endif
2635      enddo
2636!
2637      totflg = .true.
2638      do i=its,ite
2639        totflg = totflg .and. (.not. cnvflg(i))
2640      enddo
2641      if(totflg) return
2642!
2643!  assume the detrainment rate for the updrafts to be same as
2644!  the entrainment rate at cloud base
2645!
2646      do i = its,ite
2647        if(cnvflg(i)) then
2648          xlamud(i) = xlamue(i,kbcon(i))
2649        endif
2650      enddo
2651!
2652!  determine updraft mass flux for the subcloud layers
2653!
2654       do k = km1, kts, -1
2655        do i = its,ite
2656          if (cnvflg(i)) then
2657            if(k.lt.kbcon(i).and.k.ge.kb(i)) then
2658              dz       = zi(i,k+1) - zi(i,k)
2659              ptem     = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
2660              eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
2661            endif
2662          endif
2663        enddo
2664      enddo
2665!
2666!  compute mass flux above cloud base
2667!
2668      do k = kts+1, km1
2669        do i = its,ite
2670         if(cnvflg(i))then
2671           if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
2672              dz       = zi(i,k) - zi(i,k-1)
2673              ptem     = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
2674              eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
2675           endif
2676         endif
2677        enddo
2678      enddo
2679!
2680!  compute updraft cloud property
2681!
2682      do i = its,ite
2683        if(cnvflg(i)) then
2684          indx         = kb(i)
2685          hcko(i,indx) = heo(i,indx)
2686          ucko(i,indx) = uo(i,indx)
2687          vcko(i,indx) = vo(i,indx)
2688        endif
2689      enddo
2690!
2691      do k = kts+1, km1
2692        do i = its,ite
2693          if (cnvflg(i)) then
2694            if(k.gt.kb(i).and.k.lt.kmax(i)) then
2695              dz   = zi(i,k) - zi(i,k-1)
2696              tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2697              tem1 = 0.5 * xlamud(i) * dz
2698              factor = 1. + tem - tem1
2699              ptem = 0.5 * tem + pgcon
2700              ptem1= 0.5 * tem - pgcon
2701              hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                      &
2702                          (heo(i,k)+heo(i,k-1)))/factor
2703              ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)                  &
2704                          +ptem1*uo(i,k-1))/factor
2705              vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)                  &
2706                          +ptem1*vo(i,k-1))/factor
2707              dbyo(i,k) = hcko(i,k) - heso(i,k)
2708            endif
2709          endif
2710        enddo
2711      enddo
2712!
2713!   taking account into convection inhibition due to existence of
2714!    dry layers below cloud base
2715!
2716      do i=its,ite
2717        flg(i) = cnvflg(i)
2718        kbcon1(i) = kmax(i)
2719      enddo
2720!
2721      do k = kts+1, km1
2722        do i=its,ite
2723          if (flg(i).and.k.lt.kbm(i)) then
2724            if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
2725              kbcon1(i) = k
2726              flg(i)    = .false.
2727            endif
2728          endif
2729        enddo
2730      enddo
2731!
2732      do i=its,ite
2733        if(cnvflg(i)) then
2734          if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
2735        endif
2736      enddo
2737!
2738      do i=its,ite
2739        if(cnvflg(i)) then
2740          tem = p(i,kbcon(i)) - p(i,kbcon1(i))
2741          if(tem.gt.dthk) then
2742             cnvflg(i) = .false.
2743          endif
2744        endif
2745      enddo
2746!
2747      totflg = .true.
2748      do i = its,ite
2749        totflg = totflg .and. (.not. cnvflg(i))
2750      enddo
2751      if(totflg) return
2752!
2753!  determine first guess cloud top as the level of zero buoyancy
2754!    limited to the level of sigma=0.7
2755!
2756      do i = its,ite
2757        flg(i) = cnvflg(i)
2758        if(flg(i)) ktcon(i) = kbm(i)
2759      enddo
2760!
2761      do k = kts+1, km1
2762        do i=its,ite
2763          if (flg(i).and.k .lt. kbm(i)) then
2764            if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
2765               ktcon(i) = k
2766               flg(i)   = .false.
2767            endif
2768          endif
2769        enddo
2770      enddo
2771!
2772!  specify upper limit of mass flux at cloud base
2773!
2774      do i = its,ite
2775        if(cnvflg(i)) then
2776          k = kbcon(i)
2777          dp = 1000. * del(i,k)
2778          xmbmax(i) = dp / (g_ * dt2)
2779        endif
2780      enddo
2781!
2782!  compute cloud moisture property and precipitation
2783!
2784      do i = its,ite
2785        if (cnvflg(i)) then
2786          aa1(i) = 0.
2787          qcko(i,kb(i)) = qo(i,kb(i))
2788        endif
2789      enddo
2790!
2791      do k = kts+1, km1
2792        do i = its,ite
2793          if (cnvflg(i)) then
2794            if(k.gt.kb(i).and.k.lt.ktcon(i)) then
2795              dz    = zi(i,k) - zi(i,k-1)
2796              gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2797              qrch = qeso(i,k)                                                 &
2798                   + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2799              tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2800              tem1 = 0.5 * xlamud(i) * dz
2801              factor = 1. + tem - tem1
2802              qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                      &
2803                          (qo(i,k)+qo(i,k-1)))/factor
2804              dq = eta(i,k) * (qcko(i,k) - qrch)
2805!
2806!             rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
2807!
2808!  below lfc check if there is excess moisture to release latent heat
2809!
2810              if(k.ge.kbcon(i).and.dq.gt.0.) then
2811                etah = .5 * (eta(i,k) + eta(i,k-1))
2812                if(ncloud.gt.0) then
2813                  dp = 1000. * del(i,k)
2814                  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2815                  dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2816                else
2817                  qlk = dq / (eta(i,k) + etah * c0 * dz)
2818                endif
2819                aa1(i) = aa1(i) - dz * g_ * qlk
2820                qcko(i,k)= qlk + qrch
2821                pwo(i,k) = etah * c0 * dz * qlk
2822              endif
2823            endif
2824          endif
2825        enddo
2826      enddo
2827!
2828!  calculate cloud work function
2829!
2830      do k = kts+1, km1
2831        do i = its,ite
2832          if (cnvflg(i)) then
2833            if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
2834              dz1 = zl(i,k+1) - zl(i,k)       
2835              gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2836              rfact =  1. + fv_ * cp_ * gamma                                  &
2837                      * to(i,k) / hvap_
2838              aa1(i) = aa1(i) + dz1 * (g_ / (cp_ * to(i,k)))                   &
2839                      * dbyo(i,k) / (1. + gamma) * rfact
2840              val = 0.
2841              aa1(i)=aa1(i)+ dz1 * g_ * fv_ *                                  &
2842                      max(val,(qeso(i,k) - qo(i,k)))
2843            endif
2844          endif
2845        enddo
2846      enddo
2847!
2848      do i = its,ite
2849        if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
2850      enddo
2851!
2852      totflg = .true.
2853      do i=its,ite
2854        totflg = totflg .and. (.not. cnvflg(i))
2855      enddo
2856      if(totflg) return
2857!
2858!  estimate the convective overshooting as the level
2859!    where the [aafac * cloud work function] becomes zero,
2860!    which is the final cloud top limited to the level of sigma=0.7
2861!
2862      do i = its,ite
2863        if (cnvflg(i)) then
2864          aa1(i) = aafac * aa1(i)
2865        endif
2866      enddo
2867!
2868      do i = its,ite
2869        flg(i) = cnvflg(i)
2870        ktcon1(i) = kbm(i)
2871      enddo
2872!
2873      do k = kts+1, km1
2874        do i = its,ite
2875          if (flg(i)) then
2876            if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
2877              dz1 = zl(i,k+1) - zl(i,k)
2878              gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2879              rfact =  1. + fv_ * cp_ * gamma                                  &
2880                      * to(i,k) / hvap_
2881              aa1(i) = aa1(i) +                                                &
2882                      dz1 * (g_ / (cp_ * to(i,k)))                             &
2883                      * dbyo(i,k) / (1. + gamma) * rfact
2884              if(aa1(i).lt.0.) then
2885                ktcon1(i) = k
2886                flg(i) = .false.
2887              endif
2888            endif
2889          endif
2890        enddo
2891      enddo
2892!
2893!  compute cloud moisture property, detraining cloud water
2894!    and precipitation in overshooting layers
2895!
2896      do k = kts+1, km1
2897        do i = its,ite
2898          if (cnvflg(i)) then
2899            if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
2900              dz    = zi(i,k) - zi(i,k-1)
2901              gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2902              qrch = qeso(i,k)                                                 &
2903                   + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2904              tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2905              tem1 = 0.5 * xlamud(i) * dz
2906              factor = 1. + tem - tem1
2907              qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                      &
2908                          (qo(i,k)+qo(i,k-1)))/factor
2909              dq = eta(i,k) * (qcko(i,k) - qrch)
2910!
2911!  check if there is excess moisture to release latent heat
2912!
2913              if(dq.gt.0.) then
2914                etah = .5 * (eta(i,k) + eta(i,k-1))
2915                if(ncloud.gt.0) then
2916                  dp = 1000. * del(i,k)
2917                  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2918                  dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2919                else
2920                  qlk = dq / (eta(i,k) + etah * c0 * dz)
2921                endif
2922                qcko(i,k) = qlk + qrch
2923                pwo(i,k) = etah * c0 * dz * qlk
2924              endif
2925            endif
2926          endif
2927        enddo
2928      enddo
2929!
2930! exchange ktcon with ktcon1
2931!
2932      do i = its,ite
2933        if(cnvflg(i)) then
2934          kk = ktcon(i)
2935          ktcon(i) = ktcon1(i)
2936          ktcon1(i) = kk
2937        endif
2938      enddo
2939!
2940!  this section is ready for cloud water
2941!
2942      if(ncloud.gt.0) then
2943!
2944!  compute liquid and vapor separation at cloud top
2945!
2946      do i = its,ite
2947        if(cnvflg(i)) then
2948          k = ktcon(i) - 1
2949          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2950          qrch = qeso(i,k)                                                     &
2951               + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2952          dq = qcko(i,k) - qrch
2953!
2954!  check if there is excess moisture to release latent heat
2955!
2956          if(dq.gt.0.) then
2957            qlko_ktcon(i) = dq
2958            qcko(i,k) = qrch
2959          endif
2960        endif
2961      enddo
2962!
2963      endif
2964!
2965!--- compute precipitation efficiency in terms of windshear
2966!
2967      do i = its,ite
2968        if(cnvflg(i)) then
2969          vshear(i) = 0.
2970        endif
2971      enddo
2972!
2973      do k = kts+1,kte
2974        do i = its,ite
2975          if (cnvflg(i)) then
2976            if(k.gt.kb(i).and.k.le.ktcon(i)) then
2977              shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2                             &
2978                        + (vo(i,k)-vo(i,k-1)) ** 2)
2979              vshear(i) = vshear(i) + shear
2980            endif
2981          endif
2982        enddo
2983      enddo
2984!
2985      do i = its,ite
2986        if(cnvflg(i)) then
2987          vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
2988          e1=1.591-.639*vshear(i)                                              &
2989             +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
2990          edt(i)=1.-e1
2991          val =         .9
2992          edt(i) = min(edt(i),val)
2993          val =         .0
2994          edt(i) = max(edt(i),val)
2995        endif
2996      enddo
2997!
2998!--- what would the change be, that a cloud with unit mass
2999!--- will do to the environment?
3000!
3001      do k = kts,kte
3002        do i = its,ite
3003          if(cnvflg(i) .and. k .le. kmax(i)) then
3004            dellah(i,k) = 0.
3005            dellaq(i,k) = 0.
3006            dellau(i,k) = 0.
3007            dellav(i,k) = 0.
3008          endif
3009        enddo
3010      enddo
3011!
3012!--- changed due to subsidence and entrainment
3013!
3014      do k = kts+1, km1
3015        do i = its,ite
3016          if (cnvflg(i)) then
3017            if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3018              dp = 1000. * del(i,k)
3019              dz = zi(i,k) - zi(i,k-1)
3020!
3021              dv1h = heo(i,k)
3022              dv2h = .5 * (heo(i,k) + heo(i,k-1))
3023              dv3h = heo(i,k-1)
3024              dv1q = qo(i,k)
3025              dv2q = .5 * (qo(i,k) + qo(i,k-1))
3026              dv3q = qo(i,k-1)
3027              dv1u = uo(i,k)
3028              dv2u = .5 * (uo(i,k) + uo(i,k-1))
3029              dv3u = uo(i,k-1)
3030              dv1v = vo(i,k)
3031              dv2v = .5 * (vo(i,k) + vo(i,k-1))
3032              dv3v = vo(i,k-1)
3033!
3034              tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
3035              tem1 = xlamud(i)
3036!
3037              dellah(i,k) = dellah(i,k) +                                      &
3038          ( eta(i,k)*dv1h - eta(i,k-1)*dv3h                                    &
3039         -  tem*eta(i,k-1)*dv2h*dz                                             &
3040         +  tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz                      &
3041              ) *g_/dp
3042!
3043              dellaq(i,k) = dellaq(i,k) +                                      &
3044          ( eta(i,k)*dv1q - eta(i,k-1)*dv3q                                    &
3045         -  tem*eta(i,k-1)*dv2q*dz                                             &
3046         +  tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz                      &
3047              ) *g_/dp
3048!
3049              dellau(i,k) = dellau(i,k) +                                      &
3050          ( eta(i,k)*dv1u - eta(i,k-1)*dv3u                                    &
3051         -  tem*eta(i,k-1)*dv2u*dz                                             &
3052         +  tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz                      &
3053         -  pgcon*eta(i,k-1)*(dv1u-dv3u)                                       &
3054              ) *g_/dp
3055!
3056              dellav(i,k) = dellav(i,k) +                                      &
3057          ( eta(i,k)*dv1v - eta(i,k-1)*dv3v                                    &
3058         -  tem*eta(i,k-1)*dv2v*dz                                             &
3059         +  tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz                      &
3060         -  pgcon*eta(i,k-1)*(dv1v-dv3v)                                       &
3061              ) *g_/dp
3062!
3063            endif
3064          endif
3065        enddo
3066      enddo
3067!
3068!------- cloud top
3069!
3070      do i = its,ite
3071        if(cnvflg(i)) then
3072          indx = ktcon(i)
3073          dp = 1000. * del(i,indx)
3074          dv1h = heo(i,indx-1)
3075          dellah(i,indx) = eta(i,indx-1) *                                     &
3076                          (hcko(i,indx-1) - dv1h) * g_ / dp
3077          dv1q = qo(i,indx-1)
3078          dellaq(i,indx) = eta(i,indx-1) *                                     &
3079                          (qcko(i,indx-1) - dv1q) * g_ / dp
3080          dv1u = uo(i,indx-1)
3081          dellau(i,indx) = eta(i,indx-1) *                                     &
3082                          (ucko(i,indx-1) - dv1u) * g_ / dp
3083          dv1v = vo(i,indx-1)
3084          dellav(i,indx) = eta(i,indx-1) *                                     &
3085                          (vcko(i,indx-1) - dv1v) * g_ / dp
3086!
3087!  cloud water
3088!
3089          dellal(i,indx) = eta(i,indx-1) *                                     &
3090                          qlko_ktcon(i) * g_ / dp
3091        endif
3092      enddo
3093!
3094!  mass flux at cloud base for shallow convection
3095!  (Grant, 2001)
3096!
3097      do i= its,ite
3098        if(cnvflg(i)) then
3099          k = kbcon(i)
3100          ptem = g_*sflx(i)*hpbl(i)/t1(i,1)
3101          wstar(i) = ptem**h1
3102          tem = po(i,k)*100. / (rd_*t1(i,k))
3103          xmb(i) = betaw*tem*wstar(i)
3104          xmb(i) = min(xmb(i),xmbmax(i))
3105        endif
3106      enddo
3107!
3108      do k = kts,kte
3109        do i = its,ite
3110          if (cnvflg(i) .and. k .le. kmax(i)) then
3111            qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
3112            qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
3113            val     =             1.e-8
3114            qeso(i,k) = max(qeso(i,k), val )
3115          endif
3116        enddo
3117      enddo
3118!
3119      do i = its,ite
3120        delhbar(i) = 0.
3121        delqbar(i) = 0.
3122        deltbar(i) = 0.
3123        delubar(i) = 0.
3124        delvbar(i) = 0.
3125        qcond(i) = 0.
3126      enddo
3127!
3128      do k = kts,kte
3129        do i = its,ite
3130          if (cnvflg(i)) then
3131            if(k.gt.kb(i).and.k.le.ktcon(i)) then
3132              dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
3133              t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
3134              q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
3135              tem = 1./rcs
3136              u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
3137              v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
3138              dp = 1000. * del(i,k)
3139              delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
3140              delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
3141              deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
3142              delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
3143              delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
3144            endif
3145          endif
3146        enddo
3147      enddo
3148!
3149      do k = kts,kte
3150        do i = its,ite
3151          if (cnvflg(i)) then
3152            if(k.gt.kb(i).and.k.le.ktcon(i)) then
3153              qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls &
3154                        ,psat,t0c_)
3155              qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
3156              val     =             1.e-8
3157              qeso(i,k) = max(qeso(i,k), val )
3158            endif
3159          endif
3160        enddo
3161      enddo
3162!
3163      do i = its,ite
3164        rntot(i) = 0.
3165        delqev(i) = 0.
3166        delq2(i) = 0.
3167        flg(i) = cnvflg(i)
3168      enddo
3169!
3170      do k = kte, kts, -1
3171        do i = its,ite
3172          if (cnvflg(i)) then
3173            if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3174              rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
3175            endif
3176          endif
3177        enddo
3178      enddo
3179!
3180! evaporating rain
3181!
3182      do k = kte, kts, -1
3183        do i = its,ite
3184          if (k .le. kmax(i)) then
3185            deltv(i) = 0.
3186            delq(i) = 0.
3187            qevap(i) = 0.
3188            if(cnvflg(i)) then
3189              if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3190                rain(i) = rain(i) + pwo(i,k) * xmb(i) * .001 * dt2
3191              endif
3192            endif
3193            if(flg(i).and.k.lt.ktcon(i)) then
3194              evef = edt(i) * evfact
3195              if(slimsk(i).eq.1.) evef=edt(i) * evfactl
3196              qcond(i) = evef * (q1(i,k) - qeso(i,k))                          &
3197                       / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
3198              dp = 1000. * del(i,k)
3199              if(rain(i).gt.0..and.qcond(i).lt.0.) then
3200                qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rain(i))))
3201                qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
3202                delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
3203              endif
3204              if(rain(i).gt.0..and.qcond(i).lt.0..and.                         &
3205                 delq2(i).gt.rntot(i)) then
3206                qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
3207                flg(i) = .false.
3208              endif
3209              if(rain(i).gt.0..and.qevap(i).gt.0.) then
3210                tem  = .001 * dp / g_
3211                tem1 = qevap(i) * tem
3212                if(tem1.gt.rain(i)) then
3213                  qevap(i) = rain(i) / tem
3214                  rain(i) = 0.
3215                else
3216                  rain(i) = rain(i) - tem1
3217                endif
3218                q1(i,k) = q1(i,k) + qevap(i)
3219                t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
3220                deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
3221                delq(i) =  + qevap(i)/dt2
3222                delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
3223              endif
3224              dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
3225              delqbar(i) = delqbar(i) + delq(i)*dp/g_
3226              deltbar(i) = deltbar(i) + deltv(i)*dp/g_
3227            endif
3228          endif
3229        enddo
3230      enddo
3231!
3232      do i = its,ite
3233        if(cnvflg(i)) then
3234          if(rain(i).lt.0..or..not.flg(i)) rain(i) = 0.
3235          ktop(i) = ktcon(i)
3236          kbot(i) = kbcon(i)
3237          kuo(i) = 0
3238        endif
3239      enddo
3240!
3241! cloud water
3242!
3243      if (ncloud.gt.0) then
3244!
3245      do k = kts, km1
3246        do i = its,ite
3247          if (cnvflg(i)) then
3248            if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
3249              tem  = dellal(i,k) * xmb(i) * dt2
3250              tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
3251              if (ncloud.ge.4) then
3252                qi2(i,k) = qi2(i,k) + tem * tem1            ! ice
3253                qc2(i,k) = qc2(i,k) + tem *(1.0-tem1)       ! water
3254              else
3255                qc2(i,k) = qc2(i,k) + tem
3256              endif
3257            endif
3258          endif
3259        enddo
3260      enddo
3261!
3262      endif
3263!
3264      end subroutine nscv2d
3265!-------------------------------------------------------------------------------
3266!
3267END MODULE module_cu_nsas
Note: See TracBrowser for help on using the repository browser.