

      subroutine solve(nstep,nloop2,nrec,prec,nwrite,nrst,           &
                   taptim,rsttim,cloudvar,qname,budname,qbudget,asq,bsq,   &
                   xh,rxh,uh,ruh,xf,rxf,uf,ruf,yh,vh,rvh,yf,vf,rvf,  &
                   sigma,sigmaf,tauh,taus,zh,mh,rmh,tauf,zf,mf,rmf,  &
                   pi0,rho0,prs0,thv0,th0,qv0,qc0,                   &
                   ql0,rr0,rf0,rrf0,                                 &
                   zs,gz,dzdx,dzdy,rain,sws,thflux,qvflux,rth0s,     &
                   radbcw,radbce,radbcs,radbcn,                      &
                   dum1,dum2,dum3,dum4,divx,rho,prs,                 &
                   t11,t12,t13,t22,t23,t33,                          &
                   gx,u0,rru,ua,u3d,uten,uten1,                      &
                   gy,v0,rrv,va,v3d,vten,vten1,                      &
                   rrw,wa,w3d,wten,wten1,ppi,pp3d,ppten,sten,        &
                   tha,th3d,thten,thten1,thterm,tk,qa,q3d,qten,      &
                   kmh,kmv,khh,khv,tkea,tke3d,tketen,pta,pt3d,ptten, &
                   pdata,cfb,cfa,cfc,d1,d2,pdt,deft,rhs,trans,       &
                   reqs_u,reqs_v,reqs_w,reqs_s,reqs_p,reqs_tk,reqs_q,reqs_t, &
                   ww1,ww2,we1,we2,ws1,ws2,wn1,wn2,                  &
                   pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2,                  &
                   uw31,uw32,ue31,ue32,us31,us32,un31,un32,          &
                   vw31,vw32,ve31,ve32,vs31,vs32,vn31,vn32,          &
                   ww31,ww32,we31,we32,ws31,ws32,wn31,wn32,          &
                   sw31,sw32,se31,se32,ss31,ss32,sn31,sn32,          &
                   pw31,pw32,pe31,pe32,ps31,ps32,pn31,pn32,          &
                   tkw1,tkw2,tke1,tke2,tks1,tks2,tkn1,tkn2,          &
                   qw1,qw2,qe1,qe2,qs1,qs2,qn1,qn2,                  &
                   tw1,tw2,te1,te2,ts1,ts2,tn1,tn2)
      use module_mp_thompson
      use module_mp_graupel
      implicit none

!-----------------------------------------------------------------------
!  Bryan Cloud Model (CM1) release 12  (cm1r12)
!  12 August 2008
!    Bug fix:   5 January 2012  (for stretched grids with terrain)
!  http://www.mmm.ucar.edu/people/bryan/cm1/
!-----------------------------------------------------------------------

      include 'input.incl'
      include 'constants.incl'
      include 'timestat.incl'
#ifdef MPI
      include 'mpif.h'
#endif

!-----------------------------------------------------------------------
! Arrays and variables passed into solve

      integer :: nstep,nloop2,nrec,prec,nwrite,nrst
      real*8 :: taptim,rsttim
      logical, dimension(maxq) :: cloudvar
      character*3, dimension(maxq) :: qname
      character*6, dimension(maxq) :: budname
      real*8, dimension(nbudget) :: qbudget
      real*8, dimension(numq) :: asq,bsq
      real, dimension(ib:ie) :: xh,rxh,uh,ruh
      real, dimension(ib:ie+1) :: xf,rxf,uf,ruf
      real, dimension(jb:je) :: yh,vh,rvh
      real, dimension(jb:je+1) :: yf,vf,rvf
      real, dimension(kb:ke) :: sigma
      real, dimension(kb:ke+1) :: sigmaf
      real, dimension(ib:ie,jb:je,kb:ke) :: tauh,taus,zh,mh,rmh
      real, dimension(ib:ie,jb:je,kb:ke+1) :: tauf,zf,mf,rmf
      real, dimension(ib:ie,jb:je,kb:ke) :: pi0,rho0,prs0,thv0,th0,qv0,qc0
      real, dimension(ib:ie,jb:je,kb:ke) :: ql0,rr0,rf0,rrf0
      real, dimension(itb:ite,jtb:jte) :: zs,gz,dzdx,dzdy
      real, dimension(itb:ite+1,jtb:jte,ktb:kte) :: gx
      real, dimension(itb:ite,jtb:jte+1,ktb:kte) :: gy
      real, dimension(ib:ie,jb:je,nrain) :: rain,sws
      real, dimension(ib:ie,jb:je) :: thflux,qvflux,rth0s
      real, dimension(jb:je,kb:ke) :: radbcw,radbce
      real, dimension(ib:ie,kb:ke) :: radbcs,radbcn
      real, dimension(ib:ie,jb:je,kb:ke) :: dum1,dum2,dum3,dum4,divx,rho,prs
      real, dimension(ib:ie,jb:je,kb:ke) :: t11,t12,t13,t22,t23,t33
      real, dimension(ib:ie+1,jb:je,kb:ke) :: u0,rru,ua,u3d,uten,uten1
      real, dimension(ib:ie,jb:je+1,kb:ke) :: v0,rrv,va,v3d,vten,vten1
      real, dimension(ib:ie,jb:je,kb:ke+1) :: rrw,wa,w3d,wten,wten1
      real, dimension(ib:ie,jb:je,kb:ke) :: ppi,pp3d,ppten,sten
      real, dimension(ib:ie,jb:je,kb:ke) :: tha,th3d,thten,thten1,thterm,tk
      real, dimension(ibm:iem,jbm:jem,kbm:kem,numq) :: qa,q3d,qten
      real, dimension(ibc:iec,jbc:jec,kbc:kec) :: kmh,kmv,khh,khv
      real, dimension(ibt:iet,jbt:jet,kbt:ket) :: tkea,tke3d,tketen
      real, dimension(ibp:iep,jbp:jep,kbp:kep,npt) :: pta,pt3d,ptten
      real, dimension(npvals,nparcels) :: pdata
      real, dimension(ipb:ipe,jpb:jpe,kpb:kpe) :: cfb
      real, dimension(kpb:kpe) :: cfa,cfc,d1,d2
      complex, dimension(ipb:ipe,jpb:jpe,kpb:kpe) :: pdt,deft
      complex, dimension(ipb:ipe,jpb:jpe) :: rhs,trans
      integer, dimension(rmp) :: reqs_u,reqs_v,reqs_w,reqs_s,reqs_p,reqs_tk
      integer, dimension(rmp,numq) :: reqs_q
      integer, dimension(rmp,npt) :: reqs_t
      real, dimension(jmp,kmp-1) :: ww1,ww2,we1,we2
      real, dimension(imp,kmp-1) :: ws1,ws2,wn1,wn2
      real, dimension(jmp,kmp) :: pw1,pw2,pe1,pe2
      real, dimension(imp,kmp) :: ps1,ps2,pn1,pn2
      real, dimension(cmp,jmp,kmp)   :: uw31,uw32,ue31,ue32
      real, dimension(imp+1,cmp,kmp) :: us31,us32,un31,un32
      real, dimension(cmp,jmp+1,kmp) :: vw31,vw32,ve31,ve32
      real, dimension(imp,cmp,kmp)   :: vs31,vs32,vn31,vn32
      real, dimension(cmp,jmp,kmp-1) :: ww31,ww32,we31,we32
      real, dimension(imp,cmp,kmp-1) :: ws31,ws32,wn31,wn32
      real, dimension(cmp,jmp,kmp)   :: sw31,sw32,se31,se32
      real, dimension(imp,cmp,kmp)   :: ss31,ss32,sn31,sn32
      real, dimension(cmp,jmp,kmp)   :: pw31,pw32,pe31,pe32
      real, dimension(imp,cmp,kmp)   :: ps31,ps32,pn31,pn32
      real, dimension(cmp,jmp,kmt)   :: tkw1,tkw2,tke1,tke2
      real, dimension(imp,cmp,kmt)   :: tks1,tks2,tkn1,tkn2
      real, dimension(cmp,jmp,kmp,numq) :: qw1,qw2,qe1,qe2
      real, dimension(imp,cmp,kmp,numq) :: qs1,qs2,qn1,qn2
      real, dimension(cmp,jmp,kmp,npt) :: tw1,tw2,te1,te2
      real, dimension(imp,cmp,kmp,npt) :: ts1,ts2,tn1,tn2

!-----------------------------------------------------------------------
! Arrays and variables defined inside solve

      integer i,j,k,n,nrk,icom

      logical simple_comm

      real dttmp,rtime,rdt,tem,trad,cvm,ql
      real*8 afoo,bfoo


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cc   First, get turbulence coefficients   ccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

#ifdef MPI
      nf=0
      nu=0
      nv=0
      nw=0
#endif

      afoo=0.0d0
      bfoo=0.0d0


!--------------------------------------------------------------------
!  surface scheme:

      ! get surface drag
      if(idrag.eq.1)then
        call sfcdrag(u0,v0,rf0,dum1,dum2,t13,t23,ua,va)
      endif

      ! get surface flux
      if(isfcflx.eq.1)then
        call sfcflux(xh,ruh,rvh,zh,pi0,thv0,th0,u0,v0,thflux,qvflux,    &
                     rho,dum1,dum2,ua,va,ppi,tha,qa(ib,jb,kb,nqv),qbudget(8))
      endif

