integrating recent changes from fermi-test into fermi

(i'll use "fermi" for more disruptive changes to the pipeline model such
as updating the MSHRs and getting rid of the warp tracker, ripping out DWF, etc...)

[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 7805]

1 files changed, 0 insertions, 741 deletions

diff --git a/benchmarks/CUDA/WP/body_inline.h b/benchmarks/CUDA/WP/body_inline.h
deleted file mode 100644
index 5da0941..0000000
--- a/benchmarks/CUDA/WP/body_inline.h
+++ /dev/null
@@ -1,741 +0,0 @@
-#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(t)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(q)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qc)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qi)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qr)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qs)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(den)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(p)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(delz)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(cpm)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(xl)
-}
-#endif
-
-// 585 !----------------------------------------------------------------
-// 586 !     initialize the large scale variables
-     mstep = 1 ;
-
-     ttp=t0c+0.01 ;
-     dldt=cvap-cliq ;
-     xa=-dldt/rv ;
-     xb=xa+hvap/(rv*ttp) ;
-     dldti=cvap-cice ;
-     xai=-dldti/rv ;
-     xbi=xai+hsub/(rv*ttp) ;
-
-     for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-
-       pp = p[k] ;
-       tt = t[k] ;
-       tr = ttp/tt ;
-       ltr = log(tr) ;
-
-       qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
-       qq=ep2*qq/(pp-qq) ;
-       qs1[k] = MAX(qq,qmin) ;
-       rh1[k] = MAX( q[k]/qs1[k],qmin) ;
-
-       if( tt < ttp ) {
-         qq=psat*exp(ltr*(xai)+xbi*(1.-tr)) ;
-       } else {
-         qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
-       }
-       qq = ep2 * qq / (pp - qq) ;
-       qs2[k] = MAX(qq,qmin) ;
-       rh2[k] = MAX(q[k]/qs2[k],qmin) ;
-
-     }
-
-     for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          prevp[k] = 0. ;
-          psdep[k] = 0. ;
-          praut[k] = 0. ;
-          psaut[k] = 0. ;
-          pracw[k] = 0. ;
-          psaci[k] = 0. ;
-          psacw[k] = 0. ;
-          pigen[k] = 0. ;
-          pidep[k] = 0. ;
-          pcond[k] = 0. ;
-          psmlt[k] = 0. ;
-          psevp[k] = 0. ;
-          falk1[k] = 0. ;
-          falk2[k] = 0. ;
-          fall1[k] = 0. ;
-          fall2[k] = 0. ;
-          fallc[k] = 0. ;
-          falkc[k] = 0. ;
-          xni[k] = 1.e3 ;
-      }
-
-#define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
-#define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
-// calculate mstep for this colum
-
-
-     for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-       float supcol = t0c - t[k] ;
-#ifdef DEVICEEMU
-if ( ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"ZAP t0c %25.17e\n",t0c) ;
-if ( ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"ZAP supcol %25.17e\n",supcol) ;
-#endif
-DIAGOUTPUT1(t)
-       n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
-       if ( qr[k] <= qcrmin ) {
-         rsloper[k]  = rslopermax ;
-         rslopebr[k] = rsloperbmax ;
-         rslope2r[k] = rsloper2max ;
-         rslope3r[k] = rsloper3max ;
-       } else {
-DIAGOUTPUT1(qr)
-DIAGOUTPUT1(den)
-         rsloper[k]  = 1./LAMDAR(qr[k],den[k]) ;
-DIAGOUTPUT1(rsloper)
-         rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
-         rslope2r[k] = rsloper[k] * rsloper[k] ; 
-         rslope3r[k] = rslope2r[k] * rsloper[k] ; 
-       }
-       if ( qs[k] <= qcrmin ) {
-         rslopes[k]  = rslopesmax ;
-DIAGOUTPUT1(rslopes) ;
-         rslopebs[k] = rslopesbmax ;
-DIAGOUTPUT1(rslopebs) ;
-         rslope2s[k] = rslopes2max ;
-         rslope3s[k] = rslopes3max ;
-       } else {
-DIAGOUTPUT1(qs) ;
-DIAGOUTPUT1(den) ;
-DIAGOUTPUT1(n0sfac) ;
-         rslopes[k] = 1./LAMDAS(qs[k],den[k],n0sfac[k]) ;
-DIAGOUTPUT1(rslopes) ;
-         rslopebs[k] = exp(log(rslopes[k])*bvts) ;
-DIAGOUTPUT1(rslopebs) ;
-         rslope2s[k] = rslopes[k] * rslopes[k] ; 
-         rslope3s[k] = rslope2s[k] * rslopes[k] ; 
-       }
-       denfac[k] = sqrt(den0/den[k]) ;
-       w1[k] = pvtr*rslopebr[k]*denfac[k]/delz[k] ;
-       w2[k] = pvts*rslopebs[k]*denfac[k]/delz[k] ;
-
-DIAGOUTPUT1(w1)
-DIAGOUTPUT1(rslopebr)
-DIAGOUTPUT1(w2)
-DIAGOUTPUT1(rslopebs)
-DIAGOUTPUT1(denfac)
-DIAGOUTPUT1(delz)
-
-       w = MAX(w1[k],w2[k]) ;
-       numdt = MAX(trunc(w*dtcld+.5+.5),1) ;
-       if ( numdt >= mstep ) mstep = numdt ;
-//-------------------------------------------------------------
-// Ni: ice crystal number concentration   [HDC 5c]
-//-------------------------------------------------------------
-       float temp = (den[k]*MAX(qi[k],qmin)) ;
-       temp = sqrt(sqrt(temp*temp*temp)) ;
-#ifdef DEBUGDEBUG
-       xni[k] = 1.e3 ;
-#else
-       xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
-#endif
-     }
-     rmstep = 1./mstep ;
-   
-     for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-       fall1[k] = 0. ;
-       fall2[k] = 0. ;
-     }
-
-//     diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y        ! 8.794e-5*x**1.81/y
-//     viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
-//     xka(x,y) = 1.414e3*viscos(x,y)*y
-//     diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
-//     venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))    &
-//                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
-
-#define DIFFUS(x,y) (8.794e-5 * exp(log(x)*(1.81)) / (y))
-#define VISCOS(x,y) (1.496e-6 * ((x)*sqrt(x)) /((x)+120.)/(y))
-#define XKA(x,y) (1.414e3*VISCOS((x),(y))*(y))
-#define DIFFAC(a,b,c,d,e) ((d)*(a)*(a)/(XKA((c),(d))*rv*(c)*(c))+1./((e)*DIFFUS((c),(b))))
-#define VENFAC(a,b,c) (exp(log((VISCOS((b),(c))/DIFFUS((b),(a))))*((.3333333)))/sqrt(VISCOS((b),(c)))*sqrt(sqrt(den0/(c))))
-#define CONDEN(a,b,c,d,e) ((MAX((b),qmin)-(c))/(1.+(d)*(d)/(rv*(e))*(c)/((a)*(a))))
-
-     for ( n = 1 ; n <= mstep ; n++ ) {
-       k = kpe - 1 ;
-       falk1[k] = den[k]*qr[k]*w1[k]*rmstep ;
-       fall1[k] += falk1[k] ;
-       falk2[k] = den[k]*qs[k]*w2[k]*rmstep ;
-       fall2[k] += falk2[k] ;
-       qr[k] = MAX(qr[k]-falk1[k]*dtcldden,0.) ;
-       qs[k] = MAX(qs[k]-falk2[k]*dtcldden,0.) ;
-       for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
-         falk1[k] = den[k]*qr[k]*w1[k]*rmstep ;
-         fall1[k] += falk1[k] ;
-         falk2[k] = den[k]*qs[k]*w2[k]*rmstep ;
-         fall2[k] += falk2[k] ;
-         dtcldden = dtcld/den[k] ;
-         rdelz = 1./delz[k] ;
-DIAGOUTPUT1i(loop) ;
-DIAGOUTPUT1i(mstep) ;
-DIAGOUTPUT1i(n) ;
-DIAGOUTPUT1(qr) ;
-DIAGOUTPUT1(falk1) ;
-DIAGOUTPUT11(falk1) ;
-DIAGOUTPUT1(delz) ;
-DIAGOUTPUT11(delz) ;
-         qr[k] = MAX(qr[k]-
-                            (falk1[k]-falk1[k+1]*delz[k+1]*rdelz)*
-                            dtcldden,0.) ;
-DIAGOUTPUT1(qr) ;
-DIAGOUTPUT1(qs) ;
-DIAGOUTPUT1(falk2) ;
-DIAGOUTPUT1(w2) ;
-DIAGOUTPUT11(falk2) ;
-         qs[k] = MAX(qs[k]-
-                            (falk2[k]-falk2[k+1]*delz[k+1]*rdelz)*
-                            dtcldden,0.) ;
-DIAGOUTPUT1(qs) ;
-       }
-
-       for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
-DIAGOUTPUT1(t) ;
-DIAGOUTPUT1(qs) ;
-         if ( t[k] > t0c && qs[k] > 0.) {
-           xlf = xlf0 ;
-           w3[k] = VENFAC(p[k],t[k],den[k]) ;
-           coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[2]) ;
-           psmlt[k] = XKA(t[k],den[k])/xlf*(t0c-t[k])*pi/2.
-                     *n0sfac[k]*(precs1*rslope2s[k]+precs2
-                     *w3[k]*coeres) ;
-           psmlt[k] = MIN(MAX(psmlt[k]*dtcld*rmstep,-qs[k]*rmstep),0.) ;
-           qs[k] += psmlt[k] ;
-DIAGOUTPUT1i(mstep) ;
-DIAGOUTPUT1i(n) ;
-DIAGOUTPUT1(qr) ;
-DIAGOUTPUT1(psmlt) ;
-           qr[k] -= psmlt[k] ;
-DIAGOUTPUT1(qr) ;
-
-           t[k] += xlf/CPMCAL(q[k])*psmlt[k] ;
-         }
-       }
-     }
-
-//---------------------------------------------------------------
-// Vice [ms-1] : fallout of ice crystal [HDC 5a]
-//---------------------------------------------------------------
-     mstep = 1 ;
-     numdt = 1 ;
-     for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
-       if (qi[k] <= 0.) {
-         w2[k] = 0. ;
-       } else {
-         xmi = den[k]*qi[k]/xni[k] ;
-         diameter  = MAX(MIN(dicon * sqrt(xmi),dimax), 1.e-25) ;
-         w1[k] = 1.49e4*exp(log(diameter)*(1.31)) ;
-         w2[k] = w1[k]/delz[k] ;
-       }
-       numdt = MAX( trunc(w2[k]*dtcld+.5+.5),1) ;
-       if(numdt > mstep) mstep = numdt ;
-     }
-     rmstep = 1./mstep ;
-
-     for ( n = 1 ; n <= mstep ; n++ ) {
-       k = kpe - 1 ;
-       falkc[k] = den[k]*qi[k]*w2[k]*rmstep ;
-       fallc[k] = fallc[k]+falkc[k] ;
-       qi[k] = MAX(qi[k]-falkc[k]*dtcld/den[k],0.) ;
-       for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
-         falkc[k] = den[k]*qi[k]*w2[k]*rmstep ;
-         fallc[k] = fallc[k]+falkc[k] ;
-         qi[k] = MAX(qi[k]-(falkc[k]-falkc[k+1]
-                 *delz[k+1]/delz[k])*dtcld/den[k],0.) ;
-       }
-     }
-     fallsum = fall1[1]+fall2[1]+fallc[1] ;
-     fallsum_qsi = fall2[1]+fallc[1] ;
-     rainncv = 0. ;
-     if(fallsum > 0.) {
-       rainncv = fallsum*delz[1]/denr*dtcld*1000. ;
-       rain = fallsum*delz[1]/denr*dtcld*1000. + rain ;
-     }
-     snowncv = 0. ;
-     if(fallsum_qsi > 0.) {
-       snowncv = fallsum_qsi*delz[0]/denr*dtcld*1000. ;
-       snow = fallsum_qsi*delz[0]/denr*dtcld*1000. + snow ;
-     }
-     sr = 0. ;
-     if ( fallsum > 0. ) sr = fallsum_qsi*delz[0]/denr*dtcld*1000./(rainncv+1.e-12) ;
-
-//---------------------------------------------------------------
-// pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
-//       (T>T0: I->C)
-//---------------------------------------------------------------
-     for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-       float supcol = t0c-t[k] ;
-       xlf = xls-xl[k] ;
-       if( supcol < 0. ) xlf = xlf0 ;
-       if( supcol < 0 && qi[k] > 0. ) {
-         qc[k] = qc[k] + qi[k] ;
-         t[k] = t[k] - xlf/cpm[k]*qi[k] ;
-         qi[k] = 0. ;
-       }
-//---------------------------------------------------------------
-// pihmf: homogeneous freezing of cloud water below -40c [HL A45]
-//        (T<-40C: C->I)
-//---------------------------------------------------------------
-       if( supcol > 40. && qc[k] > 0. ) {
-         qi[k] = qi[k] + qc[k] ;
-         t[k] = t[k] + xlf/cpm[k]*qc[k] ;
-         qc[k] = 0. ;
-       }
-//---------------------------------------------------------------
-// pihtf: heterogeneous freezing of cloud water [HL A44]
-//        (T0>T>-40C: C->I)
-//---------------------------------------------------------------
-       if ( supcol > 0. && qc[k] > 0.) {
-         float pfrzdtc = MIN(pfrz1*(exp(pfrz2*supcol)-1.)
-           *den[k]/denr/xncr*qc[k]*qc[k]*dtcld,qc[k]) ;
-         qi[k] = qi[k] + pfrzdtc ;
-         t[k] = t[k] + xlf/cpm[k]*pfrzdtc ;
-         qc[k] = qc[k]-pfrzdtc ;
-       }
-//---------------------------------------------------------------
-// psfrz: freezing of rain water [HL A20] [LFO 45]
-//        (T<T0, R->S)
-//---------------------------------------------------------------
-       if( supcol > 0. && qr[k] > 0. ) {
-         float temp = rsloper[k] ;
-         temp = temp*temp*temp*temp*temp*temp*temp ;
-         float pfrzdtr = MIN(20.*(pi*pi)*pfrz1*n0r*denr/den[k]
-               *(exp(pfrz2*supcol)-1.)*temp*dtcld,
-               qr[k]) ;
-         qs[k] = qs[k] + pfrzdtr ;
-         t[k] = t[k] + xlf/cpm[k]*pfrzdtr ;
-         qr[k] = qr[k]-pfrzdtr ;
-DIAGOUTPUT1(qr) ;
-       }
-     }
-
-//----------------------------------------------------------------
-//     rsloper: reverse of the slope parameter of the rain(m)
-//     xka:    thermal conductivity of air(jm-1s-1k-1)
-//     work1:  the thermodynamic term in the denominator associated with
-//             heat conduction and vapor diffusion
-//             (ry88, y93, h85)
-//     work2: parameter associated with the ventilation effects(y93)
-
-     for ( k = kps-1 ; k <= kpe - 1 ; k++ ) {
-       float supcol = t0c - t[k] ;
-       n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
-       if ( qr[k] <= qcrmin ) {
-         rsloper[k]  = rslopermax ;
-         rslopebr[k] = rsloperbmax ;
-         rslope2r[k] = rsloper2max ;
-         rslope3r[k] = rsloper3max ;
-       } else {
-         rsloper[k] = 1./(sqrt(sqrt(pidn0r/((qr[k])*(den[k]))))) ;
-DIAGOUTPUT1(rsloper) ;
-DIAGOUTPUT1(qr) ;
-DIAGOUTPUT1(den) ;
-         rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
-         rslope2r[k] = rsloper[k] * rsloper[k] ;
-         rslope3r[k] = rslope2r[k] * rsloper[k] ;
-       }
-       if ( qs[k] <= qcrmin ) {
-         rslopes[k]  = rslopesmax ;
-         rslopebs[k] = rslopesbmax ;
-         rslope2s[k] = rslopes2max ;
-         rslope3s[k] = rslopes3max ;
-       } else {
-         rslopes[k] = 1./(sqrt(sqrt(pidn0s*(n0sfac[k])/((qs[k])*(den[k]))))) ;
-         rslopebs[k] = exp(log(rslopes[k])*bvts) ;
-         rslope2s[k] = rslopes[k] * rslopes[k] ;
-         rslope3s[k] = rslope2s[k] * rslopes[k] ;
-       }
-    }
-
-    for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-       w1[k] = DIFFAC(xl[k],p[k],t[k],den[k],qs1[k]) ;
-       w2[k] = DIFFAC(xls,p[k],t[k],den[k],qs2[k]) ;
-       w3[k] = VENFAC(p[k],t[k],den[k]) ;
-    }
-
-//
-//===============================================================
-//
-// warm rain processes
-//
-// - follows the processes in RH83 and LFO except for autoconcersion
-//
-//===============================================================
-//
-    for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-      float supsat = MAX(q[k],qmin)-qs1[k] ;
-      float satdt = supsat/dtcld ;
-//---------------------------------------------------------------
-// praut: auto conversion rate from cloud to rain [HDC 16]
-//        (C->R)
-//---------------------------------------------------------------
-      if(qc[k] > qc0) {
-        praut[k] = qck1*exp(log(qc[k])*((7./3.))) ;
-        praut[k] = MIN(praut[k],qc[k]/dtcld) ;
-      }
-//---------------------------------------------------------------
-// pracw: accretion of cloud water by rain [HL A40] [LFO 51]
-//        (C->R)
-//---------------------------------------------------------------
-      if(qr[k] > qcrmin && qc[k] > qmin) {
-        pracw[k] = MIN(pacrr*rslope3r[k]*rslopebr[k]
-                   *qc[k]*denfac[k],qc[k]/dtcld) ;
-      }
-//---------------------------------------------------------------
-// prevp: evaporation/condensation rate of rain [HDC 14]
-//        (V->R or R->V)
-//---------------------------------------------------------------
-      if(qr[k] > 0.) {
-        coeres = rslope2r[k]*sqrt(rsloper[k]*rslopebr[k]) ;
-        prevp[k] = (rh1[k]-1.)*(precr1*rslope2r[k]
-                     +precr2*w3[k]*coeres)/w1[k] ;
-DIAGOUTPUT1(prevp) ;
-DIAGOUTPUT1(qr) ;
-DIAGOUTPUT1(rsloper) ;
-DIAGOUTPUT1(rslope2r) ;
-DIAGOUTPUT1(rslopebr) ;
-DIAGOUTPUT1(w1) ;
-DIAGOUTPUT1(rh1) ;
-        if(prevp[k] < 0.) {
-          prevp[k] = MAX(prevp[k],-qr[k]/dtcld) ;
-          prevp[k] = MAX(prevp[k],satdt/2) ;
-        } else {
-          prevp[k] = MIN(prevp[k],satdt/2) ;
-        }
-      }
-    }
-
-//
-//===============================================================
-//
-// cold rain processes
-//
-// - follows the revised ice microphysics processes in HDC
-// - the processes same as in RH83 and RH84  and LFO behave
-//   following ice crystal hapits defined in HDC, inclduing
-//   intercept parameter for snow (n0s), ice crystal number
-//   concentration (ni), ice nuclei number concentration
-//   (n0i), ice diameter (d)
-//
-//===============================================================
-//
-    rdtcld = 1./dtcld ;
-    for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          float supcol = t0c-t[k] ;
-          float supsat = MAX(q[k],qmin)-qs2[k] ;
-          float satdt = supsat/dtcld ;
-          int ifsat = 0 ;
-//-------------------------------------------------------------
-// Ni: ice crystal number concentraiton   [HDC 5c]
-//-------------------------------------------------------------
-          float temp = (den[k]*MAX(qi[k],qmin)) ;
-          temp = sqrt(sqrt(temp*temp*temp)) ;
-          xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
-          float eacrs = exp(0.07*(-supcol)) ;
-//
-          if(supcol > 0) {
-            if(qs[k] > qcrmin && qi[k] > qmin) {
-              xmi = den[k]*qi[k]/xni[k] ;
-              diameter  = MIN(dicon * sqrt(xmi),dimax) ;
-              vt2i = 1.49e4*pow(diameter,(float)1.31) ;
-              vt2s = pvts*rslopebs[k]*denfac[k] ;
-//-------------------------------------------------------------
-// psaci: Accretion of cloud ice by rain [HDC 10]
-//        (T<T0: I->S)
-//-------------------------------------------------------------
-              acrfac = 2.*rslope3s[k]+2.*diameter*rslope2s[k]
-                      +diameter*diameter*rslopes[k] ;
-              psaci[k] = pi*qi[k]*eacrs*n0s*n0sfac[k]
-                           *abs(vt2s-vt2i)*acrfac*.25 ;
-            }
-//-------------------------------------------------------------
-// psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
-//        (T<T0: C->S, and T>=T0: C->R)
-//-------------------------------------------------------------
-            if(qs[k] > qcrmin && qc[k] > qmin) {
-              psacw[k] = MIN(pacrc*n0sfac[k]*rslope3s[k] 
-                           *rslopebs[k]*qc[k]*denfac[k]
-                           ,qc[k]*rdtcld) ;
-            }
-//-------------------------------------------------------------
-// pidep: Deposition/Sublimation rate of ice [HDC 9]
-//       (T<T0: V->I or I->V)
-//-------------------------------------------------------------
-            if(qi[k] > 0 && ifsat != 1) {
-              xmi = den[k]*qi[k]/xni[k] ;
-              diameter = dicon * sqrt(xmi) ;
-              pidep[k] = 4.*diameter*xni[k]*(rh2[k]-1.)/w2[k] ;
-              supice = satdt-prevp[k] ;
-              if(pidep[k] < 0.) {
-                pidep[k] = MAX(MAX(pidep[k],satdt*.5),supice) ;
-                pidep[k] = MAX(pidep[k],-qi[k]*rdtcld) ;
-              } else {
-                pidep[k] = MIN(MIN(pidep[k],satdt*.5),supice) ;
-              }
-              if(abs(prevp[k]+pidep[k]) >= abs(satdt)) ifsat = 1 ;
-            }
-          }
-//-------------------------------------------------------------
-// psdep: deposition/sublimation rate of snow [HDC 14]
-//        (V->S or S->V)
-//-------------------------------------------------------------
-          if( qs[k] > 0. && ifsat != 1) {
-            coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[k]) ;
-            psdep[k] = (rh2[k]-1.)*n0sfac[k]
-                         *(precs1*rslope2s[k]+precs2
-                         *w3[k]*coeres)/w2[k] ;
-            supice = satdt-prevp[k]-pidep[k] ;
-            if(psdep[k] < 0.) {
-              psdep[k] = MAX(psdep[k],-qs[k]*rdtcld) ;
-              psdep[k] = MAX(MAX(psdep[k],satdt*.5),supice) ;
-            } else {
-              psdep[k] = MIN(MIN(psdep[k],satdt*.5),supice) ;
-            }
-            if(abs(prevp[k]+pidep[k]+psdep[k]) >= abs(satdt))
-              ifsat = 1 ;
-          }
-//-------------------------------------------------------------
-// pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
-//       (T<T0: V->I)
-//-------------------------------------------------------------
-          if(supcol > 0) {
-            if(supsat > 0 && ifsat != 1) {
-              supice = satdt-prevp[k]-pidep[k]-psdep[k] ; 
-              xni0 = 1.e3*exp(0.1*supcol) ;
-              roqi0 = 4.92e-11*exp(log(xni0)*(1.33));
-              pigen[k] = MAX(0.,(roqi0/den[k]-MAX(qi[k],0.))
-                         *rdtcld) ;
-              pigen[k] = MIN(MIN(pigen[k],satdt),supice) ;
-            }
-//
-//-------------------------------------------------------------
-// psaut: conversion(aggregation) of ice to snow [HDC 12]
-//       (T<T0: I->S)
-//-------------------------------------------------------------
-            if(qi[k] > 0.) {
-              qimax = roqimax/den[k] ;
-              psaut[k] = MAX(0.,(qi[k]-qimax)*rdtcld) ;
-            }
-          }
-//-------------------------------------------------------------
-// psevp: Evaporation of melting snow [HL A35] [RH83 A27]
-//       (T>T0: S->V)
-//-------------------------------------------------------------
-          if(supcol < 0.) {
-            if(qs[k] > 0. && rh1[k] < 1.) {
-              psevp[k] = psdep[k]*w2[k]/w1[k] ;
-            }  // asked Jimy about this, 11.6.07, JM
-            psevp[k] = MIN(MAX(psevp[k],-qs[k]*rdtcld),0.) ;
-          }
-      }
-
-
-//
-//
-//----------------------------------------------------------------
-//     check mass conservation of generation terms and feedback to the
-//     large scale
-//
-      for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          if(t[k]<=t0c) {
-//
-//     cloud water
-//
-            value = MAX(qmin,qc[k]) ;
-            source = (praut[k]+pracw[k]+psacw[k])*dtcld ;
-            if (source > value) {
-              factor = value/source ;
-              praut[k] = praut[k]*factor ;
-              pracw[k] = pracw[k]*factor ;
-              psacw[k] = psacw[k]*factor ;
-            }
-//
-//     cloud ice
-//
-            value = MAX(qmin,qi[k]) ;
-            source = (psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld ;
-            if (source > value) {
-              factor = value/source ;
-              psaut[k] = psaut[k]*factor ;
-              psaci[k] = psaci[k]*factor ;
-              pigen[k] = pigen[k]*factor ;
-              pidep[k] = pidep[k]*factor ;
-            }
-//
-            w3[k]=-(prevp[k]+psdep[k]+pigen[k]+pidep[k]) ;
-//     update
-DIAGOUTPUT1(q) ;
-DIAGOUTPUT1(prevp) ;
-DIAGOUTPUT1(psdep) ;
-DIAGOUTPUT1(pigen) ;
-DIAGOUTPUT1(pidep) ;
-            q[k] = q[k]+w3[k]*dtcld ;
-DIAGOUTPUT1(q) ;
-            qc[k] = MAX(qc[k]-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
-            qr[k] = MAX(qr[k]+(praut[k]+pracw[k]+prevp[k])*dtcld,0.) ;
-            qi[k] = MAX(qi[k]-(psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld,0.) ;
-DIAGOUTPUT1(qs)
-            qs[k] = MAX(qs[k]+(psdep[k]+psaut[k]+psaci[k]+psacw[k])*dtcld,0.) ;
-DIAGOUTPUT1(qs)
-            xlf = xls-xl[k] ;
-            xlwork2 = -xls*(psdep[k]+pidep[k]+pigen[k])-xl[k]*prevp[k]-xlf*psacw[k] ;
-            t[k] = t[k]-xlwork2/cpm[k]*dtcld ;
-          } else {
-//
-//     cloud water
-//
-            value = MAX(qmin,qc[k]) ;
-            source=(praut[k]+pracw[k]+psacw[k])*dtcld ;
-            if (source > value) {
-              factor = value/source ;
-              praut[k] = praut[k]*factor ;
-              pracw[k] = pracw[k]*factor ;
-              psacw[k] = psacw[k]*factor ;
-            }
-//
-//     snow
-//
-            value = MAX(qcrmin,qs[k]) ;
-            source=(-psevp[k])*dtcld ;
-            if (source > value) {
-              factor = value/source ;
-              psevp[k] = psevp[k]*factor ;
-            }
-            w3[k]=-(prevp[k]+psevp[k]) ;
-//     update
-DIAGOUTPUT1(q) ;
-DIAGOUTPUT1(prevp) ;
-DIAGOUTPUT1(psdep) ;
-DIAGOUTPUT1(pigen) ;
-DIAGOUTPUT1(pidep) ;
-            q[k] = q[k]+w3[k]*dtcld ;
-DIAGOUTPUT1(q) ;
-            qc[k] = MAX(qc[k]-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
-            qr[k] = MAX(qr[k]+(praut[k]+pracw[k]+prevp[k] +psacw[k])*dtcld,0.) ;
-DIAGOUTPUT1(qs)
-DIAGOUTPUT1(psevp)
-
-#ifdef DEVICEEMU
-if (ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"%8s %25.17e\n","ZAP p*dt",psevp[k]*dtcld) ;
-if (ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"%8s %25.17e\n","ZAP q+p*dt",qs[k]+psevp[k]*dtcld) ;
-#endif
-            qs[k] = MAX(qs[k]+psevp[k]*dtcld,0.) ;
-DIAGOUTPUT1(qs)
-            xlf = xls-xl[k] ;
-            xlwork2 = -xl[k]*(prevp[k]+psevp[k]) ;
-            t[k] = t[k]-xlwork2/cpm[k]*dtcld ;
-          }
-      }
-DIAGOUTPUT2(qs)
-//
-// Inline expansion for fpvs
-      cvap = cpv ;
-      ttp=t0c+0.01 ;
-      dldt=cvap-cliq ;
-      xa=-dldt/rv ;
-      xb=xa+hvap/(rv*ttp) ;
-      dldti=cvap-cice ;
-      xai=-dldti/rv ;
-      xbi=xai+hsub/(rv*ttp) ;
-      for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          tr=ttp/t[k] ;
-          qs1[k]=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) ;
-          qs1[k] = ep2 * qs1[k] / (p[k] - qs1[k]) ;
-          qs1[k] = MAX(qs1[k],qmin) ;
-      }
-//
-//----------------------------------------------------------------
-//  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
-//     if there exists additional water vapor condensated/if
-//     evaporation of cloud water is not enough to remove subsaturation
-//
-      for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          w1[k] = ((MAX(q[k],qmin)-(qs1[k])) /
-            (1.+(xl[k])*(xl[k])/(rv*(cpm[k]))*(qs1[k])/((t[k])*(t[k])))) ;
-          // w3[k] = qc[k]+w1[k] ;   NOT USED
-          pcond[k] = MIN(MAX(w1[k]/dtcld,0.),MAX(q[k],0.)/dtcld) ;
-          if(qc[k] > 0. && w1[k] < 0.) {
-            pcond[k] = MAX(w1[k],-qc[k])/dtcld ;
-          }
-DIAGOUTPUT1(q) ;
-DIAGOUTPUT1(pcond) ;
-DIAGOUTPUT1(qs1) ;
-          q[k] = q[k]-pcond[k]*dtcld ;
-DIAGOUTPUT1(q) ;
-          qc[k] = MAX(qc[k]+pcond[k]*dtcld,0.) ;
-          t[k] = t[k]+pcond[k]*xl[k]/cpm[k]*dtcld ;
-      }
-//
-//
-//----------------------------------------------------------------
-//     padding for small values
-//
-      for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-          if(qc[k] <= qmin) qc[k] = 0.0 ;
-          if(qi[k] <= qmin) qi[k] = 0.0 ;
-      }
-
-//////////// end of loop ////////////////
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(t)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(q)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qc)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qi)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qr)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(qs)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(den)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(p)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(delz)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(cpm)
-}
-for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
-kDIAGOUTPUT1(xl)
-}