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