!--------------------------------------------------------------------
!                 For turbulence section only:
!  dum1 = Brunt-Vaisala frequency (N_m^2) (nm)
!  dum2 = Vertical deformation terms (S_v^2) (defsq)
!  dum3 = Horizontal deformation terms (S_h^2) (defh)
!
!  Arrays available for temporary storage:
!  dum4,divx,ppten,sten,thterm

      IF(iturb.ge.1)THEN

        call calcnm(mf,pi0,thv0,th0,cloudvar,dum1,dum2,dum3,dum4,divx,   &
                    prs,ppi,tha,qa)

        call calcdef(xh,rxh,uh,xf,rxf,uf,vh,vf,mh,mf,dum2,dum3,   &
                     divx,ppten,ua,va,wa,t11,t12,t13,t22,t23,t33,gx,gy)

      ENDIF

!--------------------------------------------------------------------
!  tke scheme

      IF(iturb.eq.1)THEN
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nkt
        do j=1,nj
        do i=1,ni
          tketen(i,j,k)=0.0
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_misc=time_misc+mytime()

        call turbtke(ruh,rvh,rmh,mf,rmf,th0,thflux,qvflux,rth0s,   &
                     dum1,dum2,dum3,dum4,divx,ppten,sten,thterm,   &
                     kmh,kmv,khh,khv,tkea,tketen,t13,t23,ua,va)

        call turbt(xh,rxh,uh,xf,uf,vh,vf,mh,mf,rho0,rr0,rf0,rrf0,   &
                   dum1,dum2,dum3,dum4,divx,tkea,tketen,kmh,kmv,gx,gy)

        if(idiff.eq.1)then
          if(difforder.eq.2)then
            call diff2w(1,0,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,tkea,tketen)
          elseif(difforder.eq.6)then
            call diff6w(dum1,dum2,dum3,tkea,tketen)
          endif
        endif

!-------------------------------------------------
!  Smagorinsky scheme

      ELSEIF(iturb.eq.2)THEN

        call turbsmag(ruh,rvh,rmh,mf,rmf,th0,thflux,qvflux,rth0s,  &
                      dum1,dum2,dum3,dum4,divx,sten,               &
                      kmh,kmv,khh,khv,t13,t23,ua,va)

!-------------------------------------------------
!  Smagorinsky scheme for axisymmetric simulations

      ELSEIF(iturb.eq.3)THEN

        call turbsmagaxi(ruh,rvh,rmh,mf,rmf,th0,thflux,qvflux,rth0s,  &
                         dum1,dum2,dum3,kmh,kmv,khh,khv,t13,t23,ua,va)

      ENDIF

!-------------------------------------------------
!  Get turbulent stresses:

      IF(iturb.ge.1)THEN

        call gettau(xf,rxf,rho0,rf0,kmh,kmv,t11,t12,t13,t22,t23,t33,ua)

      ENDIF

!--------------------------------------------------------------------
!  Dissipative heating term:  NOTE: epsilon is stored in the thterm array

      IF(idiss.eq.1)THEN
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          thterm(i,j,k) = 0.0
        enddo
        enddo
        enddo
      ENDIF

      IF(iturb.ge.1.and.idiss.eq.1)THEN

        call getepsilon(rxh,uh,xf,rxf,uf,yh,vh,yf,vf,mh,mf,rr0,rrf0,   &
                        dum1,dum2,dum3,dum4,divx,ppten,sten,thterm,    &
                        t11,t12,t13,t22,t23,t33,ua,va,wa,gx,gy)

      ENDIF

!--------------------------------------------------------------------


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CC   Pre-RK calculations   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

 
!-------------------------------------------------------------------
!  Parcel update

      if(iprcl.eq.1)then
        ! rtime valid at beginning of time step
        rtime=float(nstep-1)*dtl
        call parcel_driver(prec,xh,uh,ruh,yh,vh,rvh,zh,mh,rmh,mf,        &
                           pi0,thv0,th0,dum1,dum2,dum3,dum4,divx,prs,    &
                           ua,va,wa,ppi,ppten,tha,qa,khv,pdata,rtime)
      endif

!--------------------------------------------------------------------
!  radbc
 
      if(irbc.eq.1)then

        if(ibw.eq.1 .or. ibe.eq.1) call radbcew(radbcw,radbce,ua)
 
        if(ibs.eq.1 .or. ibn.eq.1) call radbcns(radbcs,radbcn,va)

      endif

!--------------------------------------------------------------------
!  set ppten to zero

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj
      do i=1,ni
        ppten(i,j,k)=0.0
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()

!--------------------------------------------------------------------
!  U-equation
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj
      do i=1,ni+1
        uten1(i,j,k)=0.
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()
 
      if(iturb.ge.1)then
        call turbu(xh,xf,rxf,uf,vh,mh,rho0,dum1,dum2,dum3,dum4,uten1,t11,t12,t13,t22,gx,gy)
      endif

      if(idiff.ge.1)then
        if(difforder.eq.2)then
          call diff2u(1,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,ua,uten1)
        elseif(difforder.eq.6)then
          call diff6u(u0,dum1,dum2,dum3,ua,uten1)
        endif
      endif

      if(irdamp.eq.1)then
        call rdampu(tauh,u0,ua,uten1)
      endif

      if(dns.eq.1)then
        call diff2u(2,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,ua,uten1)
      endif

!--------------------------------------------------------------------
!  V-equation
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj+1
      do i=1,ni
        vten1(i,j,k)=0.
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()
 
      if(iturb.ge.1)then
        call turbv(xh,rxh,uh,xf,vf,mh,rho0,dum1,dum2,dum3,dum4,vten1,t12,t22,t23,gx,gy)
      endif

      if(idiff.ge.1)then
        if(difforder.eq.2)then
          call diff2v(1,xh,uh,rxf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,va,vten1)
        elseif(difforder.eq.6)then
          call diff6v(v0,dum1,dum2,dum3,va,vten1)
        endif
      endif

      if(irdamp.eq.1)then
        call rdampv(tauh,v0,va,vten1)
      endif

      if(dns.eq.1)then
        call diff2v(2,xh,uh,rxf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,va,vten1)
      endif
 
!--------------------------------------------------------------------
!  W-equation
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk+1
      do j=1,nj
      do i=1,ni
        wten1(i,j,k)=0.0
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()
 
      if(iturb.ge.1)then
        call turbw(xh,rxh,uh,xf,vh,mf,rrf0,dum1,dum2,dum3,dum4,wten1,t13,t23,t33,gx,gy)
      endif

      if(idiff.ge.1)then
        if(difforder.eq.2)then
          call diff2w(1,1,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,wa,wten1)
        elseif(difforder.eq.6)then
          call diff6w(dum1,dum2,dum3,wa,wten1)
        endif
      endif

      if(irdamp.eq.1)then
        call rdampw(tauf,wa,wten1)
      endif

      if(dns.eq.1)then
        call diff2w(2,1,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,dum4,thterm,wa,wten1)
      endif

!--------------------------------------------------------------------
!  THETA-equation
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj
      do i=1,ni
        thten1(i,j,k)=0.0
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()

      if( (iturb.ge.1).or.(dns.eq.1) )then
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=0,nk+1
        do j=0,nj+1
        do i=0,ni+1
          sten(i,j,k)=th0(i,j,k)+tha(i,j,k)
        enddo
        enddo
        enddo
      endif

      if(iturb.ge.1)then
        call turbs(1,xh,rxh,uh,xf,uf,vh,vf,mh,mf,rho0,rr0,rf0,thflux,   &
                   dum1,dum2,dum3,dum4,divx,sten,thten1,khh,khv,gx,gy)
      endif

      if(idiff.eq.1)then
        if(difforder.eq.2)then
          call diff2s(1,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,tha,thten1)
        elseif(difforder.eq.6)then
          call diff6s(ql0,dum1,dum2,dum3,tha,thten1)
        endif
      endif

      if(irdamp.eq.1)then
        call rdamps(taus,th0,tha,thten1)
      endif

      if(dns.eq.1)then
        call diff2s(2,rxh,uh,xf,uf,vh,vf,mh,mf,dum1,dum2,dum3,sten,thten1)
      endif

      IF(idiss.eq.1)THEN
        ! use thterm array to store epsilon
        if(imoist.eq.1.and.neweqts.ge.1)then
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            dum1(i,j,k)=qa(i,j,k,nqv)
            dum2(i,j,k)=0.0
            dum3(i,j,k)=0.0
          enddo
          enddo
          enddo
          do n=nql1,nql2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              dum2(i,j,k)=dum2(i,j,k)+qa(i,j,k,n)
            enddo
            enddo
            enddo
          enddo
          if(iice.eq.1)then
          do n=nqs1,nqs2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              dum3(i,j,k)=dum3(i,j,k)+qa(i,j,k,n)
            enddo
            enddo
            enddo
          enddo
          endif
!$omp parallel do default(shared)  &
!$omp private(i,j,k,cvm)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            cvm=cv+cvv*dum1(i,j,k)+cpl*dum2(i,j,k)+cpi*dum3(i,j,k)
            thten1(i,j,k)=thten1(i,j,k)   &
                        +thterm(i,j,k)/( cpdcv*cvm*(pi0(i,j,k)+ppi(i,j,k)) )
            ppten(i,j,k)=ppten(i,j,k)+thterm(i,j,k)*rovcp/( cvm*(th0(i,j,k)+tha(i,j,k)) )
          enddo
          enddo
          enddo
        elseif(imoist.eq.1.and.neweqts.eq.0)then
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            thten1(i,j,k)=thten1(i,j,k)   &
                        +thterm(i,j,k)/( cp*(pi0(i,j,k)+ppi(i,j,k)) )
          enddo
          enddo
          enddo
        else
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            thten1(i,j,k)=thten1(i,j,k)   &
                        +thterm(i,j,k)/( cp*(pi0(i,j,k)+ppi(i,j,k)) )
            ppten(i,j,k)=ppten(i,j,k)+thterm(i,j,k)*rovcp/( cv*(th0(i,j,k)+tha(i,j,k)) )
          enddo
          enddo
          enddo
        endif
      ENDIF

