updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / phys / module_cu_ksas.F
blobf6c92c507d2a42d0d2c183c27a17a6b4d72c0be6
5 MODULE module_cu_ksas
6 CONTAINS
7 !-------------------------------------------------------------------------------
8    subroutine cu_ksas(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                      hpbl_hold,znu,                                            &
15                      mp_physics,dx_factor_nsas,                                &
16                      p_qc,p_qi,p_first_scalar,                                 &
17                      pgcon,                                                    &
18                      cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,                      &
19                      cice,xls,psat,f_qi,f_qc,                                  &
20                      ids,ide, jds,jde, kds,kde,                                &
21                      ims,ime, jms,jme, kms,kme,                                &
22                      its,ite, jts,jte, kts,kte)
23 !-------------------------------------------------------------------------------
24    implicit none
25 !-------------------------------------------------------------------------------
26 !-- dt          time step (s)
27 !-- dx          grid interval (m)
28 !-- p3di        3d pressure (pa) at interface level
29 !-- p3d         3d pressure (pa)
30 !-- pi3d        3d exner function (dimensionless)
31 !-- z           height above sea level (m)
32 !-- qc3d        cloud water mixing ratio (kg/kg)
33 !-- qi3d        cloud ice mixing ratio (kg/kg)
34 !-- qv3d        3d water vapor mixing ratio (kg/kg)
35 !-- t3d         temperature (k)
36 !-- raincv      cumulus scheme precipitation (mm)
37 !-- w           vertical velocity (m/s)
38 !-- dz8w        dz between full levels (m)
39 !-- u3d         3d u-velocity interpolated to theta points (m/s)
40 !-- v3d         3d v-velocity interpolated to theta points (m/s)
41 !-- ids         start index for i in domain 
42 !-- ide         end index for i in domain
43 !-- jds         start index for j in domain
44 !-- jde         end index for j in domain
45 !-- kds         start index for k in domain
46 !-- kde         end index for k in domain
47 !-- ims         start index for i in memory
48 !-- ime         end index for i in memory
49 !-- jms         start index for j in memory
50 !-- jme         end index for j in memory
51 !-- kms         start index for k in memory
52 !-- kme         end index for k in memory 
53 !-- its         start index for i in tile
54 !-- ite         end index for i in tile
55 !-- jts         start index for j in tile
56 !-- jte         end index for j in tile
57 !-- kts         start index for k in tile
58 !-- kte         end index for k in tile
59 !-------------------------------------------------------------------------------
60    integer,  intent(in   )   ::       ids,ide, jds,jde, kds,kde,               &
61                                       ims,ime, jms,jme, kms,kme,               &
62                                       its,ite, jts,jte, kts,kte,               &
63                                       itimestep, stepcu,                       &
64                                       p_qc,p_qi,p_first_scalar
65    real,     intent(in   )   ::      cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,      &
66                                      cice,xls,psat
67    real,     intent(in   )   ::      dt,dx
68    real,     optional, intent(in ) :: pgcon
69    real,     dimension( ims:ime, kms:kme, jms:jme ),optional                  ,&
70              intent(inout)   ::                                       rthcuten,&
71                                                                        rucuten,&
72                                                                        rvcuten,&
73                                                                       rqccuten,&
74                                                                       rqicuten,&
75                                                                       rqvcuten
76    logical, optional ::                                              F_QC,F_QI
77    real,     dimension( ims:ime, kms:kme, jms:jme )                           ,&
78              intent(in   )   ::                                           qv3d,&
79                                                                           qc3d,&
80                                                                           qi3d,&
81                                                                          rho3d,&
82                                                                            p3d,&
83                                                                           pi3d,&
84                                                                            t3d
85    real,     dimension( ims:ime, kms:kme, jms:jme )                           ,&
86              intent(in   )   ::                                           p3di
87    real,     dimension( ims:ime, kms:kme, jms:jme )                           ,&
88              intent(in   )   ::                                           dz8w,&  
89                                                                              w
90    real,     dimension( ims:ime, jms:jme )                                    ,&
91              intent(inout) ::                                           raincv,&
92                                                                         pratec
93    real,     dimension( ims:ime, jms:jme )                                    ,&
94              intent(out) ::                                               hbot,&
95                                                                           htop
97    real,     dimension( ims:ime, jms:jme )                                    ,&
98              intent(in   ) ::                                            xland
100    real,     dimension( ims:ime, kms:kme, jms:jme )                           ,&
101               intent(in   )   ::                                           u3d,&
102                                                                            v3d
103    logical,  dimension( ims:ime, jms:jme )                                    ,&
104              intent(inout) ::                                      cu_act_flag
106    real,     dimension( ims:ime, jms:jme )                                    ,&
107               intent(in   )   ::                                          hpbl,&
108                                                                      hpbl_hold,&
109                                                                            hfx,&
110                                                                            qfx
111    real,     dimension( kms:kme )                                             ,&
112               intent(in   )   ::                                           znu 
113    integer,   intent(in   )   ::                                    mp_physics
114    integer,   intent(in   )   ::                                dx_factor_nsas 
115    integer :: ncloud
117 !local
119    real,  dimension( its:ite, jts:jte )  ::                            raincv1,&
120                                                                        raincv2,&
121                                                                        pratec1,&
122                                                                        pratec2
123    real,   dimension( its:ite, kts:kte )  ::                               del,&
124                                                                          prsll,&
125                                                                            dot,&
126                                                                             u1,&
127                                                                             v1,&
128                                                                             t1,&
129                                                                            q1, &
130                                                                            qc2,&
131                                                                            qi2
132    real,   dimension( its:ite, kts:kte+1 )  ::                           prsii,&
133                                                                            zii
134    real,   dimension( its:ite, kts:kte )  ::                               zll 
135    real,   dimension( its:ite)  ::                                         rain
136    real ::                                                          delt,rdelt
137    integer, dimension (its:ite)  ::                                       kbot,&
138                                                                           ktop,&
139                                                                           icps
140    real :: pgcon_use
141    integer ::  i,j,k,kp, kbmax,kbm,kmax
143 ! microphysics scheme --> ncloud 
145    if (mp_physics .eq. 0) then
146      ncloud = 0
147    elseif ( mp_physics .eq. 1 .or. mp_physics .eq. 3 ) then
148      ncloud = 1
149    else
150      ncloud = 2
151    endif
153    if(present(pgcon)) then
154      pgcon_use = pgcon
155    else
156 !    pgcon_use  = 0.7     ! Gregory et al. (1997, QJRMS)
157      pgcon_use  = 0.55    ! Zhang & Wu (2003,JAS)
158      ! 0.55 is a physically-based value used by GFS
159      ! HWRF uses 0.2, for model tuning purposes
160    endif
162    do j = jts,jte
163      do i = its,ite
164        cu_act_flag(i,j)=.TRUE.
165      enddo
166    enddo
167    delt=dt*stepcu
168    rdelt=1./delt
170    kbmax = kte
171    kbm   = kte
172    kmax  = kte
173    do k = kts,kte
174      if(znu(k).gt.0.45) kbmax = k + 1
175      if(znu(k).gt.0.70) kbm   = k + 1
176      if(znu(k).gt.0.05) kmax  = k + 1
177    enddo
179 ! outer most J_loop
181    do j = jts,jte
182      do k = kts,kte
183        kp = k+1
184        do i = its,ite
185          dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
186          prsll(i,k)=p3d(i,k,j)*0.001
187          prsii(i,k)=p3di(i,k,j)*0.001
188        enddo
189      enddo
191      do i = its,ite
192        prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
193      enddo
195      do i = its,ite
196        zii(i,1)=0.0
197      enddo     
199      do k = kts,kte                                            
200        do i = its,ite
201          zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
202        enddo
203      enddo
205      do k = kts,kte                
206        do i = its,ite                                                  
207          zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
208        enddo                                                         
209      enddo
211      do k = kts,kte
212        do i = its,ite
213          del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
214          u1(i,k)=u3d(i,k,j)
215          v1(i,k)=v3d(i,k,j)
216          q1(i,k)=qv3d(i,k,j)
217 !        q1(i,k)=qv3d(i,k,j)/(1.+qv3d(i,k,j))
218          t1(i,k)=t3d(i,k,j)
219          qi2(i,k) = qi3d(i,k,j)
220          qc2(i,k) = qc3d(i,k,j)
221        enddo
222      enddo
224 ! NCEP SAS 
226      call nsas2d(delt=delt,delx=dx,del=del(its,kts),                           &
227               prsl=prsll(its,kts),prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),   &
228               zl=zll(its,kts),                                                 &
229               ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts),                 &
230               q1=q1(its,kts),t1=t1(its,kts),rain=rain(its),                    &
231               kbot=kbot(its),ktop=ktop(its),                                   &
232               icps=icps(its),                                                  &
233               lat=j,slimsk=xland(ims,j),dot=dot(its,kts),                      &
234               u1=u1(its,kts), v1=v1(its,kts),                                  &
235               cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv,                      &
236               rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2,                               &
237               cice=cice,xls=xls,psat=psat,                                     &
238               dx_factor_nsas=dx_factor_nsas,                                   &
239               hpbl=hpbl(ims,j),hpbl_hold=hpbl_hold(ims,j),                     &
240               kbmax=kbmax,kbm=kbm,kmax=kmax,                                   &
241               ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,               &
242               ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,               &
243               its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
245      do i = its,ite
246        pratec1(i,j)=rain(i)*1000./(stepcu*dt)
247        raincv1(i,j)=rain(i)*1000./(stepcu)
248      enddo
250      do i = its,ite
251        raincv(i,j) = raincv1(i,j)
252        pratec(i,j) = pratec1(i,j)
253        hbot(i,j) = kbot(i)
254        htop(i,j) = ktop(i)
255      enddo
257      IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
258        do k = kts,kte
259          do i = its,ite
260            rthcuten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
261            rqvcuten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
262          enddo
263        enddo
264      ENDIF
266      IF(PRESENT(rucuten).AND.PRESENT(rvcuten)) THEN
267        do k = kts,kte
268          do i = its,ite
269            rucuten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
270            rvcuten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
271          enddo
272        enddo
273      ENDIF
275      IF(PRESENT( rqicuten )) THEN
276        IF ( F_QI ) THEN
277          do k = kts,kte
278            do i = its,ite
279              rqicuten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
280            enddo
281          enddo
282        ENDIF
283      ENDIF
285      IF(PRESENT( rqccuten )) THEN
286        IF ( F_QC ) THEN
287          do k = kts,kte
288            do i = its,ite
289              rqccuten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
290            enddo
291          enddo
292        ENDIF
293      ENDIF
295    enddo ! outer most J_loop
297    return
298    end subroutine cu_ksas
300 !-------------------------------------------------------------------------------
301 ! NCEP SAS (Deep Convection Scheme)
302 !-------------------------------------------------------------------------------
303    subroutine nsas2d(delt,delx,del,prsl,prsi,prslk,zl,                         &
304             ncloud,                                                            &
305             qc2,qi2,                                                           &
306             q1,t1,rain,kbot,ktop,                                              &
307             icps,                                                              &
308             lat,slimsk,dot,u1,v1,cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2,     &
309             cice,xls,psat,                                                     &
310             dx_factor_nsas,                                                    &
311             hpbl,hpbl_hold,                                                    &
312             kbmax,kbm,kmax,                                                    &
313             ids,ide, jds,jde, kds,kde,                                         &
314             ims,ime, jms,jme, kms,kme,                                         &
315             its,ite, jts,jte, kts,kte)
316 !-------------------------------------------------------------------------------
318 ! subprogram:    phys_cps_sas      computes convective heating and moistening
319 !                                                      and momentum transport
321 ! abstract: computes convective heating and moistening using a one
322 !   cloud type arakawa-schubert convection scheme originally developed
323 !   by georg grell. the scheme has been revised at ncep since initial 
324 !   implementation in 1993. it includes updraft and downdraft effects.
325 !   the closure is the cloud work function. both updraft and downdraft
326 !   are assumed to be saturated and the heating and moistening are
327 !   accomplished by the compensating environment. convective momentum transport
328 !   is taken into account. the name comes from "simplified arakawa-schubert
329 !   convection parameterization".
331 ! developed by hua-lu pan, wan-shu wu, songyou hong, and jongil han
332 !   implemented into wrf by kyosun sunny lim and songyou hong
333 !   module with cpp-based options is available in grims 
335 ! program history log:
336 !   92-03-01  hua-lu pan       operational development
337 !   96-03-01  song-you hong    revised closure, and trigger 
338 !   99-03-01  hua-lu pan       multiple clouds
339 !   06-03-01  young-hwa byun   closure based on moisture convergence (optional)
340 !   09-10-01  jung-eun kim     f90 format with standard physics modules
341 !   10-07-01  jong-il han      revised cloud model,trigger, as in gfs july 2010
342 !   10-12-01  kyosun sunny lim wrf compatible version
343 !   14-01-09  song-you hong    dx dependent trigger, closure, and mass flux
344 !   15-02-26  song-you hong    negative qv generation suppressed
345 !   15-06-01  jongil han       gfs sas
346 !   15-06-07  song-you hong    scale-aware cps
347 !   15-06-28  song-you hong    diurnal cycle with cldwrk in pbl
348 !   15-10-01  ji-young han     bug fix (dz to dz1)
349 !   16-01-01  ji-young han     enhanced entrainment at lower RH
350 !   16-06-01  ji-young han     moisture-based trigger threshold
351 !   16-09-01  ji-young han     bug fix & eliminate inconsistency
352 !   16-09-23  ji-young han     code clean-up
353 !   16-11-01  ji-young han     revised pgcon & bug fix in vshear
354 !   17-02-23  ji-young han     revised xlamb
355 !   18-03-01  ji-young han     kim sas
356 !   19-03-01  ji-yeon jang     revised triggering ftn for higher resolution<10km
357 !             yong-hee lee
358 !             young cheol kwon
359 !             eun-jeong lee    bug fix in c0 for overshooting layer
360 !             seung-bu park    bug fix in xmb trigger
362 ! usage:    call phys_cps_sas(delt,delx,del,prsl,prsi,prslk,prsik,zl,          &
363 !                             q2,q1,t1,u1,v1,rcs,slimsk,dot,cldwrk,rain,       &
364 !                             jcap,ncloud,lat,kbot,ktop,icps,                  &
365 !                             ids,ide, jds,jde, kds,kde,                       &
366 !                             ims,ime, jms,jme, kms,kme,                       &
367 !                             its,ite, jts,jte, kts,kte)
369 !   delt     - real model integration time step
370 !   delx     - real model grid interval
371 !   del      - real (kms:kme) sigma layer thickness
372 !   prsl     - real (ims:ime,kms:kme) pressure values
373 !   prsi     - real (ims:ime,kms:kme) pressure values at interface level
374 !   prslk    - real (ims:ime,kms:kme) pressure values to the kappa
375 !   prsik    - real (ims:ime,kms:kme) pressure values to the kappa at interface lev.
376 !   zl       - real (ims:ime,kms:kme) height above sea level
377 !   zi       - real (ims:ime,kms:kme) height above sea level at interface level
378 !   rcs      - real
379 !   slimsk   - real (ims:ime) land(1),sea(0), ice(2) flag
380 !   dot      - real (ims:ime,kms:kme) vertical velocity
381 !   jcap     - integer wave number 
382 !   ncloud   - integer no_cloud(0),no_ice(1),cloud+ice(2) 
383 !   lat      - integer  current latitude index
385 ! output argument list:
386 !   q2       - real (ims:ime,kms:kme) detrained hydrometeors in kg/kg
387 !            - in case of the  --> qc2(cloud), qi2(ice)
388 !   q1       - real (ims:ime,kms:kme) adjusted specific humidity in kg/kg
389 !   t1       - real (ims:ime,kms:kme) adjusted temperature in kelvin
390 !   u1       - real (ims:ime,kms:kme) adjusted zonal wind in m/s
391 !   v1       - real (ims:ime,kms:kme) adjusted meridional wind in m/s
392 !   cldwrk   - real (ims:ime) cloud work function
393 !   rain     - real (ims:ime) convective rain in meters
394 !   kbot     - integer (ims:ime) cloud bottom level
395 !   ktop     - integer (ims:ime) cloud top level
396 !   icps     - integer (ims:ime) bit flag indicating deep convection
398 ! subprograms called:
399 !   fpvs     - function to compute saturation vapor pressure
401 ! remarks: function fpvs is inlined by fpp.
402 !          nonstandard automatic arrays are used.
404 ! references :
405 !   pan and wu    (1995, ncep office note)
406 !   hong and pan  (1998, mon wea rev)
407 !   park and hong (2007,jmsj)
408 !   byun and hong (2007, mon wea rev)
409 !   han and pan   (2011, wea. forecasting)
410 !   lim et al.    (2014, wea. forecasting)
411 !   han et al.    (2016, mon wea rev)
412 !   kwon and hong (2017, mon wea rev)
413 !   han and hong  (2018, wea. forecasting)
415 !-------------------------------------------------------------------------------
416    implicit none
417 !-------------------------------------------------------------------------------
419 ! model tunable parameters 
421    real,parameter  ::  betal  = 0.05,   betas  = 0.05
422    real,parameter  ::  c0     = 0.002,  c1     = 0.002
423    real,parameter  ::  xlamdd = 1.0e-4, xlamde = 1.0e-4
424    real,parameter  ::  clam   = 0.1,    cxlamu = 1.0e-3
425    real,parameter  ::  aafac  = 0.1
426    real,parameter  ::  dthk   = 25.
427    real,parameter  ::  cinpcrmx = 240.,cinpcrmn = 120.
428    real,parameter  ::  cinacrmx = -120.
429    real,parameter  ::  bet1 = 1.875, cd1 = 0.506
430    real,parameter  ::  f1   = 2.0,  gam1 = 0.5,  tfac = 1.0
431    real,parameter  ::  dx1km = 1000., dx5km = 5000., dx250m = 250.
432    real,parameter  ::  edtmaxl = 0.3, edtmaxs = 0.3
433    real,parameter  ::  evfacts = 0.3, evfactl = 0.3
434    real,parameter  ::  tf=233.16,tcr=273.16,tcrf=1.0/(tcr-tf)
436 !  passing variables
438    real            ::  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
439    real            ::  pi_,qmin_,t0c_,cice,xlv0,xls,psat
440    integer         ::  dx_factor_nsas
441    integer         ::  lat,                                                    &
442                        ncloud,                                                 &
443                        ids,ide, jds,jde, kds,kde,                              &
444                        ims,ime, jms,jme, kms,kme,                              &
445                        its,ite, jts,jte, kts,kte
447    real            ::  delt,rcs
448    real            ::  del(its:ite,kts:kte),                                   &
449                        prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme),           &
450                        prsi(its:ite,kts:kte+1),                                &
451                        zl(its:ite,kts:kte),                                    &
452                        q1(its:ite,kts:kte),t1(its:ite,kts:kte),                &
453                        u1(its:ite,kts:kte),v1(its:ite,kts:kte),                &
454                        qci(its:ite,kts:kte),qrs(its:ite,kts:kte),              &
455                        dot(its:ite,kts:kte)
456    real            ::  qi2(its:ite,kts:kte)
457    real            ::  qc2(its:ite,kts:kte)
459    real            ::  rain(its:ite)
460    real            ::  hpbl(ims:ime),hpbl_hold(ims:ime)
461    integer         ::  kbot(its:ite),ktop(its:ite),icps(its:ite)
462    real            ::  slimsk(ims:ime)
464 !  local variables and arrays
466    integer         ::  i,k,kmax,kbmax,kbm,jmn,indx,     kts1,kte1,kmax1,kk
467    real            ::  p(its:ite,kts:kte),pdot(its:ite),acrtfct(its:ite)
468    real            ::  pden(its:ite)
469    real            ::  zi(its:ite,kts:kte+1)
470    real            ::  uo(its:ite,kts:kte),vo(its:ite,kts:kte)
471    real            ::  to(its:ite,kts:kte),qo(its:ite,kts:kte)
472    real            ::  hcko(its:ite,kts:kte)
473    real            ::  qcko(its:ite,kts:kte),eta(its:ite,kts:kte)
474    real            ::  etad(its:ite,kts:kte)
475    real            ::  qrcdo(its:ite,kts:kte)
476    real            ::  pwo(its:ite,kts:kte),pwdo(its:ite,kts:kte)
477    real            ::  c0t(its:ite,kts:kte)
478    real            ::  c1t(its:ite,kts:kte)
479    real            ::  bb1, bb2, wucb
480    real            ::  cinpcr
481    real            ::  cinacr
483 ! for updraft velocity calculation
485    real            ::  po1(its:ite,kts:kte),wu2(its:ite,kts:kte),              &
486                        buo(its:ite,kts:kte),drag(its:ite,kts:kte)
487    real            ::  wbar(its:ite),wc(its:ite),clear(its:ite)
489    real            ::  sigma, sigma_con
490    real            ::  mbdt
491    real            ::  dtconv(its:ite)
492    real            ::  deltv(its:ite),acrt(its:ite)
493    real            ::  apbl(its:ite),dtpbl(its:ite)
494    real            ::  qeso(its:ite,kts:kte)
495    real            ::  tvo(its:ite,kts:kte),dbyo(its:ite,kts:kte)
496    real            ::  heo(its:ite,kts:kte),heso(its:ite,kts:kte)
497    real            ::  qrcd(its:ite,kts:kte)
498    real            ::  dellah(its:ite,kts:kte),dellaq(its:ite,kts:kte)
500    integer         ::  kb(its:ite),kbcon(its:ite)
501    integer         ::  kbcon1(its:ite)
502    real            ::  hmax(its:ite),delq(its:ite)
503    real            ::  hkbo(its:ite),qkbo(its:ite)
504    integer         ::  lmin(its:ite),jmin(its:ite)
505    integer         ::  ktcon(its:ite)
506    integer         ::  ktcon1(its:ite)
507    integer         ::  kbdtr(its:ite)
508    real            ::  hmin(its:ite),pwavo(its:ite)
509    real            ::  aa1(its:ite),vshear(its:ite)
510    real            ::  qevap(its:ite)
511    real            ::  edt(its:ite)
512    real            ::  edt_s(its:ite)
513    real            ::  edto(its:ite),pwevo(its:ite)
514    real            ::  qcond(its:ite)
515    real            ::  hcdo(its:ite,kts:kte)
516    real            ::  qcdo(its:ite,kts:kte)
517    real            ::  xhkb(its:ite),xqkb(its:ite)
518    real            ::  xpwav(its:ite),xpwev(its:ite),xhcd(its:ite,kts:kte)
519    real            ::  xaa0(its:ite),f(its:ite),xk(its:ite)
520    real            ::  xmb(its:ite)
521    real            ::  edtx(its:ite),xqcd(its:ite,kts:kte)
522    real            ::  hsbar(its:ite),xmbmax(its:ite)
523    real            ::  xlamb(its:ite,kts:kte),xlamd(its:ite)
524    real            ::  cina(its:ite)
525    real            ::  delhbar(its:ite),delqbar(its:ite),deltbar(its:ite)
526    real            ::  qcirs(its:ite,kts:kte)
527    real            ::  dellal(its:ite,kts:kte)
528    real            ::  rntot(its:ite),delqev(its:ite),delq2(its:ite) 
530    real            ::  fent1(its:ite,kts:kte),fent2(its:ite,kts:kte)
531    real            ::  frh(its:ite,kts:kte)
532    real            ::  xlamud(its:ite),sumx(its:ite)
533    real            ::  frh_sum(its:ite),cinpcri(its:ite)
534    real            ::  aa2(its:ite)
535    real            ::  ucko(its:ite,kts:kte),vcko(its:ite,kts:kte)
536    real            ::  ucdo(its:ite,kts:kte),vcdo(its:ite,kts:kte)
537    real            ::  dellau(its:ite,kts:kte),dellav(its:ite,kts:kte)
538    real            ::  delubar(its:ite),delvbar(its:ite)
539    real            ::  qlko_ktcon(its:ite)
541    real            ::  alpha,beta,                                             &
542                        dt2,dtmin,dtmax,                                        &
543                        el2orc,eps,fact1,fact2,                                 &
544                        tem,tem1
545    real            ::  dz,dp,es,pprime,qs,                                     &
546                        dqsdp,desdt,dqsdt,gamma,                                &
547                             c0fac,alpha1,beta1,ccn_f,                          &
548                        dt,dq,po,     delx,                                     &
549                        factor,onemf,dz1,qrch,etah,qlk,qc,rfact,shear,          &
550                        e1,dh,                 edtmax,dhh,dg,aup,adw,           &
551                        dv1,dv2,dv3,dv1q,dv2q,dv3q,                             &
552                        dv1u,dv2u,dv3u,dv1v,dv2v,dv3v,                          &
553                        dellat,xdby,xqrch,    xpw,xpwd,                         &
554                        qrsk(its:ite,kts:kte),evef,ptem,ptem1
556    logical         ::  totflg, cnvflg(its:ite),flg(its:ite)
557    real            ::  pgcon(its:ite,kts:kte)
559 !-----------------------------------------------------------------------
561 ! define miscellaneous values
563    pi_   = 3.14159
564    qmin_ = 1.0e-30
565    t0c_ = 273.15
566    xlv0 = hvap_
567    rcs  = 1.
568    el2orc = hvap_*hvap_/(rv_*cp_)
569    eps    = rd_/rv_
570    fact1  = (cvap_-cliq_)/rv_
571    fact2  = hvap_/rv_-fact1*t0c_
572    kts1 = kts + 1
573    kte1 = kte - 1
574    dt2    = delt
575    dtmin  = max(dt2,600.)
576    dtmax  = max(dt2,10800.)
577    mbdt   = dt2
578    sigma_con = tan(0.4*pi_)/(dx5km-dx1km)                     ! 7.7 e-4 m-1
579    sigma  = (1.-1./pi_*(atan(sigma_con*(delx-dx5km))+pi_/2.)) ! 1(1km),0.1(10km)
581    if (delx.lt.dx5km) then
582      sigma = min(sigma - 0.01684 * delx/1000. + 0.0842, 1.0)
583    endif
585    cinpcr = (cinpcrmn + 0.5*(cinpcrmx-cinpcrmn)) * (1.-sigma)
587 !  initialize arrays
589    do i = its,ite
590      rain(i)    = 0.0
591      kbot(i)    = kte+1
592      ktop(i)    = 0
593      icps(i)    = 0
594      cnvflg(i)  = .true.
595      dtconv(i)  = 3600.
596      pdot(i)    = 0.0
597      edto(i)    = 0.0
598      edtx(i)    = 0.0
599      xmbmax(i)  = 0.3
600      aa2(i)     = 0.0
601      qlko_ktcon(i) = 0.0
602      lmin(i)    = 1
603      jmin(i)    = 1
604      edt(i)     = 0.0
605      cina(i)    = 0.0
606      frh_sum(i) = 0.0
607      cinpcri(i) = 0.0
608      apbl(i)    = 0.0
609      dtpbl(i)   = 0.0
610      do k = kts,kte
611        pgcon(i,k) = 0.5
612      enddo
613    enddo
615 ! Define top layer for search of the downdraft originating layer
616 ! and the maximum thetae for updraft
618    kmax = min(kmax,kte)
619    kmax1 = kmax - 1
620    kbm = min(kbm,kte)
622 ! convert surface pressure to mb from cb
624    do k = kts,kte
625      do i = its,ite
626        c0t(i,k)  = c0
627        c1t(i,k)  = c1 * sigma
628        qcirs(i,k)= 0.0
629        qci(i,k)  = 0.0
630        qrs(i,k)  = 0.0
631        qrsk(i,k) = 0.0
632        wu2(i,k)  = 0.0
633        buo(i,k)  = 0.0
634        drag(i,k) = 0.0
635        pwo(i,k)  = 0.0
636        pwdo(i,k) = 0.0
637        dellal(i,k) = 0.0
638        hcko(i,k) = 0.0
639        qcko(i,k) = 0.0
640        hcdo(i,k) = 0.0
641        qcdo(i,k) = 0.0
642      enddo
643    enddo
645    do k = kts,kmax
646      do i = its,ite
647        p(i,k) = prsl(i,k) * 10.
648        to(i,k) = t1(i,k)
649        qo(i,k) = q1(i,k)
650        dbyo(i,k) = 0.0
651        fent1(i,k) = 1.0
652        fent2(i,k) = 1.0
653        frh(i,k) = 0.0
654        ucko(i,k) = 0.0
655        vcko(i,k) = 0.0
656        ucdo(i,k) = 0.0
657        vcdo(i,k) = 0.0
658        uo(i,k) = u1(i,k) * rcs
659        vo(i,k) = v1(i,k) * rcs
660      enddo
661    enddo
663 ! column variables
665 !  p is pressure of the layer (mb)
666 !  t is temperature at t-dt (k)..tn
667 !  q is mixing ratio at t-dt (kg/kg)..qn
668 !  to is temperature at t+dt (k)... this is after advection and turbulan
669 !  qo is mixing ratio at t+dt (kg/kg)..q1
671    do k = kts,kmax
672      do i = its,ite
673        qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
674        qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
675        qeso(i,k) = max(qeso(i,k),qmin_)
676        qo(i,k)   = max(qo(i,k), 1.e-10 )
677 !      tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
678      enddo
679    enddo
681 ! compute moist static energy
683    do k = kts,kmax
684      do i = its,ite
685        heo(i,k)  = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qo(i,k)
686        heso(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qeso(i,k)
687      enddo
688    enddo
690 ! Determine level with largest moist static energy
691 ! This is the level where updraft starts
693    do i = its,ite
694      hmax(i) = heo(i,1)
695      kb(i) = 1
696    enddo
698    do k = kts1,kbm
699      do i = its,ite
700        if(heo(i,k).gt.hmax(i)) then
701          kb(i) = k
702          hmax(i) = heo(i,k)
703        endif
704      enddo
705    enddo
707    do i = its,ite
708      if (qo(i,kb(i)).lt.qmin_) cnvflg(i) = .false.
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
770    do i = its,ite
771      if(kbcon(i).eq.kmax) cnvflg(i) = .false.
772    enddo
774    totflg = .true.
775    do i = its,ite
776      totflg = totflg .and. (.not. cnvflg(i))
777    enddo
778    if(totflg) return
780    do k = kts,kmax1
781      do i = its,ite
782        if (cnvflg(i)) then
783          if (k.ge.kb(i).and.k.le.kbcon(i)) then
784            frh_sum(i) = frh_sum(i) + (1-frh(i,k))
785          endif
786        endif
787      enddo
788    enddo
790    do i = its,ite
791      if(cnvflg(i)) then
792        tem1 = p(i,kb(i)) - p(i,kbcon(i))
793        cinpcri(i) = cinpcr * frh_sum(i)/(kbcon(i)-kb(i)+1)
794        if (tem1.gt.cinpcri(i)) cnvflg(i) = .false.
795      endif
796    enddo
798    totflg = .true.
799    do i = its,ite
800      totflg = totflg .and. (.not. cnvflg(i))
801    enddo
802    if(totflg) return
804    do k = kts1,kte
805      do i = its,ite
806        zi(i,k) = 0.5*(zl(i,k-1)+zl(i,k))
807      enddo
808    enddo
810    do k = kts,kte1
811      do i = its,ite
812        xlamb(i,k) = clam / zi(i,k+1) 
813      enddo
814    enddo
816 !  assume that updraft entrainment rate above cloud base is
817 !  same as that at cloud base
819    do k = kts1,kmax1
820      do i = its,ite
821        if(cnvflg(i).and.(k.gt.kbcon(i))) then
822          xlamb(i,k) = xlamb(i,kbcon(i))
823        endif
824      enddo
825    enddo
827 !  assume the detrainment rate for the updrafts to be same as
828 !  the entrainment rate at cloud base
830    do i = its,ite
831      if(cnvflg(i)) then
832        xlamud(i) = xlamb(i,kbcon(i))
833      endif
834    enddo
836 !  functions rapidly decreasing with height, mimicking a cloud ensemble
837 !    (Bechtold et al., 2008)
839    do k = kts1,kmax1
840      do i = its,ite
841        if(cnvflg(i).and.(k.gt.kbcon(i))) then
842          tem = qeso(i,k)/qeso(i,kbcon(i))
843          fent1(i,k) = tem**2
844          fent2(i,k) = tem**3
845        endif
846      enddo
847    enddo
849 !  final entrainment rate as the sum of turbulent part and organized entrainment
850 !    depending on the environmental relative humidity
851 !    (Bechtold et al., 2008)
853    do k = kts1,kmax1
854      do i = its,ite
855        if(cnvflg(i).and.(k.ge.kbcon(i))) then
856           tem = cxlamu * frh(i,k) * fent2(i,k)
857           xlamb(i,k) = xlamb(i,k)*fent1(i,k) + tem
858        endif
859      enddo
860    enddo
862 ! determine updraft mass flux
864    do k = kts,kte
865      do i = its,ite
866       if(cnvflg(i)) then
867          eta(i,k) = 1.
868        endif
869      enddo
870    enddo
872    do k = kbmax,kts1,-1
873      do i = its,ite
874        if(cnvflg(i).and.k.lt.kbcon(i).and.k.ge.kb(i)) then
875          dz = zi(i,k+2) - zi(i,k+1)
876          ptem     = 0.5*(xlamb(i,k)+xlamb(i,k+1))-xlamud(i)
877          eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
878        endif
879      enddo
880    enddo
881    do k = kts1,kmax1
882      do i = its,ite
883        if(cnvflg(i).and.k.gt.kbcon(i)) then
884          dz  = zi(i,k+1) - zi(i,k)
885          ptem     = 0.5*(xlamb(i,k)+xlamb(i,k-1))-xlamud(i)
886          eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
887        endif
888      enddo
889    enddo
890    do i = its,ite
891      if(cnvflg(i)) then
892        dz = zi(i,3) - zi(i,2)
893        ptem     = 0.5*(xlamb(i,1)+xlamb(i,2))-xlamud(i)
894        eta(i,1) = eta(i,2) / (1. + ptem * dz)
895      endif
896    enddo
898 ! work up updraft cloud properties
900    do i = its,ite
901      if(cnvflg(i)) then
902        indx = kb(i)
903        hcko(i,indx) = hkbo(i)
904        qcko(i,indx) = qkbo(i)
905        ucko(i,indx) = uo(i,indx)
906        vcko(i,indx) = vo(i,indx)
907        pwavo(i) = 0.
908      endif
909    enddo
911 ! cloud property below cloud base is modified by the entrainment proces
913    do k = kts1,kmax1
914      do i = its,ite
915        if(cnvflg(i).and.k.gt.kb(i)) then
916          dz   = zi(i,k+1) - zi(i,k)
917          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
918          tem1 = 0.5 * xlamud(i) * dz
919          factor = 1. + tem - tem1
920          hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                           &
921                      (heo(i,k)+heo(i,k-1)))/factor
922          dbyo(i,k) = hcko(i,k) - heso(i,k)
923        endif
924      enddo
925    enddo
927 !   taking account into convection inhibition due to existence of
928 !    dry layers below cloud base
930    do i = its,ite
931      flg(i) = cnvflg(i)
932      kbcon1(i) = kmax
933    enddo
935    do k = kts1,kmax
936      do i = its,ite
937        if(flg(i).and.k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
938          kbcon1(i) = k
939          flg(i) = .false.
940        endif
941      enddo
942    enddo
944    do i = its,ite
945      if(cnvflg(i)) then
946        if(kbcon1(i).eq.kmax) cnvflg(i) = .false.
947      endif
948    enddo
950    do i = its,ite
951      if(cnvflg(i)) then
952        tem = p(i,kbcon(i)) - p(i,kbcon1(i))
953        if(tem.gt.dthk) then
954           cnvflg(i) = .false.
955        endif
956      endif
957    enddo
959    totflg = .true.
960    do i = its,ite
961      totflg = totflg .and. (.not. cnvflg(i))
962    enddo
963    if(totflg) return
965 ! calculate convective inhibition
967    do k = kts1,kmax1
968      do i = its,ite
969        if (cnvflg(i)) then
970          if (k.gt.kb(i).and.k.lt.kbcon1(i)) then
971            dz1 = (zi(i,k+1) - zi(i,k))
972            gamma = el2orc * qeso(i,k) / (to(i,k)**2)
973            rfact =  1. + fv_ * cp_ * gamma * to(i,k) / hvap_
974            cina(i) = cina(i) + dz1 * (g_ / (cp_ * to(i,k)))                    &
975                    * dbyo(i,k) / (1. + gamma) * rfact
976            cina(i) = cina(i) + dz1 * g_ * fv_ * max(0.,(qeso(i,k) - qo(i,k)))
977          endif
978        endif
979      enddo
980    enddo
982    do i = its,ite
983      if (cnvflg(i)) then
984        cinacr = cinacrmx
985        if (cina(i).lt.cinacr) cnvflg(i) = .false.
986      endif
987    enddo
989    totflg = .true.
990    do i = its,ite
991      totflg = totflg .and. (.not. cnvflg(i))
992    enddo
993    if (totflg) return
995 !  determine cloud top
997    do i = its,ite
998      flg(i) = cnvflg(i)
999      ktcon(i) = 1
1000    enddo
1002 !   check inversion
1004    do k = kts1,kmax1
1005      do i = its,ite
1006        if(dbyo(i,k).lt.0..and.flg(i).and.k.gt. kbcon1(i)) then
1007          ktcon(i) = k
1008          flg(i)   = .false.
1009        endif
1010      enddo
1011    enddo
1013 ! check cloud depth
1015    do i = its,ite
1016      if(cnvflg(i).and.(p(i,kbcon(i)) - p(i,ktcon(i))).lt.150.)                 &
1017      cnvflg(i) = .false.
1018    enddo
1020    totflg = .true.
1021    do i = its,ite
1022      totflg = totflg .and. (.not. cnvflg(i))
1023    enddo
1024    if(totflg) return
1026 !  search for downdraft originating level above theta-e minimum
1028    do i = its,ite 
1029      if(cnvflg(i)) then
1030        hmin(i) = heo(i,kbcon1(i))
1031        lmin(i) = kbmax
1032        jmin(i) = kbmax
1033     endif
1034    enddo
1036    do k = kts1,kbmax 
1037      do i = its,ite 
1038        if(cnvflg(i).and.k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
1039          lmin(i) = k + 1
1040          hmin(i) = heo(i,k)
1041        endif
1042      enddo
1043    enddo
1045 ! make sure that jmin is within the cloud
1047    do i = its,ite
1048      if(cnvflg(i)) then
1049        jmin(i) = min(lmin(i),ktcon(i)-1)
1050        jmin(i) = max(jmin(i),kbcon1(i)+1)
1051        if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
1052        if(jmin(i).le.kbcon(i)) cnvflg(i) = .false.
1053      endif
1054    enddo
1056 !  specify upper limit of mass flux at cloud base
1058    do i = its,ite
1059      if(cnvflg(i)) then
1060        k = kbcon(i)
1061        dp = 1000. * del(i,k)
1062        xmbmax(i) = dp / (g_ * dt2)
1063      endif
1064    enddo
1066 ! compute cloud moisture property and precipitation
1068    do k = kts1,kmax
1069      do i = its,ite
1070        if (cnvflg(i).and.k.gt.kb(i)) then
1071          alpha1 = min((-0.7*log(100.)+24.)*0.0001,c0)
1072          beta1 = 0.07
1074          if (to(i,k).gt.t0c_) then
1075            c0fac = alpha1
1076          else
1077            c0fac = alpha1*exp(beta1*(to(i,k)-t0c_))
1078          endif
1080          c0fac = max(0.0,c0fac)
1081          c0t(i,k) = c0fac
1082        endif
1084        if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1085          dz1 = (zi(i,k+1) - zi(i,k))
1086          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1087          qrch = qeso(i,k)                                                      &
1088               + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1089          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz1
1090          tem1 = 0.5 * xlamud(i) * dz1
1091          factor = 1. + tem - tem1
1092          qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                           &
1093                     (qo(i,k)+qo(i,k-1)))/factor
1094          qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1096 ! check if there is excess moisture to release latent heat
1098          if(qcirs(i,k).gt.0. .and. k.ge.kbcon(i)) then
1099            etah = .5 * (eta(i,k) + eta(i,k-1))
1100            if(ncloud.gt.0..and.k.gt.jmin(i)) then
1101              dp = 1000. * del(i,k)
1102              ptem = c0t(i,k) + c1t(i,k)
1103              qlk = qcirs(i,k) / (eta(i,k) + etah * ptem * dz1)
1104              dellal(i,k) = etah * c1t(i,k) * dz1 * qlk * g_ / dp
1105            else
1106              qlk = qcirs(i,k) / (eta(i,k) + etah * c0t(i,k) * dz1)
1107            endif
1108            pwo(i,k) = etah * c0t(i,k) * dz1 * qlk
1109            qc = qlk + qrch
1110            qcko(i,k) = qc
1111            pwavo(i) = pwavo(i) + pwo(i,k)
1112            buo(i,k) = buo(i,k) - g_ * qlk
1114 ! compute buoyancy and drag for updraft velocity
1116            if (k.ge.kbcon(i)) then
1117              rfact = 1. + fv_ * cp_ * gamma                                    &
1118                      * to(i,k) / hvap_
1119              buo(i,k) = buo(i,k) + (g_ / (cp_ * to(i,k)))                      &
1120                       * dbyo(i,k) / (1. + gamma)                               &
1121                       * rfact
1122              buo(i,k) = buo(i,k) + g_ * fv_ *                                  &
1123                         max(0.,(qeso(i,k) - qo(i,k)))
1124              drag(i,k) = max(xlamb(i,k),xlamud(i))
1125            endif
1126          endif
1127        endif
1128      enddo
1129    enddo
1131 ! calculate cloud work function at t+dt
1133    do i = its,ite
1134      if (cnvflg(i)) then
1135        aa1(i) = 0.
1136      endif
1137    enddo
1139    do k = kts1,kmax
1140      do i = its,ite
1141        if (cnvflg(i)) then
1142          if (k.ge.kbcon(i) .and. k.lt.ktcon(i)) then
1143            dz1 = zl(i,k+1) - zl(i,k)
1144            aa1(i) = aa1(i) + buo(i,k) * dz1
1145          endif
1146        endif
1147      enddo
1148    enddo
1150    do i = its,ite
1151      if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1152    enddo
1154    totflg = .true.
1155    do i = its,ite
1156      totflg = totflg .and. (.not. cnvflg(i))
1157    enddo
1158    if(totflg) return
1160 !    estimate the convective overshooting as the level
1161 !    where the [aafac * cloud work function] becomes zero,
1162 !    which is the final cloud top
1164    do i = its,ite
1165      if (cnvflg(i)) then
1166        aa2(i) = aafac * aa1(i)
1167      endif
1168    enddo
1170    do i = its,ite
1171      flg(i) = cnvflg(i)
1172      ktcon1(i) = kmax1
1173    enddo
1175    do k = kts1,kmax
1176      do i = its, ite
1177        if (flg(i)) then
1178          if(k.ge.ktcon(i).and.k.lt.kmax) then
1179            dz1 = zl(i,k+1) - zl(i,k)
1180            gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1181            rfact =  1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1182            aa2(i) = aa2(i) +dz1 * (g_ / (cp_ * to(i,k)))                       &
1183                        * dbyo(i,k) / (1. + gamma)* rfact
1184            if(aa2(i).lt.0.) then
1185              ktcon1(i) = k
1186              flg(i) = .false.
1187            endif
1188          endif
1189        endif
1190      enddo
1191    enddo
1193 !  compute cloud moisture property, detraining cloud water
1194 !  and precipitation in overshooting layers
1196    do k = kts1,kmax
1197      do i = its,ite
1198        if (cnvflg(i)) then
1199          if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1200            dz = (zi(i,k+1) - zi(i,k))
1201            gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1202            qrch = qeso(i,k)+ gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1203            tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1204            tem1 = 0.5 * xlamud(i) * dz
1205            factor = 1. + tem - tem1
1206            qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                         &
1207                       (qo(i,k)+qo(i,k-1)))/factor
1208            qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1210 !  check if there is excess moisture to release latent heat
1212            if(qcirs(i,k).gt.0.) then
1213              etah = .5 * (eta(i,k) + eta(i,k-1))
1214              if(ncloud.gt.0.) then
1215                dp = 1000. * del(i,k)
1216                qlk = qcirs(i,k) / (eta(i,k) + etah * (c0t(i,k) + c1t(i,k)) * dz)
1217                dellal(i,k) = etah * c1t(i,k) * dz * qlk * g_ / dp
1218              else
1219                qlk = qcirs(i,k) / (eta(i,k) + etah * c0t(i,k) * dz)
1220              endif
1221              pwo(i,k) = etah * c0t(i,k) * dz * qlk
1222              qc = qlk + qrch
1223              qcko(i,k) = qc
1224              pwavo(i) = pwavo(i) + pwo(i,k)
1225            endif
1226          endif
1227        endif
1228      enddo
1229    enddo
1231 ! compute updraft velocity square(wu2)
1233    bb1 = 2. * (1.+bet1*cd1)
1234    bb2 = 2. / (f1*(1.+gam1))
1236 !  bb1 = 12.0
1237 !  bb2 = 0.67
1239    do i = its,ite
1240      if (cnvflg(i)) then
1241        k = kbcon1(i)
1242        po = .5 * (p(i,k) + p(i,k+1))
1243        tem = po / (rd_ * to(i,k))
1244        wucb = -10.*dot(i,k) / (tem * g_)
1245        if (wucb.gt.0.) then
1246          wu2(i,k) = wucb * wucb
1247        else
1248          wu2(i,k) = 0.
1249        endif
1250      endif
1251    enddo
1253    do k = kts1,kmax
1254      do i = its,ite
1255        if (cnvflg(i)) then
1256          if (k.gt.kbcon1(i) .and. k.lt.ktcon(i)) then
1257            dz   = zi(i,k+1) - zi(i,k)
1258            tem  = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
1259            tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
1260            ptem = (1. - tem) * wu2(i,k-1)
1261            ptem1 = 1. + tem
1262            wu2(i,k) = (ptem + tem1) / ptem1
1263            wu2(i,k) = max(wu2(i,k), 0.)
1264          endif
1265        endif
1266      enddo
1267    enddo
1269 ! compute mean updraft velocity and mean grid-scale vertical velocity
1271    wc = 0. ; wbar = 0. ; sumx = 0.
1273    ptem = -0.5 * rd_ / g_
1274    do k = kts1,kmax
1275      do i = its,ite
1276        po1(i,k) = .5 * (p(i,k) + p(i,k+1))
1277      enddo
1278    enddo
1280    do k = kts1,kmax
1281      do i = its,ite
1282        if (cnvflg(i)) then
1283          if (k.gt.kbcon1(i) .and. k.lt.ktcon(i)) then
1284            dz = zi(i,k+1) - zi(i,k)
1285            tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1286            wc(i) = wc(i) + tem * dz
1287            tem  = 10. * dot(i,k)   * to(i,k)   / po1(i,k)
1288            tem1 = 10. * dot(i,k-1) * to(i,k-1) / po1(i,k-1)
1289            wbar(i) = wbar(i) + ptem * (tem + tem1) * dz
1290            sumx(i) = sumx(i) + dz
1291          endif
1292        endif
1293      enddo
1294    enddo
1296    do i = its,ite
1297      if (cnvflg(i)) then
1298        if (sumx(i) == 0.) then
1299          cnvflg(i) = .false.
1300        else
1301          wc(i) = wc(i) / sumx(i)
1302          wbar(i) = wbar(i) / sumx(i)
1303        endif
1304        if (wc(i).lt.1.e-4) cnvflg(i) = .false.
1305      endif
1306    enddo
1308 ! compute mean cloud core fraction
1309 ! assume mean cloud core fraction to be the ratio of
1310 ! mean grid-scale vertical velocity (wbar) and mean updraft velocity
1312    do i = its,ite
1313      if (cnvflg(i)) then
1314        tem = wbar(i) / wc(i)
1315        tem = max(tem, 0.)
1316        clear(i) = 1. - tem
1317        clear(i) = max(min(clear(i), 1.0), 0.)
1318        if (wbar(i).gt.0. .and. wbar(i).gt.wc(i)) cnvflg(i) = .false.
1319      endif
1320    enddo
1322 ! exchange ktcon with ktcon1
1324    do i = its,ite
1325      if(cnvflg(i)) then
1326        kk = ktcon(i)
1327        ktcon(i) = ktcon1(i)
1328        ktcon1(i) = kk
1329      endif
1330    enddo
1332    do i = its,ite
1333      if (ktcon(i).le.kb(i)) cnvflg(i) = .false.
1334    enddo
1336    do k = kts,kte
1337      do i = its,ite
1338        if (cnvflg(i)) then
1339          if (k.le.ktcon(i)) then
1340            pgcon(i,k) = 0.5+0.5*exp(3.*(zl(i,k)-zl(i,ktcon(i)))/zl(i,ktcon(i)))
1341          else
1342            pgcon(i,k) = 1.0
1343          endif
1344        endif
1345        pgcon(i,k) = min(pgcon(i,k), 1.0)
1346        pgcon(i,k) = max(pgcon(i,k), 0.5)
1347      enddo
1348    enddo
1350    do k = kts1,kmax1
1351      do i = its,ite
1352        if (cnvflg(i) .and. k.gt.kb(i)) then
1353          dz   = zi(i,k+1) - zi(i,k)
1354          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1355          tem1 = 0.5 * xlamud(i) * dz
1356          factor = 1. + tem - tem1
1357          ptem  = 0.5 * tem + pgcon(i,k)
1358          ptem1 = 0.5 * tem - pgcon(i,k)
1359          ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)                       &
1360                      +ptem1*uo(i,k-1))/factor
1361          vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)                       &
1362                      +ptem1*vo(i,k-1))/factor
1363        endif
1364      enddo
1365    enddo
1367 ! this section is ready for cloud water
1369    if (ncloud.gt.0) then
1371 !  compute liquid and vapor separation at cloud top
1373      do i = its,ite
1374        if(cnvflg(i)) then
1375          k = ktcon(i)-1
1376          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1377          qrch = qeso(i,k)                                                      &
1378                 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1379          dq = qcko(i,k) - qrch
1381 !  check if there is excess moisture to release latent heat
1383          if(dq.gt.0.) then
1384            qlko_ktcon(i) = dq * sigma
1385            qcko(i,k) = qrch + dq * (1.-sigma)
1386          endif
1387        endif
1388      enddo
1389    endif
1391 ! ..... downdraft calculations .....
1393 ! determine downdraft strength in terms of wind shear
1395    do i = its,ite
1396      if(cnvflg(i)) then
1397        vshear(i) = 0.
1398      endif
1399    enddo
1401    do k = kts1,kmax
1402      do i = its,ite
1403        if(k.gt.kb(i).and.k.le.ktcon(i).and.cnvflg(i)) then
1404          shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2                                  &
1405                    + (vo(i,k)-vo(i,k-1)) ** 2)
1406          vshear(i) = vshear(i) + shear
1407        endif
1408      enddo
1409    enddo
1411    do i = its,ite
1412      if(cnvflg(i)) then
1413        vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1414        e1 = 1.591-.639*vshear(i)                                               &
1415            +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1416        edt(i)  = 1.-e1
1418        ccn_f = 1.0
1419        edt_s(i) = edt(i)
1420        edt_s(i) = min(edt_s(i),.9)
1421        edt_s(i) = max(edt_s(i),.0)
1422        edt(i)  = min(edt(i),.9)
1423        edt(i)  = max(edt(i),.0)
1424        edt(i)  = edt(i) * ccn_f
1425        edto(i) = edt(i)
1426        edtx(i) = edt(i)
1427      endif
1428    enddo
1430 ! determine detrainment rate between 1 and kbdtr
1432    do i = its,ite
1433      if(cnvflg(i)) then
1434        sumx(i) = 0.
1435      endif
1436    enddo
1438    do k = kts,kmax1
1439      do i = its,ite
1440        if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1441          dz = zi(i,k+2) - zi(i,k+1)
1442          sumx(i) = sumx(i) + dz
1443        endif
1444      enddo
1445    enddo
1447    do i = its,ite
1448      beta = betas
1449      if(slimsk(i).eq.1.) beta = betal
1450      if(cnvflg(i)) then
1451        kbdtr(i) = kbcon(i)
1452        kbdtr(i) = max(kbdtr(i),1)
1453        dz =(sumx(i)+zi(i,2))/float(kbcon(i))
1454        tem = 1./float(kbcon(i))
1455        xlamd(i) = (1.-beta**tem)/dz
1456      endif
1457    enddo
1459 ! determine downdraft mass flux
1461    do k = kts,kmax
1462      do i = its,ite
1463        if(cnvflg(i)) then
1464          etad(i,k) = 1.
1465        endif
1466        qrcdo(i,k) = 0.
1467        qrcd(i,k) = 0.
1468      enddo
1469    enddo
1471    do k = kmax1,kts,-1
1472      do i = its,ite
1473        if(cnvflg(i)) then
1474          if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
1475            dz = (zi(i,k+2) - zi(i,k+1))
1476            ptem = xlamdd-xlamde
1477            etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1478          elseif(k.lt.kbcon(i)) then
1479            dz = (zi(i,k+2) - zi(i,k+1))
1480            ptem = xlamd(i)+xlamdd-xlamde
1481            etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1482          endif
1483        endif
1484      enddo
1485    enddo
1487 ! downdraft moisture properties
1489    do i = its,ite
1490      if(cnvflg(i)) then
1491       pwevo(i) = 0.
1492      endif
1493    enddo
1495    do i = its,ite
1496      if(cnvflg(i)) then 
1497        jmn = jmin(i)
1498        hcdo(i,jmn) = heo(i,jmn)
1499        qcdo(i,jmn) = qo(i,jmn)
1500        qrcdo(i,jmn) = qeso(i,jmn)
1501        ucdo(i,jmn) = uo(i,jmn)
1502        vcdo(i,jmn) = vo(i,jmn)
1503      endif
1504    enddo
1506    do k = kmax1,kts,-1 
1507      do i = its,ite 
1508        if (cnvflg(i) .and. k.lt.jmin(i)) then
1509          dz = zi(i,k+2) - zi(i,k+1)
1510          if(k.ge.kbcon(i)) then
1511            tem  = xlamde * dz
1512            tem1 = 0.5 * xlamdd * dz
1513          else
1514            tem  = xlamde * dz
1515            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1516          endif
1517           factor = 1. + tem - tem1
1518           ptem  = 0.5 * tem - pgcon(i,k)
1519           ptem1 = 0.5 * tem + pgcon(i,k)
1520           hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*                          &
1521                       (heo(i,k)+heo(i,k+1)))/factor
1522           ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)                    &
1523                      +ptem1*uo(i,k))/factor
1524           vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)                    &
1525                      +ptem1*vo(i,k))/factor
1526           dbyo(i,k) = hcdo(i,k) - heso(i,k)
1527        endif
1528      enddo
1529    enddo
1531    do k = kmax1,kts,-1
1532      do i = its,ite
1533        if(cnvflg(i).and.k.lt.jmin(i)) then
1534          dq = qeso(i,k)
1535          dt = to(i,k)
1536          gamma = el2orc * dq / dt**2
1537          qrcdo(i,k) = dq+(1./hvap_)*(gamma/(1.+gamma))*dbyo(i,k)
1538          dz = zi(i,k+2) - zi(i,k+1)
1539          if(k.ge.kbcon(i)) then
1540            tem  = xlamde * dz
1541            tem1 = 0.5 * xlamdd * dz
1542          else
1543            tem  = xlamde * dz
1544            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1545          endif
1546          factor = 1. + tem - tem1
1547          qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5*                           &
1548                      (qo(i,k)+qo(i,k+1)))/factor
1549          pwdo(i,k) = etad(i,k+1) * qcdo(i,k) -etad(i,k+1) * qrcdo(i,k)
1550          qcdo(i,k) = qrcdo(i,k)
1551          pwevo(i) = pwevo(i) + pwdo(i,k)
1552        endif
1553      enddo
1554    enddo
1556 ! final downdraft strength dependent on precip
1557 ! efficiency (edt), normalized condensate (pwav), and
1558 ! evaporate (pwev)
1560    do i = its,ite
1561      edtmax = edtmaxl
1562      if(slimsk(i).eq.2.) edtmax = edtmaxs
1563      if(cnvflg(i)) then
1564        if(pwevo(i).lt.0.) then
1565          edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1566          edto(i) = min(edto(i),edtmax)
1567        else
1568          edto(i) = 0.
1569        endif
1570      endif
1571    enddo
1573 ! downdraft cloudwork functions
1575    do k = kmax1,kts,-1
1576      do i = its,ite
1577        if(cnvflg(i).and.k.lt.jmin(i)) then
1578          gamma = el2orc * qeso(i,k) / to(i,k)**2
1579          dhh = hcdo(i,k)
1580          dt = to(i,k)
1581          dg = gamma
1582          dh = heso(i,k)
1583          dz = -1.*(zl(i,k+1)-zl(i,k))
1584          aa1(i) = aa1(i)+edto(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))           &
1585                 *(1.+fv_*cp_*dg*dt/hvap_)
1586          aa1(i) = aa1(i)+edto(i)*dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1587        endif
1588      enddo
1589    enddo
1591    do i = its,ite
1592      if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1593    enddo
1595    totflg = .true.
1596    do i = its,ite
1597      totflg = totflg .and. (.not. cnvflg(i))
1598    enddo
1599    if(totflg) return
1601 ! what would the change be, that a cloud with unit mass
1602 ! will do to the environment?
1604    do k = kts,kmax
1605      do i = its,ite
1606        if(cnvflg(i)) then
1607          dellah(i,k) = 0.
1608          dellaq(i,k) = 0.
1609          dellau(i,k) = 0.
1610          dellav(i,k) = 0.
1611        endif
1612      enddo
1613    enddo
1615    do i = its,ite
1616      if(cnvflg(i)) then
1617        dp = 1000. * del(i,1)
1618        dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)                          &
1619                    - heo(i,1)) * g_ / dp
1620        dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1)                          &
1621                    - qo(i,1)) * g_ / dp
1622        dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)                          &
1623                    - uo(i,1)) * g_ / dp
1624        dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)                          &
1625                    - vo(i,1)) * g_ / dp
1626      endif
1627    enddo
1629 ! changed due to subsidence and entrainment
1631    do k = kts1,kmax1
1632      do i = its,ite
1633        if(cnvflg(i).and.k.lt.ktcon(i)) then
1634          aup = 1.
1635          if(k.le.kb(i)) aup = 0.
1636          adw = 1.
1637          if(k.gt.jmin(i)) adw = 0.
1638          dv1 = heo(i,k)
1639          dv2 = .5 * (heo(i,k) + heo(i,k-1))
1640          dv3 = heo(i,k-1)
1641          dv1q = qo(i,k)
1642          dv2q = .5 * (qo(i,k) + qo(i,k-1))
1643          dv3q = qo(i,k-1)
1644          dv1u = uo(i,k)
1645          dv2u = .5 * (uo(i,k) + uo(i,k-1))
1646          dv3u = uo(i,k-1)
1647          dv1v = vo(i,k)
1648          dv2v = .5 * (vo(i,k) + vo(i,k-1))
1649          dv3v = vo(i,k-1)
1650          dp = 1000. * del(i,k)
1651          dz = zi(i,k+1) - zi(i,k)
1652          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1))
1653          tem1 = xlamud(i)
1654          if(k.le.kbcon(i)) then
1655            ptem  = xlamde
1656            ptem1 = xlamd(i)+xlamdd
1657          else
1658            ptem  = xlamde
1659            ptem1 = xlamdd
1660          endif
1662          dellah(i,k) = dellah(i,k) +                                           &
1663              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1               &
1664          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3               &
1665          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2*dz              &
1666          +  aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz                  &
1667          +  adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz) *g_/dp
1669          dellaq(i,k) = dellaq(i,k) +                                           &
1670              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1q              &
1671          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3q              &
1672          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz             &
1673          +  aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz                  &
1674          +  adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz) *g_/dp
1676          dellau(i,k) = dellau(i,k) +                                           &
1677              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1u              &
1678          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3u              &
1679          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz             &
1680          +  aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz                  &
1681          +  adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz          &
1682          - pgcon(i,k)*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u))*g_/dp
1684          dellav(i,k) = dellav(i,k) +                                           &
1685              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1v              &
1686          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3v              &
1687          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz             &
1688          +  aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz                  &
1689          +  adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz          &
1690          - pgcon(i,k)*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v))*g_/dp
1691        endif
1692      enddo
1693    enddo
1695 ! cloud top
1697    do i = its,ite
1698      if(cnvflg(i)) then
1699        indx = ktcon(i)
1700        dp = 1000. * del(i,indx)
1701        dv1 = heo(i,indx-1)
1702        dellah(i,indx) = eta(i,indx-1) *                                        &
1703                         (hcko(i,indx-1) - dv1) * g_ / dp
1704        dv1q = qo(i,indx-1)
1705        dellaq(i,indx) = eta(i,indx-1) *                                        &
1706                         (qcko(i,indx-1) - dv1q) * g_ / dp
1707        dv1u = uo(i,indx-1)
1708        dellau(i,indx) = eta(i,indx-1) *                                        &
1709                         (ucko(i,indx-1) - dv1u) * g_ / dp
1710        dv1v = vo(i,indx-1)
1711        dellav(i,indx) = eta(i,indx-1) *                                        &
1712                         (vcko(i,indx-1) - dv1v) * g_ / dp
1714 !  cloud water
1716        dellal(i,indx) = eta(i,indx-1) * qlko_ktcon(i) * g_ / dp
1717      endif
1718    enddo
1720 ! final changed variable per unit mass flux
1722    do k = kts,kmax
1723      do i = its,ite
1724        if(cnvflg(i).and.k.gt.ktcon(i)) then
1725          qo(i,k) = q1(i,k)
1726          to(i,k) = t1(i,k)
1727        endif
1728        if(cnvflg(i).and.k.le.ktcon(i)) then
1729          qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1730          dellat  = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1731          to(i,k) = dellat * mbdt + t1(i,k)
1732          qo(i,k) = max(qo(i,k),1.0e-10)
1733        endif
1734      enddo
1735    enddo
1737 !------------------------------------------------------------------------------
1739 ! the above changed environment is now used to calulate the
1740 ! effect the arbitrary cloud (with unit mass flux)
1741 ! which then is used to calculate the real mass flux,
1742 ! necessary to keep this change in balance with the large-scale
1743 ! destabilization.
1745 ! environmental conditions again
1747 !------------------------------------------------------------------------------
1749    do k = kts,kmax
1750      do i = its,ite
1751        if(cnvflg(i)) then
1752          qeso(i,k)=0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1753          qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1754          qeso(i,k) = max(qeso(i,k),qmin_)
1755 !        tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
1756        endif
1757      enddo
1758    enddo
1760 ! moist static energy
1762    do k = kts,kmax1
1763      do i = its,ite
1764        if(cnvflg(i)) then
1765          dz = .5 * (zl(i,k+1) - zl(i,k))
1766          dp = .5 * (p(i,k+1) - p(i,k))
1767          es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1768          pprime = p(i,k+1) + (eps-1.) * es
1769          qs = eps * es / pprime
1770          dqsdp = - qs / pprime
1771          desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1772          dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
1773          gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1774          dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
1775          dq = dqsdt * dt + dqsdp * dp
1776          to(i,k) = to(i,k+1) + dt
1777          qo(i,k) = qo(i,k+1) + dq
1778          po = .5 * (p(i,k) + p(i,k+1))
1779          qeso(i,k) = 0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1780          qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
1781          qeso(i,k) = max(qeso(i,k),qmin_)
1782          qo(i,k)   = max(qo(i,k), 1.0e-10)
1783          heo(i,k)  = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                         &
1784                      cp_ * to(i,k) + hvap_ * qo(i,k)
1785          heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                         &
1786                      cp_ * to(i,k) + hvap_ * qeso(i,k)
1787        endif
1788      enddo
1789    enddo
1791    k = kmax
1792    do i = its,ite
1793      if(cnvflg(i)) then
1794        heo(i,k)  = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qo(i,k)
1795        heso(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qeso(i,k)
1796      endif
1797    enddo
1799    do i = its,ite
1800      if(cnvflg(i)) then
1801        xaa0(i) = 0.
1802        xpwav(i) = 0.
1803        indx = kb(i)
1804        xhkb(i) = heo(i,indx)
1805        xqkb(i) = qo(i,indx)
1806        hcko(i,indx) = xhkb(i)
1807        qcko(i,indx) = xqkb(i)
1808      endif
1809    enddo
1811 ! ..... static control .....
1813 ! moisture and cloud work functions
1815    do k = kts1,kmax1
1816      do i = its,ite
1817        if(cnvflg(i).and.k.gt.kb(i).and.k.le.ktcon(i)) then
1818          dz = zi(i,k+1) - zi(i,k)
1819          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1820          tem1 = 0.5 * xlamud(i) * dz
1821          factor = 1. + tem - tem1
1822          hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                           &
1823                     (heo(i,k)+heo(i,k-1)))/factor
1824        endif
1825      enddo
1826    enddo
1828    do k = kts1,kmax1
1829      do i = its,ite
1830        if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1831          alpha1 = min((-0.7*log(100.)+24.)*0.0001,c0)
1832          beta1 = 0.07
1834          if (to(i,k).gt.t0c_) then
1835            c0fac = alpha1
1836          else
1837            c0fac = alpha1*exp(beta1*(to(i,k)-t0c_))
1838          endif
1840          c0fac = max(0.0,c0fac)
1841          c0t(i,k) = c0fac
1843          dz = zi(i,k+1) - zi(i,k)
1844          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1845          xdby = hcko(i,k) - heso(i,k)
1846          xqrch = qeso(i,k)                                                     &
1847               + gamma * xdby / (hvap_ * (1. + gamma))
1848          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1849          tem1 = 0.5 * xlamud(i) * dz
1850          factor = 1. + tem - tem1
1851          qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*(qo(i,k)+qo(i,k-1)))/factor
1852          dq = eta(i,k) * qcko(i,k) - eta(i,k) * xqrch
1853          if(k.ge.kbcon(i).and.dq.gt.0.) then
1854            etah = .5 * (eta(i,k) + eta(i,k-1))
1855            if(ncloud.gt.0..and.k.gt.jmin(i)) then
1856              qlk = dq / (eta(i,k) + etah * (c0t(i,k) + c1t(i,k)) * dz)
1857            else
1858              qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1859            endif
1860            xpw = etah * c0t(i,k) * dz * qlk
1861            if(k.lt.ktcon1(i)) then
1862              xaa0(i) = xaa0(i) - (zl(i,k+1) - zl(i,k)) * g_ * qlk
1863            endif
1864            qcko(i,k) = qlk + xqrch
1865            xpwav(i) = xpwav(i) + xpw
1866          endif
1867        endif
1869        if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1870          dz1 = zl(i,k+1) - zl(i,k)
1871          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1872          rfact =  1. + fv_ * cp_ * gamma                                       &
1873                   * to(i,k) / hvap_
1874          xdby = hcko(i,k) - heso(i,k)
1875          xaa0(i) = xaa0(i)                                                     &
1876                  + dz1 * (g_ / (cp_ * to(i,k)))                                &
1877                  * xdby / (1. + gamma)                                         &
1878                  * rfact
1879          xaa0(i)=xaa0(i)+                                                      &
1880                   dz1 * g_ * fv_ *                                             &
1881                   max(0.,(qeso(i,k) - qo(i,k)))
1882        endif
1883      enddo
1884    enddo
1886 ! ..... downdraft calculations .....
1888 ! downdraft moisture properties
1890    do i = its,ite
1891      xpwev(i) = 0.
1892    enddo
1894    do i = its,ite
1895      if(cnvflg(i)) then
1896        jmn = jmin(i)
1897        xhcd(i,jmn) = heo(i,jmn)
1898        xqcd(i,jmn) = qo(i,jmn)
1899        qrcd(i,jmn) = qeso(i,jmn)
1900      endif
1901    enddo
1903    do k = kmax1,kts,-1
1904      do i = its,ite
1905        if(cnvflg(i).and.k.lt.jmin(i)) then
1906          dz = zi(i,k+2) - zi(i,k+1)
1907          if(k.ge.kbcon(i)) then
1908             tem  = xlamde * dz
1909             tem1 = 0.5 * xlamdd * dz
1910          else
1911             tem  = xlamde * dz
1912             tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1913          endif
1914          factor = 1. + tem - tem1
1915          xhcd(i,k) = ((1.-tem1)*xhcd(i,k+1)+tem*0.5*                           &
1916                     (heo(i,k)+heo(i,k+1)))/factor
1917        endif
1918      enddo
1919    enddo
1921    do k = kmax1,kts,-1
1922      do i = its,ite
1923        if(cnvflg(i).and.k.lt.jmin(i)) then
1924          dq = qeso(i,k)
1925          dt = to(i,k)
1926          gamma = el2orc * dq / dt**2
1927          dh = xhcd(i,k) - heso(i,k)
1928          qrcd(i,k) = dq+(1./hvap_)*(gamma/(1.+gamma))*dh
1929          dz = zi(i,k+2) - zi(i,k+1)
1930          if(k.ge.kbcon(i)) then
1931            tem  = xlamde * dz
1932            tem1 = 0.5 * xlamdd * dz
1933          else
1934            tem  = xlamde * dz
1935            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1936          endif
1937          factor = 1. + tem - tem1
1938          xqcd(i,k) = ((1.-tem1)*xqcd(i,k+1)+tem*0.5*                           &
1939                    (qo(i,k)+qo(i,k+1)))/factor
1940          xpwd     = etad(i,k+1) * (xqcd(i,k) - qrcd(i,k))
1941          xqcd(i,k)= qrcd(i,k)
1942          xpwev(i) = xpwev(i) + xpwd
1943        endif
1944      enddo
1945    enddo
1947    do i = its,ite
1948      edtmax = edtmaxl
1949      if(slimsk(i).eq.2.) edtmax = edtmaxs
1950      if(cnvflg(i)) then
1951        if(xpwev(i).ge.0.) then
1952          edtx(i) = 0.
1953        else
1954          edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1955          edtx(i) = min(edtx(i),edtmax)
1956        endif
1957      endif
1958    enddo
1960 ! downdraft cloudwork functions
1962    do k = kmax1,kts,-1
1963      do i = its,ite
1964        if(cnvflg(i).and.k.lt.jmin(i)) then
1965          gamma = el2orc * qeso(i,k) / to(i,k)**2
1966          dhh = xhcd(i,k)
1967          dt = to(i,k)
1968          dg = gamma
1969          dh = heso(i,k)
1970          dz =-1.*(zl(i,k+1)-zl(i,k))
1971          xaa0(i) = xaa0(i)+edtx(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))         &
1972                  *(1.+fv_*cp_*dg*dt/hvap_)
1973          xaa0(i) = xaa0(i)+edtx(i)*                                            &
1974                    dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1975        endif
1976      enddo
1977    enddo
1979 ! calculate critical cloud work function
1981    do i = its,ite
1982      tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1983      dtconv(i) = tfac * tem / wc(i)
1984      dtconv(i) = max(dtconv(i),dtmin)
1985      dtconv(i) = min(dtconv(i),dtmax)
1986 !    dtconv(i) = max(1800., dt2)
1987    enddo
1989 ! large scale forcing
1991 ! compute mean velocity between kb kbcon
1993    do i = its,ite
1994      if (cnvflg(i).and.slimsk(i).eq.2.) then
1995        wc(i) = 0.
1996        sumx(i) = 0.
1997      endif
1998    enddo
2000    do k = kts,kbmax
2001      do i = its,ite
2002        if (cnvflg(i).and.slimsk(i).eq.2.) then
2003          if (k.gt.kb(i) .and. k.le.kbcon(i)) then
2004            dz = zi(i,k+1) - zi(i,k)
2005            tem = sqrt(uo(i,k)**2 + vo(i,k)**2)
2006            wc(i) = wc(i) + tem * dz
2007            sumx(i) = sumx(i) + dz
2008          endif
2009        endif
2010      enddo
2011    enddo
2013    do i = its,ite
2014      if (cnvflg(i).and.slimsk(i).eq.2.) then
2015        wc(i) = wc(i) / sumx(i)
2016        dtpbl(i) = zl(i,kbcon(i))/wc(i)
2017      else
2018        dtpbl(i) = dtconv(i)
2019      endif
2020    enddo
2022    do i = its,ite
2023      if (cnvflg(i)) then
2024        apbl(i) = g_ * (hpbl(i)-hpbl_hold(i))/dt2 * dtpbl(i)
2025        apbl(i) = min(max(apbl(i),0.0),aa1(i))
2026      endif
2027    enddo
2029    do i = its,ite
2030      if(cnvflg(i)) then
2031        f(i) = (aa1(i) - apbl(i)) / dtconv(i)
2032        if(f(i).le.0.) cnvflg(i) = .false.
2033      endif
2034      if(cnvflg(i)) then
2035        xk(i) = (xaa0(i) - aa1(i)) / mbdt
2036        if(xk(i).ge.0.) cnvflg(i) = .false.
2037      endif
2039 ! kernel, cloud base mass flux
2041      if(cnvflg(i)) then
2042        xmb(i) = -f(i) / xk(i)
2043        xmb(i) = xmb(i) * clear(i) * (1.-sigma)
2044        xmb(i) = min(xmb(i),xmbmax(i))
2045      endif
2046      pden(i) = p(i,kbcon(i))/to(i,kbcon(i))/rd_
2047      pdot(i) = 10.* dot(i,kbcon(i))
2048 !    if (pden(i)*pdot(i).gt.xmb(i)) cnvflg(i) = .false.
2049      if (xmb(i).lt.pdot(i)/g_) cnvflg(i) = .false.
2050    enddo
2051    totflg = .true.
2052    do i = its,ite
2053      totflg = totflg .and. (.not. cnvflg(i))
2054    enddo
2055    if(totflg) return
2057 ! restore t0 and qo to t1 and q1 in case convection stops
2059    do k = kts,kmax
2060      do i = its,ite
2061        if (cnvflg(i)) then
2062          to(i,k) = t1(i,k)
2063          qo(i,k) = q1(i,k)
2064          uo(i,k) = u1(i,k)
2065          vo(i,k) = v1(i,k)
2066          qeso(i,k) = 0.01*fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2067          qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
2068          qeso(i,k) = max(qeso(i,k),qmin_)
2069        endif
2070      enddo
2071    enddo
2073 ! feedback: simply the changes from the cloud with unit mass flux
2074 !           multiplied by  the mass flux necessary to keep the
2075 !           equilibrium with the larger-scale.
2077    do i = its,ite
2078      delhbar(i) = 0.
2079      delqbar(i) = 0.
2080      deltbar(i) = 0.
2081      qcond(i) = 0.
2082      delubar(i) = 0.
2083      delvbar(i) = 0.
2084    enddo
2086    do k = kts,kmax
2087      do i = its,ite
2088        if (cnvflg(i).and.k.le.ktcon(i).and.dellaq(i,k).le.0.) then
2089          if (q1(i,k).gt.0.) then
2090            tem = dellaq(i,k) * xmb(i) * dt2
2091            dellaq(i,k) = max(tem,-q1(i,k))/(xmb(i)*dt2)
2092          else
2093            dellaq(i,k) = 0.0
2094          endif
2095        endif
2096      enddo
2097    enddo
2099    do k = kts,kmax
2100      do i = its,ite
2101        if(cnvflg(i).and.k.le.ktcon(i)) then
2102          dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
2103          t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2104          q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2105          tem=1./rcs
2106          u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
2107          v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem 
2108          dp = 1000. * del(i,k)
2109          delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
2110          delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
2111          deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
2112          delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
2113          delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
2114        endif
2115      enddo
2116    enddo
2118    do k = kts,kmax 
2119      do i = its,ite 
2120        if (cnvflg(i) .and. k.le.ktcon(i)) then
2121          qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2122          qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
2123          qeso(i,k) = max(qeso(i,k), qmin_)
2124        endif
2125      enddo
2126    enddo
2128    do i = its,ite 
2129      rntot(i) = 0.
2130      delqev(i) = 0.
2131      delq2(i) = 0.
2132      flg(i) = cnvflg(i) 
2133    enddo
2135 !  comptute rainfall  
2137    do k = kmax,kts,-1
2138      do i = its,ite
2139        if(cnvflg(i).and.k.lt.ktcon(i)) then
2140          aup = 1.
2141          if(k.le.kb(i)) aup = 0.
2142          adw = 1.
2143          if(k.ge.jmin(i)) adw = 0.
2144          rntot(i) = rntot(i)                                                   &
2145                + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
2146                * xmb(i) * .001 * dt2
2147        endif
2148      enddo
2149    enddo
2151 !  conversion rainfall (m) and compute the evaporation of falling raindrops 
2153    do k = kmax,kts,-1
2154      do i = its,ite
2155        delq(i) = 0.0
2156        deltv(i) = 0.0
2157        qevap(i) = 0.0
2158        if(cnvflg(i).and.k.lt.ktcon(i)) then
2159          aup = 1.
2160          if(k.le.kb(i)) aup = 0.
2161          adw = 1.
2162          if(k.ge.jmin(i)) adw = 0.
2163          rain(i) = rain(i)                                                     &
2164                + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
2165                * xmb(i) * .001 * dt2
2166          qrsk(i,k) = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2167        endif
2169        if(cnvflg(i).and.flg(i).and.k.lt.ktcon(i)) then
2170          ccn_f = 1.0
2171          evef = edt_s(i) * evfacts * ccn_f
2172          if(slimsk(i).eq.1.) evef = edt_s(i) * evfactl * ccn_f
2173          qcond(i) = evef * (q1(i,k) - qeso(i,k)) / (1. + el2orc *              &
2174                   qeso(i,k) / t1(i,k)**2)
2175          dp = 1000. * del(i,k)
2176          if(rain(i).gt.0..and.qcond(i).lt.0.) then
2177            qevap(i) = -qcond(i) * (1. - exp(-.32 * sqrt(dt2 * rain(i))))
2178            qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
2179            delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
2180            if (delq2(i).gt.rntot(i)) then
2181              qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
2182              flg(i) = .false.
2183            endif 
2184          endif
2185          if(rain(i).gt.0..and.qevap(i).gt.0.) then
2186            q1(i,k) = q1(i,k) + qevap(i)
2187            t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
2188            rain(i) = rain(i) - .001 * qevap(i) * dp / g_
2189            delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
2190            deltv(i) =  - (hvap_/cp_)*qevap(i)/dt2
2191            delq(i) =  + qevap(i)/dt2
2192          endif
2193          dellaq(i,k) = dellaq(i,k) + delq(i)/xmb(i)
2194          delqbar(i)  = delqbar(i) + delq(i)*dp/g_
2195          deltbar(i)  = deltbar(i) + deltv(i)*dp/g_
2196        endif
2198        if (cnvflg(i).and.k.lt.ktcon(i)) then
2199          qrs(i,k) = max(qrsk(i,k) - max(qevap(i),0.),0.)
2200          qci(i,k) = max(qcirs(i,k) - aup*pwo(i,k),0.)
2201        endif
2202      enddo
2203    enddo
2205 ! consider the negative rain in the event of rain evaporation and downdrafts
2207    do i = its,ite
2208      if(cnvflg(i)) then
2209        if(rain(i).lt.0..and..not.flg(i)) rain(i) = 0.
2210        if(rain(i).le.0.) then
2211          rain(i) = 0.
2212        else
2213          ktop(i) = ktcon(i)
2214          kbot(i) = kbcon(i)
2215          icps(i) = 1
2216        endif
2217      endif
2218    enddo
2220    do k = kts,kmax
2221      do i = its,ite
2222        if(cnvflg(i).and.rain(i).le.0.) then
2223           t1(i,k) = to(i,k)
2224           q1(i,k) = qo(i,k)
2225           u1(i,k) = uo(i,k)
2226           v1(i,k) = vo(i,k)
2227        endif
2228      enddo
2229    enddo
2231 !  detrainment of cloud water and ice
2233    if (ncloud.gt.0) then
2234      do k = kts,kmax 
2235        do i = its,ite 
2236          if (cnvflg(i) .and. rain(i).gt.0.) then
2237            if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
2238              tem  = dellal(i,k) * xmb(i) * dt2
2239              tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2240              if (ncloud.ge.2) then
2241                qi2(i,k) = qi2(i,k) + tem * tem1            ! ice
2242                qc2(i,k) = qc2(i,k) + tem *(1.0-tem1)       ! water
2243              else
2244                qc2(i,k) = qc2(i,k) + tem
2245              endif
2246            endif
2247          endif
2248        enddo
2249      enddo
2250    endif
2252    end subroutine nsas2d
2253 !-------------------------------------------------------------------------------
2254    REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2255 !-------------------------------------------------------------------------------
2256    IMPLICIT NONE
2257 !-------------------------------------------------------------------------------
2258    REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
2259            xai,xbi,ttp,tr
2260    INTEGER :: ice
2262    ttp=t0c+0.01
2263    dldt=cvap-cliq
2264    xa=-dldt/rv
2265    xb=xa+hvap/(rv*ttp)
2266    dldti=cvap-cice
2267    xai=-dldti/rv
2268    xbi=xai+hsub/(rv*ttp)
2269    tr=ttp/t
2270    if(t.lt.ttp.and.ice.eq.1) then
2271      fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2272    else
2273      fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2274    endif
2276    if (t.lt.180.) then
2277      tr=ttp/180.
2278      if(t.lt.ttp.and.ice.eq.1) then
2279        fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2280      else
2281        fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2282      endif
2283    endif
2285    if (t.ge.330.) then
2286      tr=ttp/330
2287      if(t.lt.ttp.and.ice.eq.1) then
2288        fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2289      else
2290        fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2291      endif
2292    endif
2294    END FUNCTION fpvs
2295 !-------------------------------------------------------------------------------
2297 !-------------------------------------------------------------------------------
2298    subroutine nsasinit(rthcuten,rqvcuten,rqccuten,rqicuten,                    &
2299                       rucuten,rvcuten,                                         &  
2300                       restart,p_qc,p_qi,p_first_scalar,                        &
2301                       allowed_to_read,                                         &
2302                       ids, ide, jds, jde, kds, kde,                            &
2303                       ims, ime, jms, jme, kms, kme,                            &
2304                       its, ite, jts, jte, kts, kte                  )
2305 !-------------------------------------------------------------------------------
2306    implicit none
2307 !-------------------------------------------------------------------------------
2308    logical , intent(in)           ::  allowed_to_read,restart
2309    integer , intent(in)           ::  ids, ide, jds, jde, kds, kde,            &
2310                                       ims, ime, jms, jme, kms, kme,            &
2311                                       its, ite, jts, jte, kts, kte
2312    integer , intent(in)           ::  p_first_scalar, p_qi, p_qc
2313    real,     dimension( ims:ime , kms:kme , jms:jme ) , intent(out) ::         &
2314                                                               rthcuten,        &
2315                                                               rqvcuten,        &
2316                                                                rucuten,        &
2317                                                                rvcuten,        &
2318                                                               rqccuten,        &
2319                                                               rqicuten
2320    integer :: i, j, k, itf, jtf, ktf
2322    jtf=min0(jte,jde-1)
2323    ktf=min0(kte,kde-1)
2324    itf=min0(ite,ide-1)
2326    if(.not.restart)then
2327      do j = jts,jtf
2328        do k = kts,ktf
2329          do i = its,itf
2330            rthcuten(i,k,j)=0.
2331            rqvcuten(i,k,j)=0.
2332            rucuten(i,k,j)=0.   
2333            rvcuten(i,k,j)=0.   
2334          enddo
2335        enddo
2336      enddo
2338      if (p_qc .ge. p_first_scalar) then
2339        do j = jts,jtf
2340          do k = kts,ktf
2341            do i = its,itf
2342              rqccuten(i,k,j)=0.
2343            enddo
2344          enddo
2345        enddo
2346      endif
2348      if (p_qi .ge. p_first_scalar) then
2349        do j = jts,jtf
2350          do k = kts,ktf
2351            do i = its,itf
2352              rqicuten(i,k,j)=0.
2353            enddo
2354          enddo
2355        enddo
2356      endif
2357    endif
2359    end subroutine nsasinit
2360 !-------------------------------------------------------------------------------
2362 END MODULE module_cu_ksas