!---- Rotunno-Emanuel "radiation" term
!---- (currently capped at 2 K/day ... see RE87 p 546)

      IF(rterm.eq.1)THEN
!$omp parallel do default(shared)  &
!$omp private(i,j,k,trad)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          trad = -tha(i,j,k)/(12.0*3600.0)
          if( tha(i,j,k).gt. 1.0 ) trad = -1.0/(12.0*3600.0)
          if( tha(i,j,k).lt.-1.0 ) trad =  1.0/(12.0*3600.0)
          thten1(i,j,k)=thten1(i,j,k)+trad
          ppten(i,j,k)=ppten(i,j,k)+(pi0(i,j,k)+ppi(i,j,k))*(     &
                         rddcv*trad/(th0(i,j,k)+tha(i,j,k)) )
        enddo
        enddo
        enddo
      ENDIF
      if(timestats.ge.1) time_misc=time_misc+mytime()

!-------------------------------------------------------------------
!  Passive Tracers

      if(iptra.eq.1)then
        do n=1,npt
          call gettenq(0,xh,rxh,uh,xf,uf,vh,vf,mh,mf,ql0,rho0,rr0,rf0,qvflux,   &
                       dum1,dum2,dum3,dum4,divx,                           &
                       pta(ib,jb,kb,n),ptten(ib,jb,kb,n),khh,khv,gx,gy)
        enddo
      endif

!-------------------------------------------------------------------
!  Moisture

      if(imoist.eq.1)then
        DO n=1,numq
          if(n.eq.nqv)then
            call gettenq(1,xh,rxh,uh,xf,uf,vh,vf,mh,mf,qv0,rho0,rr0,rf0,qvflux,   &
                         dum1,dum2,dum3,dum4,divx,                      &
                         qa(ib,jb,kb,n),qten(ib,jb,kb,n),khh,khv,gx,gy)
          else
            call gettenq(0,xh,rxh,uh,xf,uf,vh,vf,mh,mf,ql0,rho0,rr0,rf0,qvflux,   &
                         dum1,dum2,dum3,dum4,divx,                      &
                         qa(ib,jb,kb,n),qten(ib,jb,kb,n),khh,khv,gx,gy)
          endif
        ENDDO
      endif

!-------------------------------------------------------------------


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CC   Begin RK section   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

!--------------------------------------------------------------------
! RK2 begin

      DO NRK=1,3

        dttmp=dtl/float(4-nrk)

        IF(.not.terrain_flag)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=2,nk
          do j=0,nj+1
          do i=0,ni+1
            rrw(i,j,k)=rf0(i,j,k)*w3d(i,j,k)
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_advs=time_advs+mytime()

        ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=2,nk
          do j=1,nj
          do i=1,ni
            rrw(i,j,k)=rf0(i,j,k)*( w3d(i,j,k)*gz(i,j)                 &
                 +0.25*(sigmaf(k)-zt)/(zt-zs(i,j))*(                   &
                       ( (u3d(i+1,j,k-1)+u3d(i  ,j,k-1))               &
                        +(u3d(i+1,j,k  )+u3d(i  ,j,k  )) )*dzdx(i,j)   &
                      +( (v3d(i,j+1,k-1)+v3d(i,j  ,k-1))               &
                        +(v3d(i,j+1,k  )+v3d(i,j  ,k  )) )*dzdy(i,j) ) )
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_advs=time_advs+mytime()

          call bcw(rrw,0)
#ifdef MPI
          call comm_1w_start(rrw,ww1,ww2,we1,we2,   &
                                 ws1,ws2,wn1,wn2,reqs_w)
#endif

        ENDIF

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+1
        do i=0,ni+2
          rru(i,j,k)=0.5*(rho0(i-1,j,k)+rho0(i,j,k))*u3d(i,j,k)
        enddo
        enddo
        enddo

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+2
        do i=0,ni+1
          rrv(i,j,k)=0.5*(rho0(i,j-1,k)+rho0(i,j,k))*v3d(i,j,k)
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_advs=time_advs+mytime()

#ifdef MPI
        if(terrain_flag)then
          call comm_1w_end(rrw,ww1,ww2,we1,we2,   &
                               ws1,ws2,wn1,wn2,reqs_w)
        endif
#endif

        IF(.not.terrain_flag)THEN

!------------
        IF(axisymm.eq.0)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=0,ni+1
            divx(i,j,k)=(rru(i+1,j,k)-rru(i,j,k))*rdx*uh(i)        &
                       +(rrv(i,j+1,k)-rrv(i,j,k))*rdy*vh(j)        &
                       +(rrw(i,j,k+1)-rrw(i,j,k))*rdz*mh(i,j,k)
          enddo
          enddo
          enddo

        ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=0,ni+1
            divx(i,j,k)=(xf(i+1)*rru(i+1,j,k)-xf(i)*rru(i,j,k))*rdx*uh(i)*rxh(i)   &
                       +(rrw(i,j,k+1)-rrw(i,j,k))*rdz*mh(i,j,k)
          enddo
          enddo
          enddo

        ENDIF
!------------

        ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=0,ni+1
            divx(i,j,k)=(rru(i+1,j,k)-rru(i,j,k))*rdx*uh(i)        &
                       +(rrv(i,j+1,k)-rrv(i,j,k))*rdy*vh(j)        &
                       +(rrw(i,j,k+1)-rrw(i,j,k))*rdz*mh(i,j,k)/gz(i,j)
          enddo
          enddo
          enddo

        ENDIF
        if(timestats.ge.1) time_divx=time_divx+mytime()

!--------------------------------------------------------------------
!  TKE advection
 
        if(iturb.eq.1)then
          call integtkerk(xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,mf,           &
                          rho0,rr0,rf0,rrf0,                              &
                          dum1,dum2,dum3,dum4,divx,t12,rru,rrv,rrw,       &
                          wten,tkea,tke3d,tketen,dttmp,nrk                &
#ifdef MPI
              ,tkw1,tkw2,tke1,tke2,tks1,tks2,tkn1,tkn2,reqs_tk            &
#endif
                          )
        endif

!--------------------------------------------------------------------
!  Passive Tracers

        if(iptra.eq.1)then
          DO n=1,npt
            call integqrk(2,bfoo,xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,mf,        &
                          rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,t12,     &
                          rru,rrv,rrw,sten,pta(ib,jb,kb,n),pt3d(ib,jb,kb,n),  &
                          ptten(ib,jb,kb,n),dttmp,nrk                         &
#ifdef MPI
                     ,tw1(1,1,1,n),tw2(1,1,1,n),te1(1,1,1,n),te2(1,1,1,n)     &
                     ,ts1(1,1,1,n),ts2(1,1,1,n),tn1(1,1,1,n),tn2(1,1,1,n)     &
                     ,reqs_t(1,n)                                             &
#endif
                          )
          ENDDO
        endif

!--------------------------------------------------------------------
!  Calculate misc. variables
!
!    These arrays store variables that are used later in the
!    SOUND subroutine.  Do not modify t11 or t22 until after sound!
!    dum1 = vapor
!    dum2 = all liquid
!    dum3 = all solid

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+1
        do i=0,ni+1
          dum1(i,j,k)=0.0
          dum2(i,j,k)=0.0
          dum3(i,j,k)=0.0
        enddo
        enddo
        enddo

        IF(imoist.eq.1)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=0,ni+1
            dum1(i,j,k)=max(q3d(i,j,k,nqv),0.0)
          enddo
          enddo
          enddo

          DO n=nql1,nql2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=0,nj+1
            do i=0,ni+1
              dum2(i,j,k)=dum2(i,j,k)+q3d(i,j,k,n)
            enddo
            enddo
            enddo
          ENDDO

          IF(iice.eq.1)THEN
            DO n=nqs1,nqs2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
              do k=1,nk
              do j=0,nj+1
              do i=0,ni+1
                dum3(i,j,k)=dum3(i,j,k)+q3d(i,j,k,n)
              enddo
              enddo
              enddo
            ENDDO
          ENDIF

        ENDIF

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+1
        do i=0,ni+1
          t11(i,j,k)=(th0(i,j,k)+th3d(i,j,k))*(1.0+reps*dum1(i,j,k))     &
                     /(1.0+dum1(i,j,k)+max(0.0,dum2(i,j,k))+max(0.0,dum3(i,j,k)))
        enddo
        enddo
        enddo

        IF(thsmall.eq.1.and.imoist.eq.1)THEN
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            t12(i,j,k)=dum1(i,j,k)
            t13(i,j,k)=max(0.0,dum2(i,j,k))+max(0.0,dum3(i,j,k))
          enddo
          enddo
          enddo
        ENDIF

        if(timestats.ge.1) time_buoyan=time_buoyan+mytime()

!--------------------------------------------------------------------
!  THETA-equation
 
        icom=1
        if(imoist.eq.1.and.neweqts.eq.2) icom=0

        call integthrk(icom,xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,mf,          &
                       th0,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,          &
                       divx,t23,t22,gx,rru,u3d,gy,rrv,v3d,rrw,w3d,         &
                       tha,th3d,thten1,thten,thterm,tk,q3d,dttmp,nrk       &
#ifdef MPI
         ,sw31,sw32,se31,se32,ss31,ss32,sn31,sn32,reqs_s                   &
#endif
                       )

!--------------------------------------------------------------------
! Moisture

        if(imoist.eq.1)then
          DO n=1,numq
            icom=1
            if(neweqts.eq.2.and.(n.eq.nqv.or.n.eq.2.or.n.eq.4)) icom=0
            call integqrk(icom,bsq(n),xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,mf,   &
                          rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,t23,     &
                          rru,rrv,rrw,sten,qa(ib,jb,kb,n),q3d(ib,jb,kb,n),    &
                          qten(ib,jb,kb,n),dttmp,nrk                          &
#ifdef MPI
                     ,qw1(1,1,1,n),qw2(1,1,1,n),qe1(1,1,1,n),qe2(1,1,1,n)     &
                     ,qs1(1,1,1,n),qs2(1,1,1,n),qn1(1,1,1,n),qn2(1,1,1,n)     &
                     ,reqs_q(1,n)                                             &
#endif
                          )
          ENDDO
        endif

!--------------------------------------------------------------------
!  Pressure equation

      IF(psolver.le.3)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          sten(i,j,k)=ppten(i,j,k)
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_misc=time_misc+mytime()

        if(hadvorder.eq.5)then
          call adv5s(0,bfoo,xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,         &
                     rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,t23,   &
                     rru,rrv,rrw,ppi,pp3d,sten,0,dttmp)
        elseif(hadvorder.eq.6)then
          call adv6s(0,bfoo,xh,rxh,uh,ruh,xf,vh,rvh,gz,mh,rmh,         &
                     rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,t23,   &
                     rru,rrv,rrw,ppi,pp3d,sten,0,dttmp)
        endif

        IF(imoist.eq.1.and.neweqts.eq.2)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            t33(i,j,k)=ppi(i,j,k)+dttmp*sten(i,j,k)
          enddo
          enddo
          enddo

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            t23(i,j,k)=t33(i,j,k)
          enddo
          enddo
          enddo

          IF(thsmall.eq.1)THEN
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              th3d(i,j,k)=tha(i,j,k)+dttmp*thten(i,j,k)
            enddo
            enddo
            enddo
          ENDIF
          if(timestats.ge.1) time_satadj=time_satadj+mytime()

          call calcprs(pi0,prs,t33)
          call calcrho(pi0,th0,rho,prs,t33,th3d,q3d)
          if(ptype.eq.1)then
            call satadj(nrk,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,   &
                        rho,dum3,t33,prs,th3d,                      &
                        q3d(ib,jb,kb,1),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3))
          elseif(ptype.eq.2)then
            call satadj_ice(nrk,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,   &
                            rho,dum3,t33,prs,th3d,                      &
                q3d(ib,jb,kb,1),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3),   &
                q3d(ib,jb,kb,4),q3d(ib,jb,kb,5),q3d(ib,jb,kb,6))
          elseif(ptype.eq.6)then
            call satadj(nrk,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,   &
                        rho,dum3,t33,prs,th3d,                      &
                        q3d(ib,jb,kb,1),q3d(ib,jb,kb,2),ql0)
          endif

          IF(thsmall.eq.1)THEN
            rdt=1.0/dttmp
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              thten(i,j,k)=(th3d(i,j,k)-tha(i,j,k))*rdt
            enddo
            enddo
            enddo
            if(timestats.ge.1) time_satadj=time_satadj+mytime()
          ENDIF

          if(nrk.lt.3)then
            if(thsmall.eq.0)then
              call bcs(th3d)
#ifdef MPI
              call comm_3s_start(th3d,sw31,sw32,se31,se32,   &
                                      ss31,ss32,sn31,sn32,reqs_s)
#endif
            endif
            call bcs(q3d(ib,jb,kb,nqv))
#ifdef MPI
            call comm_3s_start(q3d(ib,jb,kb,nqv),                            &
              qw1(1,1,1,nqv),qw2(1,1,1,nqv),qe1(1,1,1,nqv),qe2(1,1,1,nqv),   &
              qs1(1,1,1,nqv),qs2(1,1,1,nqv),qn1(1,1,1,nqv),qn2(1,1,1,nqv),   &
              reqs_q(1,nqv))
#endif
            call bcs(q3d(ib,jb,kb,2))
#ifdef MPI
            call comm_3s_start(q3d(ib,jb,kb,2),                            &
              qw1(1,1,1,2),qw2(1,1,1,2),qe1(1,1,1,2),qe2(1,1,1,2),   &
              qs1(1,1,1,2),qs2(1,1,1,2),qn1(1,1,1,2),qn2(1,1,1,2),   &
              reqs_q(1,2))
#endif
          endif
          if(nrk.lt.3.and.iice.eq.1)then
            call bcs(q3d(ib,jb,kb,4))
#ifdef MPI
            call comm_3s_start(q3d(ib,jb,kb,4),                            &
              qw1(1,1,1,4),qw2(1,1,1,4),qe1(1,1,1,4),qe2(1,1,1,4),   &
              qs1(1,1,1,4),qs2(1,1,1,4),qn1(1,1,1,4),qn2(1,1,1,4),   &
              reqs_q(1,4))
#endif
          endif

          rdt=1.0/dttmp

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            sten(i,j,k)=sten(i,j,k)+(t33(i,j,k)-t23(i,j,k))*rdt
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_satadj=time_satadj+mytime()

        ENDIF

      ENDIF

!--------------------------------------------------------------------
!  U-equation

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni+1
          uten(i,j,k)=uten1(i,j,k)
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_misc=time_misc+mytime()
 
        if(icor.eq.1)then

        if(pertcor.eq.1)then

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj+1
          do i=0,ni+1
            dum1(i,j,k)=v3d(i,j,k)-v0(i,j,k)
          enddo
          enddo
          enddo

        else

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj+1
          do i=0,ni+1
            dum1(i,j,k)=v3d(i,j,k)
          enddo
          enddo
          enddo

        endif

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni+1
            uten(i,j,k)=uten(i,j,k)+fcor*             &
             0.25*( (dum1(i  ,j,k)+dum1(i  ,j+1,k))   &
                   +(dum1(i-1,j,k)+dum1(i-1,j+1,k)) )
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_cor=time_cor+mytime()

        endif

        if(axisymm.eq.1)then

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=2,ni+1
            uten(i,j,k)=uten(i,j,k)+0.5*(   &
                 ( v3d(i-1,j,k)**2)*rxh(i-1)+(v3d(i,j,k)**2)*rxh(i) )
          enddo
          enddo
          enddo

        endif

        if(hadvorder.eq.5)then
          call adv5u(xf,rxf,uf,vh,gz,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,u3d,uten,rrv,rrw)
        elseif(hadvorder.eq.6)then
          call adv6u(xf,rxf,uf,vh,gz,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,u3d,uten,rrv,rrw)
        endif

!--------------------------------------------------------------------
!  V-equation
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=1,nj+1
        do i=1,ni
          vten(i,j,k)=vten1(i,j,k)
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_misc=time_misc+mytime()
 
        if(icor.eq.1)then

!--------------
          IF(axisymm.eq.0)THEN

        if(pertcor.eq.1)then

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=1,ni+1
            dum1(i,j,k)=u3d(i,j,k)-u0(i,j,k)
          enddo
          enddo
          enddo

        else

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=0,nj+1
          do i=1,ni+1
            dum1(i,j,k)=u3d(i,j,k)
          enddo
          enddo
          enddo

        endif

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj+1
          do i=1,ni
            vten(i,j,k)=vten(i,j,k)-fcor*             &
             0.25*( (dum1(i,j  ,k)+dum1(i+1,j  ,k))   &
                   +(dum1(i,j-1,k)+dum1(i+1,j-1,k)) )
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_cor=time_cor+mytime()

          ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            vten(i,j,k)=vten(i,j,k)-fcor*0.5*(xf(i)*u3d(i,j,k)+xf(i+1)*u3d(i+1,j,k))*rxh(i)
          enddo
          enddo
          enddo

          ENDIF
!--------------

        endif

        if(axisymm.eq.1)then

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            vten(i,j,k)=vten(i,j,k)-(v3d(i,j,k)*rxh(i))*0.5*(xf(i)*u3d(i,j,k)+xf(i+1)*u3d(i+1,j,k))*rxh(i)
          enddo
          enddo
          enddo

        endif

        if(hadvorder.eq.5)then
          call adv5v(xh,rxh,uh,xf,vf,gz,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,rrv,v3d,vten,rrw)
        elseif(hadvorder.eq.6)then
          call adv6v(xh,rxh,uh,xf,vf,gz,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,rrv,v3d,vten,rrw)
        endif

!--------------------------------------------------------------------
!  W-equation

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk+1
        do j=1,nj
        do i=1,ni
          wten(i,j,k)=wten1(i,j,k)
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_misc=time_misc+mytime()
 
        if( (thsmall.eq.0) .or. (thsmall.eq.1.and.imoist.eq.1) )   &
        call buoyan(thv0,qv0,qc0,dum2,wten,t11,t12,t13)

        if(hadvorder.eq.5)then
          call adv5w(xh,rxh,uh,xf,vh,gz,mf,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,rrv,rrw,w3d,wten)
        elseif(hadvorder.eq.6)then
          call adv6w(xh,rxh,uh,xf,vh,gz,mf,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,divx,  &
                     rru,rrv,rrw,w3d,wten)
        endif


#ifdef MPI
!--------------------------------------------------------------------
!  Finish comm

        if(iturb.eq.1)then
          call comm_3t_end(tke3d,tkw1,tkw2,tke1,tke2,   &
                                 tks1,tks2,tkn1,tkn2,reqs_tk)
        endif
        if(iptra.eq.1)then
          do n=1,npt
            call comm_3s_end(pt3d(ib,jb,kb,n),                           &
                  tw1(1,1,1,n),tw2(1,1,1,n),te1(1,1,1,n),te2(1,1,1,n),   &
                  ts1(1,1,1,n),ts2(1,1,1,n),tn1(1,1,1,n),tn2(1,1,1,n),   &
                  reqs_t(1,n))
          enddo
        endif
        IF(nrk.lt.3.or.imoist.eq.0)THEN
          icom=1
          if(imoist.eq.1.and.neweqts.eq.2) icom=0
          if(icom.eq.1)then
            call comm_3s_end(th3d,sw31,sw32,se31,se32,   &
                                  ss31,ss32,sn31,sn32,reqs_s)
          endif
        ENDIF
        IF(nrk.lt.3.and.imoist.eq.1)THEN
          DO n=1,numq
            icom=1
            if(neweqts.eq.2.and.(n.eq.nqv.or.n.eq.2.or.n.eq.4)) icom=0
            if(icom.eq.1)then
              call comm_3s_end(q3d(ib,jb,kb,n),                            &
                    qw1(1,1,1,n),qw2(1,1,1,n),qe1(1,1,1,n),qe2(1,1,1,n),   &
                    qs1(1,1,1,n),qs2(1,1,1,n),qn1(1,1,1,n),qn2(1,1,1,n),   &
                    reqs_q(1,n))
            endif
          ENDDO
        ENDIF
        icom=1
        if(imoist.eq.1.and.neweqts.eq.2) icom=0
        IF(nrk.lt.3)THEN
          icom=1
          if(imoist.eq.1.and.neweqts.eq.2) icom=0
          if(icom.eq.0.and.thsmall.eq.0)then
            call comm_3s_end(th3d,sw31,sw32,se31,se32,   &
                                  ss31,ss32,sn31,sn32,reqs_s)
          endif
          if(imoist.eq.1)then
            DO n=1,numq
              icom=1
              if(neweqts.eq.2.and.(n.eq.nqv.or.n.eq.2.or.n.eq.4)) icom=0
              if(icom.eq.0)then
                call comm_3s_end(q3d(ib,jb,kb,n),                            &
                      qw1(1,1,1,n),qw2(1,1,1,n),qe1(1,1,1,n),qe2(1,1,1,n),   &
                      qs1(1,1,1,n),qs2(1,1,1,n),qn1(1,1,1,n),qn2(1,1,1,n),   &
                      reqs_q(1,n))
              endif
            ENDDO
          endif
        ENDIF
        icom=1
        if(imoist.eq.1.and.neweqts.eq.2) icom=0

#endif
!--------------------------------------------------------------------

        IF(axisymm.eq.1)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            v3d(i,j,k)=va(i,j,k)+dttmp*vten(i,j,k)
          enddo
          enddo
          enddo
          if(timestats.ge.1) time_misc=time_misc+mytime()

          call bcv(v3d)

        ENDIF

!--------------------------------------------------------------------
!  call sound

        IF(psolver.eq.1)THEN

          call soundns(xh,rxh,uh,xf,uf,vh,vf,mh,mf,pi0,thv0,       &
                       radbcw,radbce,radbcs,radbcn,                &
                       divx,ua,u3d,uten,va,v3d,vten,wa,w3d,wten,   &
                       ppi,pp3d,sten,t11,t22,dttmp,nrk,            &
                       reqs_u,reqs_v,reqs_w,reqs_s,reqs_p,         &
                       uw31,uw32,ue31,ue32,us31,us32,un31,un32,    &
                       vw31,vw32,ve31,ve32,vs31,vs32,vn31,vn32,    &
                       ww31,ww32,we31,we32,ws31,ws32,wn31,wn32,    &
                       sw31,sw32,se31,se32,ss31,ss32,sn31,sn32,    &
                       pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2)

        ELSEIF(psolver.eq.2)THEN

          call sounde(xh,rxh,uh,ruh,xf,uf,vh,rvh,vf,mh,rmh,mf,rmf,   &
                     pi0,thv0,rho0,rr0,rf0,th0,dzdx,dzdy,            &
                     radbcw,radbce,radbcs,radbcn,                    &
                     dum1,dum2,dum3,dum4,divx,                       &
                     gx,ua,u3d,uten,gy,va,v3d,vten,wa,w3d,wten,      &
                     ppi,pp3d,sten,tha,th3d,thten,thterm,tk,         &
                     t11,t22,nrk,sigma,sigmaf,                       &
                     reqs_u,reqs_v,reqs_w,reqs_s,reqs_p,             &
                     uw31,uw32,ue31,ue32,us31,us32,un31,un32,        &
                     vw31,vw32,ve31,ve32,vs31,vs32,vn31,vn32,        &
                     ww31,ww32,we31,we32,ws31,ws32,wn31,wn32,        &
                     sw31,sw32,se31,se32,ss31,ss32,sn31,sn32,        &
                     pw31,pw32,pe31,pe32,ps31,ps32,pn31,pn32,        &
                     pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2)

        ELSEIF(psolver.eq.3)THEN

          call sound(xh,rxh,uh,ruh,xf,uf,vh,rvh,vf,mh,rmh,mf,rmf,    &
                     pi0,thv0,rho0,rr0,rf0,th0,dzdx,dzdy,            &
                     radbcw,radbce,radbcs,radbcn,                    &
                     dum1,dum2,dum3,dum4,t12,t13,t23,t33,            &
                     gx,ua,u3d,uten,gy,va,v3d,vten,wa,w3d,wten,      &
                     ppi,pp3d,sten,tha,th3d,thten,thterm,tk,         &
                     t11,t22,nrk,sigma,sigmaf,                       &
                     reqs_u,reqs_v,reqs_w,reqs_s,reqs_p,             &
                     uw31,uw32,ue31,ue32,us31,us32,un31,un32,        &
                     vw31,vw32,ve31,ve32,vs31,vs32,vn31,vn32,        &
                     ww31,ww32,we31,we32,ws31,ws32,wn31,wn32,        &
                     sw31,sw32,se31,se32,ss31,ss32,sn31,sn32,        &
                     pw31,pw32,pe31,pe32,ps31,ps32,pn31,pn32,        &
                     pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2)

        ELSEIF(psolver.eq.4.or.psolver.eq.5)THEN

          call anelp(uh,uf,vh,vf,mh,rmh,mf,rmf,pi0,thv0,rho0,prs0,rf0,   &
                     radbcw,radbce,radbcs,radbcn,dum1,divx,              &
                     ua,u3d,uten,va,v3d,vten,wa,w3d,wten,                &
                     ppi,pp3d,t11,cfb,cfa,cfc,d1,d2,pdt,deft,rhs,trans,dttmp,nrk)

        ENDIF

!--------------------------------------------------------------------
!  radbc

        if(irbc.eq.4)then

          if(ibw.eq.1 .or. ibe.eq.1)then
            call radbcew4(ruf,radbcw,radbce,ua,u3d,dttmp)
          endif

          if(ibs.eq.1 .or. ibn.eq.1)then
            call radbcns4(rvf,radbcs,radbcn,va,v3d,dttmp)
          endif

        endif

!--------------------------------------------------------------------
! RK loop end

      ENDDO


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CC   End of RK section   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


!--------------------------------------------------------------------
!  Get pressure

    IF(psolver.eq.4.or.psolver.eq.5)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj
      do i=1,ni
        prs(i,j,k)=prs0(i,j,k)
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()

    ELSE

      call calcprs(pi0,prs,pp3d)

    ENDIF

!--------------------------------------------------------------------
!  Get density

    IF(psolver.eq.4.or.psolver.eq.5)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=1,nk
      do j=1,nj
      do i=1,ni
        rho(i,j,k)=rho0(i,j,k)
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_misc=time_misc+mytime()

    ELSE

      call calcrho(pi0,th0,rho,prs,pp3d,th3d,q3d)

    ENDIF

!-----------------------------------------------------------------
!  Explicit microphysics

      IF(imoist.eq.1)THEN

!-----------------------------------------------------------------------
!  store t in dum1

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          dum1(i,j,k)=(th0(i,j,k)+th3d(i,j,k))*(pi0(i,j,k)+pp3d(i,j,k))
        enddo
        enddo
        enddo

!-----------------------------------------------------------------------
!  prep for efall calculation:  store cvm in dum2

        IF(efall.eq.1)THEN
          if(neweqts.ge.1)then
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              dum2(i,j,k)=qa(i,j,k,nqv)
              dum3(i,j,k)=0.0
              dum4(i,j,k)=0.0
            enddo
            enddo
            enddo
            do n=nql1,nql2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
              do k=1,nk
              do j=1,nj
              do i=1,ni
                dum3(i,j,k)=dum3(i,j,k)+qa(i,j,k,n)
              enddo
              enddo
              enddo
            enddo
            if(iice.eq.1)then
              do n=nqs1,nqs2
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
                do k=1,nk
                do j=1,nj
                do i=1,ni
                  dum4(i,j,k)=dum4(i,j,k)+qa(i,j,k,n)
                enddo
                enddo
                enddo
              enddo
            endif
          else
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
            do k=1,nk
            do j=1,nj
            do i=1,ni
              dum2(i,j,k)=0.0
              dum3(i,j,k)=0.0
              dum4(i,j,k)=0.0
            enddo
            enddo
            enddo
          endif
!$omp parallel do default(shared)  &
!$omp private(i,j,k,n,ql)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            dum2(i,j,k)=cv+cvv*dum2(i,j,k)+cpl*dum3(i,j,k)+cpi*dum4(i,j,k)
          enddo
          enddo
          enddo
        ENDIF

!-----------------------------------------------------------------------
!  store rr in dum3

        IF(axisymm.eq.0)THEN

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            dum3(i,j,k)=rho(i,j,k)
          enddo
          enddo
          enddo

        ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            dum3(i,j,k)=rho(i,j,k)*2.0*pi*xh(i)
          enddo
          enddo
          enddo

        ENDIF

!-----------------------------------------------------------------------
!  Kessler scheme
        IF(ptype.eq.1)THEN
          simple_comm = .false.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-14,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
          call pdefq(1.0e-14,asq(3),ruh,rvh,rmh,rho,q3d(ib,jb,kb,3))
          call k_fallout(rho,q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          call geterain(cpl,lv1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          if(efall.ge.1)then
            call getefall(cpl,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          call kessler(qbudget(3),qbudget(4),qbudget(5),ruh,rvh,rmh,pi0,th0,dum1,   &
                       rho,dum3,pp3d,th3d,prs,                            &
                       q3d(ib,jb,kb,nqv),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3))
          call bcs(q3d(ib,jb,kb,3))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,3),                               &
                     qw1(1,1,1,3),qw2(1,1,1,3),qe1(1,1,1,3),qe2(1,1,1,3),   &
                     qs1(1,1,1,3),qs2(1,1,1,3),qn1(1,1,1,3),qn2(1,1,1,3),   &
                     reqs_q(1,3))
#endif
          call satadj(4,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,   &
                      rho,dum3,pp3d,prs,th3d,                   &
                      q3d(ib,jb,kb,nqv),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3))
          call bcs(q3d(ib,jb,kb,1))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,1),                               &
                     qw1(1,1,1,1),qw2(1,1,1,1),qe1(1,1,1,1),qe2(1,1,1,1),   &
                     qs1(1,1,1,1),qs2(1,1,1,1),qn1(1,1,1,1),qn2(1,1,1,1),   &
                     reqs_q(1,1))
#endif
          call bcs(q3d(ib,jb,kb,2))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,2),                               &
                     qw1(1,1,1,2),qw2(1,1,1,2),qe1(1,1,1,2),qe2(1,1,1,2),   &
                     qs1(1,1,1,2),qs2(1,1,1,2),qn1(1,1,1,2),qn2(1,1,1,2),   &
                     reqs_q(1,2))
#endif
          call bcs(th3d)
#ifdef MPI
          call comm_3s_start(th3d,sw31,sw32,se31,se32,   &
                                  ss31,ss32,sn31,sn32,reqs_s)
#endif
          call bcs(pp3d)
#ifdef MPI
          call comm_3s_start(pp3d,pw31,pw32,pe31,pe32,   &
                                  ps31,ps32,pn31,pn32,reqs_p)
#endif
!-----------------------------------------------------------------------
!  Goddard scheme
        ELSEIF(ptype.eq.2)THEN
          simple_comm = .false.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-14,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
          call pdefq(1.0e-14,asq(3),ruh,rvh,rmh,rho,q3d(ib,jb,kb,3))
          call pdefq(1.0e-14,asq(4),ruh,rvh,rmh,rho,q3d(ib,jb,kb,4))
          call pdefq(1.0e-14,asq(5),ruh,rvh,rmh,rho,q3d(ib,jb,kb,5))
          call pdefq(1.0e-14,asq(6),ruh,rvh,rmh,rho,q3d(ib,jb,kb,6))
          call goddard(qbudget(3),qbudget(4),qbudget(5),ruh,rvh,rmh,pi0,th0,             &
                       rho,dum3,prs,pp3d,th3d,                            &
     q3d(ib,jb,kb,1), q3d(ib,jb,kb,2),q3d(ib,jb,kb,3),qten(ib,jb,kb,3),   &
     q3d(ib,jb,kb,4),qten(ib,jb,kb,4),q3d(ib,jb,kb,5),qten(ib,jb,kb,5),   &
     q3d(ib,jb,kb,6),qten(ib,jb,kb,6))
          call satadj_ice(4,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,     &
                          rho,dum3,pp3d,prs,th3d,                     &
              q3d(ib,jb,kb,1),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3),   &
              q3d(ib,jb,kb,4),q3d(ib,jb,kb,5),q3d(ib,jb,kb,6))
          call bcs(q3d(ib,jb,kb,1))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,1),                               &
                     qw1(1,1,1,1),qw2(1,1,1,1),qe1(1,1,1,1),qe2(1,1,1,1),   &
                     qs1(1,1,1,1),qs2(1,1,1,1),qn1(1,1,1,1),qn2(1,1,1,1),   &
                     reqs_q(1,1))
#endif
          call bcs(q3d(ib,jb,kb,2))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,2),                               &
                     qw1(1,1,1,2),qw2(1,1,1,2),qe1(1,1,1,2),qe2(1,1,1,2),   &
                     qs1(1,1,1,2),qs2(1,1,1,2),qn1(1,1,1,2),qn2(1,1,1,2),   &
                     reqs_q(1,2))
#endif
          call geterain(cpl,lv1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          call geterain(cpi,ls1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,4),qten(ib,jb,kb,4))
          call geterain(cpi,ls1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,5),qten(ib,jb,kb,5))
          call geterain(cpi,ls1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,6),qten(ib,jb,kb,6))
          if(efall.ge.1)then
            call getefall(cpl,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
            call getefall(cpi,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,4),qten(ib,jb,kb,4))
            call getefall(cpi,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,5),qten(ib,jb,kb,5))
            call getefall(cpi,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,6),qten(ib,jb,kb,6))
          endif
          call bcs(th3d)
#ifdef MPI
          call comm_3s_start(th3d,sw31,sw32,se31,se32,   &
                                  ss31,ss32,sn31,sn32,reqs_s)
#endif
          call bcs(pp3d)
#ifdef MPI
          call comm_3s_start(pp3d,pw31,pw32,pe31,pe32,   &
                                  ps31,ps32,pn31,pn32,reqs_p)
#endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,3),qten(ib,jb,kb,3))
          call bcs(q3d(ib,jb,kb,3))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,3),                               &
                     qw1(1,1,1,3),qw2(1,1,1,3),qe1(1,1,1,3),qe2(1,1,1,3),   &
                     qs1(1,1,1,3),qs2(1,1,1,3),qn1(1,1,1,3),qn2(1,1,1,3),   &
                     reqs_q(1,3))
#endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,4),qten(ib,jb,kb,4))
          call bcs(q3d(ib,jb,kb,4))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,4),                               &
                     qw1(1,1,1,4),qw2(1,1,1,4),qe1(1,1,1,4),qe2(1,1,1,4),   &
                     qs1(1,1,1,4),qs2(1,1,1,4),qn1(1,1,1,4),qn2(1,1,1,4),   &
                     reqs_q(1,4))
#endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,5),qten(ib,jb,kb,5))
          call bcs(q3d(ib,jb,kb,5))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,5),                               &
                     qw1(1,1,1,5),qw2(1,1,1,5),qe1(1,1,1,5),qe2(1,1,1,5),   &
                     qs1(1,1,1,5),qs2(1,1,1,5),qn1(1,1,1,5),qn2(1,1,1,5),   &
                     reqs_q(1,5))
#endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,6),qten(ib,jb,kb,6))
          call bcs(q3d(ib,jb,kb,6))
#ifdef MPI
          call comm_3s_start(q3d(ib,jb,kb,6),                               &
                     qw1(1,1,1,6),qw2(1,1,1,6),qe1(1,1,1,6),qe2(1,1,1,6),   &
                     qs1(1,1,1,6),qs2(1,1,1,6),qn1(1,1,1,6),qn2(1,1,1,6),   &
                     reqs_q(1,6))
#endif
!-----------------------------------------------------------------------
!  Thompson scheme
        ELSEIF(ptype.eq.3)THEN
          simple_comm = .true.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-12,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
          call pdefq(1.0e-12,asq(3),ruh,rvh,rmh,rho,q3d(ib,jb,kb,3))
          call pdefq(1.0e-12,asq(4),ruh,rvh,rmh,rho,q3d(ib,jb,kb,4))
          call pdefq(1.0e-12,asq(5),ruh,rvh,rmh,rho,q3d(ib,jb,kb,5))
          call pdefq(1.0e-12,asq(6),ruh,rvh,rmh,rho,q3d(ib,jb,kb,6))
!!!          call pdefq(    1.0,asq(7),ruh,rvh,rmh,rho,q3d(ib,jb,kb,7))
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            dum1(i,j,k)=pi0(i,j,k)+pp3d(i,j,k)
            dum2(i,j,k)=dz*rmh(i,j,k)
          enddo
          enddo
          enddo
          call mp_gt_driver(q3d(ib,jb,kb,1),q3d(ib,jb,kb,2),q3d(ib,jb,kb,3), &
                            q3d(ib,jb,kb,4),q3d(ib,jb,kb,5),q3d(ib,jb,kb,6), &
                            q3d(ib,jb,kb,7),th0,th3d,dum1,prs,dum2,dtl,rain)
!-----------------------------------------------------------------------
!  LFOICE scheme

        ELSEIF(ptype.eq.4)THEN
          simple_comm = .true.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-14,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
          call pdefq(1.0e-14,asq(3),ruh,rvh,rmh,rho,q3d(ib,jb,kb,3))
          call pdefq(1.0e-14,asq(4),ruh,rvh,rmh,rho,q3d(ib,jb,kb,4))
          call pdefq(1.0e-14,asq(5),ruh,rvh,rmh,rho,q3d(ib,jb,kb,5))
          call pdefq(1.0e-14,asq(6),ruh,rvh,rmh,rho,q3d(ib,jb,kb,6))
          call lfo_ice_drive(mf, pi0, prs0, pp3d, prs, th0, th3d,    &
                             qv0, rho0, q3d, qten, dum1)
          do n=2,numq
            call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                         q3d(ib,jb,kb,n),qten(ib,jb,kb,n))
          enddo

!-----------------------------------------------------------------------
!  Morrison scheme

        ELSEIF(ptype.eq.5)THEN
          simple_comm = .true.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-12,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
          call pdefq(1.0e-12,asq(3),ruh,rvh,rmh,rho,q3d(ib,jb,kb,3))
          call pdefq(1.0e-12,asq(4),ruh,rvh,rmh,rho,q3d(ib,jb,kb,4))
          call pdefq(1.0e-12,asq(5),ruh,rvh,rmh,rho,q3d(ib,jb,kb,5))
          call pdefq(1.0e-12,asq(6),ruh,rvh,rmh,rho,q3d(ib,jb,kb,6))
!!!          call pdefq(    1.0,asq(7),ruh,rvh,rmh,rho,q3d(ib,jb,kb,7))
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            ! dum1 = pi
            ! dum2 = dz
            ! dum3 = theta
            dum1(i,j,k)=pi0(i,j,k)+pp3d(i,j,k)
            dum2(i,j,k)=dz*rmh(i,j,k)
            dum3(i,j,k)=th0(i,j,k)+th3d(i,j,k)
          enddo
          enddo
          enddo
          call mp_graupel(nstep,dum3,                                         &
                          q3d(ib,jb,kb, 1),q3d(ib,jb,kb, 2),q3d(ib,jb,kb, 3), &
                          q3d(ib,jb,kb, 4),q3d(ib,jb,kb, 5),q3d(ib,jb,kb, 6), &
                          q3d(ib,jb,kb, 7),q3d(ib,jb,kb, 8),q3d(ib,jb,kb, 9), &
                          q3d(ib,jb,kb,10),q3d(ib,jb,kb,11),                  &
                          dum1,prs,dtl,dum2,w3d,rain,dum4,qbudget(1),ruh,rvh)
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            th3d(i,j,k)=dum3(i,j,k)-th0(i,j,k)
          enddo
          enddo
          enddo

!-----------------------------------------------------------------------
!  RE87 approach

        ELSEIF(ptype.eq.6)THEN
          simple_comm = .true.
          call pdefq(    0.0,asq(1),ruh,rvh,rmh,rho,q3d(ib,jb,kb,1))
          call pdefq(1.0e-14,asq(2),ruh,rvh,rmh,rho,q3d(ib,jb,kb,2))
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            if(q3d(i,j,k,2).gt.0.001)then
              qten(i,j,k,2) = v_t
            else
              qten(i,j,k,2) = 0.0
            endif
          enddo
          enddo
          enddo
          call geterain(cpl,lv1,qbudget(7),ruh,rvh,dum1,rho,q3d(ib,jb,kb,2),qten(ib,jb,kb,2))
          if(efall.ge.1)then
            call getefall(cpl,ruh,rvh,mf,pi0,th0,dum1,dum2,rho,   &
                          pp3d,th3d,q3d(ib,jb,kb,2),qten(ib,jb,kb,2))
          endif
          call fallout(qbudget(6),ruh,rvh,zh,mh,mf,rain,rho,   &
                       q3d(ib,jb,kb,2),qten(ib,jb,kb,2))
          call satadj(4,qbudget(1),qbudget(2),ruh,rvh,rmh,pi0,th0,   &
                      rho,dum3,pp3d,prs,th3d,                   &
                      q3d(ib,jb,kb,nqv),q3d(ib,jb,kb,2),ql0)

!-----------------------------------------------------------------------
!  insert new microphysics schemes here
!-----------------------------------------------------------------------
!        ELSEIF(ptype.eq.7)THEN
!-----------------------------------------------------------------------
! otherwise, stop for undefined ptype
        ELSE
          print *,'  Undefined ptype!'
          call stopcm1
        ENDIF

!-----------------------------------------------------------------------
!Begin:  message passing for simple_comm
      IF(simple_comm)THEN
          call bcs(th3d)
          call bcs(pp3d)
          DO n=1,numq
            call bcs(q3d(ib,jb,kb,n))
          ENDDO
#ifdef MPI
          call comm_3s_start(th3d,sw31,sw32,se31,se32,   &
                                  ss31,ss32,sn31,sn32,reqs_s)
          call comm_3s_start(pp3d,pw31,pw32,pe31,pe32,   &
                                  ps31,ps32,pn31,pn32,reqs_p)
          DO n=1,numq
            call comm_3s_start(q3d(ib,jb,kb,n),                               &
                       qw1(1,1,1,n),qw2(1,1,1,n),qe1(1,1,1,n),qe2(1,1,1,n),   &
                       qs1(1,1,1,n),qs2(1,1,1,n),qn1(1,1,1,n),qn2(1,1,1,n),   &
                       reqs_q(1,n))
          ENDDO
#endif
      ENDIF
!Done:  message passing for simple_comm
!-----------------------------------------------------------------------

      ENDIF

!-----------------------------------------------------------------
!  Equate the two arrays

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je
      do i=ib,ie+1
        ua(i,j,k)=u3d(i,j,k)
      enddo
      enddo
      enddo
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je+1
      do i=ib,ie
        va(i,j,k)=v3d(i,j,k)
      enddo
      enddo
      enddo
 
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke+1
      do j=jb,je
      do i=ib,ie
        wa(i,j,k)=w3d(i,j,k)
      enddo
      enddo
      enddo

      if(iturb.eq.1)then
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
        do k=kbt,ket
        do j=jbt,jet
        do i=ibt,iet
          tkea(i,j,k)=tke3d(i,j,k)
        enddo
        enddo
        enddo
      endif

      if(iptra.eq.1)then
        do n=1,npt
          call pdefq(0.0,afoo,ruh,rvh,rmh,rho,pt3d(ib,jb,kb,n))
          call bcs(pt3d(ib,jb,kb,n))
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=kb,ke
          do j=jb,je
          do i=ib,ie
            pta(i,j,k,n)=pt3d(i,j,k,n)
          enddo
          enddo
          enddo
        enddo
      endif
      if(timestats.ge.1) time_integ=time_integ+mytime()
 
#ifndef MPI
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je
      do i=ib,ie
        ppi(i,j,k)=pp3d(i,j,k)
      enddo
      enddo
      enddo

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je
      do i=ib,ie
        tha(i,j,k)=th3d(i,j,k)
      enddo
      enddo
      enddo

      if(imoist.eq.1)then
 
        do n=1,numq

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=kb,ke
          do j=jb,je
          do i=ib,ie
            qa(i,j,k,n)=q3d(i,j,k,n)
          enddo
          enddo
          enddo

        enddo
 
      endif

      if(timestats.ge.1) time_integ=time_integ+mytime()
#endif

!---- finish communication of moisture ----

#ifdef MPI
  IF(imoist.eq.1)THEN
    IF(simple_comm)THEN
      call fincomm(tha,th3d,sw31,sw32,se31,se32,   &
                            ss31,ss32,sn31,sn32,reqs_s)
      call fincomm(ppi,pp3d,pw31,pw32,pe31,pe32,   &
                            ps31,ps32,pn31,pn32,reqs_p)
      do n=1,numq
        call fincomm(qa(ib,jb,kb,n),q3d(ib,jb,kb,n),                        &
                     qw1(1,1,1,n),qw2(1,1,1,n),qe1(1,1,1,n),qe2(1,1,1,n),   &
                     qs1(1,1,1,n),qs2(1,1,1,n),qn1(1,1,1,n),qn2(1,1,1,n),   &
                     reqs_q(1,n))
      enddo
    ELSE
      if(ptype.eq.1)then
        call fincomm(qa(ib,jb,kb,3),q3d(ib,jb,kb,3),                        &
                     qw1(1,1,1,3),qw2(1,1,1,3),qe1(1,1,1,3),qe2(1,1,1,3),   &
                     qs1(1,1,1,3),qs2(1,1,1,3),qn1(1,1,1,3),qn2(1,1,1,3),   &
                     reqs_q(1,3))
        call fincomm(qa(ib,jb,kb,1),q3d(ib,jb,kb,1),                        &
                     qw1(1,1,1,1),qw2(1,1,1,1),qe1(1,1,1,1),qe2(1,1,1,1),   &
                     qs1(1,1,1,1),qs2(1,1,1,1),qn1(1,1,1,1),qn2(1,1,1,1),   &
                     reqs_q(1,1))
        call fincomm(qa(ib,jb,kb,2),q3d(ib,jb,kb,2),                        &
                     qw1(1,1,1,2),qw2(1,1,1,2),qe1(1,1,1,2),qe2(1,1,1,2),   &
                     qs1(1,1,1,2),qs2(1,1,1,2),qn1(1,1,1,2),qn2(1,1,1,2),   &
                     reqs_q(1,2))
        call fincomm(tha,th3d,sw31,sw32,se31,se32,   &
                              ss31,ss32,sn31,sn32,reqs_s)
        call fincomm(ppi,pp3d,pw31,pw32,pe31,pe32,   &
                              ps31,ps32,pn31,pn32,reqs_p)
      elseif(ptype.eq.2)then
        call fincomm(qa(ib,jb,kb,1),q3d(ib,jb,kb,1),                        &
                     qw1(1,1,1,1),qw2(1,1,1,1),qe1(1,1,1,1),qe2(1,1,1,1),   &
                     qs1(1,1,1,1),qs2(1,1,1,1),qn1(1,1,1,1),qn2(1,1,1,1),   &
                     reqs_q(1,1))
        call fincomm(qa(ib,jb,kb,2),q3d(ib,jb,kb,2),                        &
                     qw1(1,1,1,2),qw2(1,1,1,2),qe1(1,1,1,2),qe2(1,1,1,2),   &
                     qs1(1,1,1,2),qs2(1,1,1,2),qn1(1,1,1,2),qn2(1,1,1,2),   &
                     reqs_q(1,2))
        call fincomm(tha,th3d,sw31,sw32,se31,se32,   &
                              ss31,ss32,sn31,sn32,reqs_s)
        call fincomm(ppi,pp3d,pw31,pw32,pe31,pe32,   &
                              ps31,ps32,pn31,pn32,reqs_p)
        call fincomm(qa(ib,jb,kb,3),q3d(ib,jb,kb,3),                        &
                     qw1(1,1,1,3),qw2(1,1,1,3),qe1(1,1,1,3),qe2(1,1,1,3),   &
                     qs1(1,1,1,3),qs2(1,1,1,3),qn1(1,1,1,3),qn2(1,1,1,3),   &
                     reqs_q(1,3))
        call fincomm(qa(ib,jb,kb,4),q3d(ib,jb,kb,4),                        &
                     qw1(1,1,1,4),qw2(1,1,1,4),qe1(1,1,1,4),qe2(1,1,1,4),   &
                     qs1(1,1,1,4),qs2(1,1,1,4),qn1(1,1,1,4),qn2(1,1,1,4),   &
                     reqs_q(1,4))
        call fincomm(qa(ib,jb,kb,5),q3d(ib,jb,kb,5),                        &
                     qw1(1,1,1,5),qw2(1,1,1,5),qe1(1,1,1,5),qe2(1,1,1,5),   &
                     qs1(1,1,1,5),qs2(1,1,1,5),qn1(1,1,1,5),qn2(1,1,1,5),   &
                     reqs_q(1,5))
        call fincomm(qa(ib,jb,kb,6),q3d(ib,jb,kb,6),                        &
                     qw1(1,1,1,6),qw2(1,1,1,6),qe1(1,1,1,6),qe2(1,1,1,6),   &
                     qs1(1,1,1,6),qs2(1,1,1,6),qn1(1,1,1,6),qn2(1,1,1,6),   &
                     reqs_q(1,6))
      endif
    ENDIF

  ELSE

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je
      do i=ib,ie
        ppi(i,j,k)=pp3d(i,j,k)
      enddo
      enddo
      enddo

!$omp parallel do default(shared)  &
!$omp private(i,j,k)
      do k=kb,ke
      do j=jb,je
      do i=ib,ie
        tha(i,j,k)=th3d(i,j,k)
      enddo
      enddo
      enddo
      if(timestats.ge.1) time_integ=time_integ+mytime()

  ENDIF
#endif

      if(imove.eq.1.and.imoist.eq.1)then
        call movesfc(uh,vh,rain(ib,jb,2),dum1(ib,jb,1),dum1(ib,jb,2),dum1(ib,jb,3))
      endif

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cc   All done   cccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

#ifdef MPI
      call MPI_BARRIER (MPI_COMM_WORLD,ierr)
      if(timestats.ge.1) time_mpb=time_mpb+mytime()
#endif

!--------------------------------------------------------------------
!  Get statistics

!$omp parallel do default(shared)  &
!$omp private(i,j,n,tem)
      do j=1,nj
      do i=1,ni
        tem = sqrt( (umove+0.5*(ua(i,j,1)+ua(i+1,j,1)))**2    &
                   +(vmove+0.5*(va(i,j,1)+va(i,j+1,1)))**2 ) 
        do n=1,nrain
          sws(i,j,n)=max(sws(i,j,n),tem)
        enddo
      enddo
      enddo

      if(imove.eq.1)then
        call movesfc(uh,vh,sws(ib,jb,2),dum1(ib,jb,1),dum1(ib,jb,2),dum1(ib,jb,3))
      endif

      if(mod(nstep,statfrq).eq.0)then
        IF(axisymm.eq.0)THEN
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            ppten(i,j,k)=rho(i,j,k)
          enddo
          enddo
          enddo
        ELSE
!$omp parallel do default(shared)  &
!$omp private(i,j,k)
          do k=1,nk
          do j=1,nj
          do i=1,ni
            ppten(i,j,k)=rho(i,j,k)*2.0*pi*xh(i)
          enddo
          enddo
          enddo
        ENDIF
        call statpack(nrec,cloudvar,qname,budname,qbudget,asq,bsq,      &
                      xh,rxh,uh,ruh,xf,uf,yh,vh,rvh,vf,zh,mh,rmh,mf,    &
                      pi0,rho0,thv0,th0,qv0,u0,v0,                      &
                      dum1,dum2,dum3,dum4,divx,ppten,prs,               &
                      ua,va,wa,ppi,tha,qa,qten,kmh,kmv,khh,khv,tkea,pta)
      endif

!--------------------------------------------------------------------
!  Writeout and stuff

      rtime=float(nstep)*dtl
      if(myid.eq.0)then
        if(timeformat.eq.1)then
          write(6,110) nstep,rtime,' sec '
        elseif(timeformat.eq.2)then
          write(6,110) nstep,rtime/60.0,' min '
        elseif(timeformat.eq.3)then
          write(6,110) nstep,rtime/3600.0,' hour'
        elseif(timeformat.eq.4)then
          write(6,110) nstep,rtime/86400.0,' day '
        else
          write(6,110) nstep,rtime,' sec'
        endif
110     format(2x,i12,4x,f18.6,a5)
      endif
      if(timestats.ge.1) time_misc=time_misc+mytime()

      if( (rtime.ge.sngl(taptim)).or.stopit )then
        nwrite=nwrite+1
      IF(output_format.eq.1)THEN
        call writeout(51,nwrite,qname,sigma,zh,pi0,prs0,rho0,th0,qv0,u0,v0,  &
                      zs,rain,sws,thflux,qvflux,dum1,dum2,              &
                      rho,prs,ua,uten,va,vten,wa,wten,ppi,tha,          &
                      qa,kmh,kmv,khh,khv,tkea,pta)
        if(terrain_flag .and. output_interp.eq.1)then
          call writeout(71,nwrite,qname,sigma,zh,pi0,prs0,rho0,th0,qv0,u0,v0,  &
                        zs,rain,sws,thflux,qvflux,dum1,dum2,              &
                        rho,prs,ua,uten,va,vten,wa,wten,ppi,tha,          &
                        qa,kmh,kmv,khh,khv,tkea,pta)
        endif
#ifdef NETCDF
      ELSEIF(output_format.eq.2)THEN
        call writeout_cdf(nwrite,qname,sigma,sigmaf,xh,xf,yh,yf,zh,zf, &
                      pi0,prs0,rho0,th0,qv0,u0,v0,                     &
                      zs,rain,sws,thflux,qvflux,dum1,dum2,             &
                      rho,prs,ua,uten,va,vten,wa,wten,ppi,tha,         &
                      qa,kmh,kmv,khh,khv,tkea,pta)
#endif
#ifdef HDFOUT

        else if(output_format.eq.3) then
        call writeout_mult_hdf(rtime,.true.,qname,sigma,xf,xh,yf,yh,zf,zh,pi0,rho0,th0,qv0,u0,v0,  &
                      zs,rain,thflux,qvflux,sws,dum1,dum2,              &
                      rho,prs,ua,uten,va,vten,wa,wten,ppi,tha,          &
                      qa,kmh,kmv,khh,khv,tkea,pta)
        else if(output_format.eq.4) then
        call writeout_mult_hdf(rtime,.false.,qname,sigma,xf,xh,yf,yh,zf,zh,pi0,rho0,th0,qv0,u0,v0,  &
                      zs,rain,thflux,qvflux,sws,dum1,dum2,              &
                      rho,prs,ua,uten,va,vten,wa,wten,ppi,tha,          &
                      qa,kmh,kmv,khh,khv,tkea,pta)
#endif
      ENDIF
        taptim=taptim+tapfrq
        if(timestats.ge.1) time_write=time_write+mytime()
      endif

      if(rtime.ge.rsttim .and. rstfrq.gt.0)then
        nrst=nrst+1
        call write_restart(nrec,prec,nwrite,nrst,                    &
                           qbudget,asq,bsq,                          &
                           rain,sws,radbcw,radbce,radbcs,radbcn,     &
                           ua,va,wa,ppi,tha,qa,tkea,pta,pdata,rtime)
        rsttim=rsttim+rstfrq
        if(timestats.ge.1) time_write=time_write+mytime()
      endif

!-------------------------------------------------------------------
!  Parcel update (final time step)

      if( (iprcl.eq.1).and.(nstep.eq.nloop2) )then
        write(outfile,*) '  Calling parcel driver for last step'
        call parcel_driver(prec,xh,uh,ruh,yh,vh,rvh,zh,mh,rmh,mf,        &
                           pi0,thv0,th0,dum1,dum2,dum3,dum4,divx,prs,    &
                           ua,va,wa,ppi,ppten,tha,qa,khv,pdata,rtime)
      endif

!-------------------------------------------------------------------

      if(stopit)then
        if(myid.eq.0)then
          print *
          print *,' Courant number has exceeded 1.5 '
          print *
          print *,' Stopping model .... '
          print *
        endif
        call stopcm1
      endif

      return
      end


