vpdsk(i) = vpdsk(i) + vkdp(i,k)*hybd(k)#ifdef HADVTEST!!jr Set etadot to zero for horizontal advection test!           etadot(i,k+1) = 0.#else           etadot(i,k+1) = etadot(i,k+1) - vpdsk(i) + hybi(k+1)*(ddpn(i)+vpdsn(i))#endif        end do     end if!! Convert eta-dot(dp/deta) to eta-dot!     tmp = etamid(k+1) - etamid(k)     do i = 1,nlon        etadot(i,k+1) = etadot(i,k+1)*tmp/(pmid(i,k+1) - pmid(i,k))     end do  end do!! Zonal and meridional derivatives of the hydrostatic term in the! momentum equations.!  do k = 1,plev     do i = 1,nlon        tmp1     = (1. + zvir*q3(i,k))        tmp2     = t3(i,k)*zvir        rtvl     = rair*( tmp1*tl(i,k) + tmp2*ql(i,k) )        rtvm     = rair*( tmp1*tm(i,k) + tmp2*qm(i,k) )!        tmp      = rpmid(i,k)*pdel(i,k)        hsl(i,k) = -rtvl*tmp        hsm(i,k) = -rtvm*tmp     end do     if(k .ge. nprlev) then        abp0 = ps0*(hyam(k)*(hybi(k+1) - hybi(k)) - hybm(k)*(hyai(k+1) - hyai(k)))        do i = 1,nlon           tmp      = rtv(i,k)*abp0/(pmid(i,k)*pmid(i,k))           hsl(i,k) = hsl(i,k) - pspsl(i)*tmp           hsm(i,k) = hsm(i,k) - pspsm(i)*tmp        end do     endif  end do!! Calculate RHS of all prognostic eqns!  do k = 1,plev     tmp  = cappa*t0(k)*hypi(plevp)/hypm(k)     tmp1 = hypd(k)/hypi(plevp)     do i = 1,nlon!! Surface Pressure eqn!        prhsl(i,k) =  div(i,k)*(tmp1 - pdel(i,k)/ps(i))        psldm(i,k) = -div(i,k)*tmp1 - ed1p(i,k+1) + ed1p(i,k)!! Temperature eqn!        trhsl (i,k) = cappa*tv(i,k)/(1. + cpvir*q3(i,k))*omga(i,k)*rpmid(i,k)        trhsl (i,k) = trhsl(i,k) - omga(i,k)/ps(i)*tmp        tsldm (i,k) = 0.5*(ed1p(i,k+1) + ed1p(i,k))*tmp!! ... horizontal advection portion of Hortal Temperature correction!        trhsl (i,k) = trhsl(i,k) - &                      rcoslat*( u3(i,k)*lpsstarl(i) + v3(i,k)*lpsstarm(i) )*hortalc(k)!! ... Ritchie damping term for Temperature eqn.!        tsld0a(i,k) = rcoslat*( u3l(i,k)*lpsstarl(i) + v3l(i,k)*lpsstarm(i) )!! U/V eqns (includes only the diagonal portion of the hydrostatic term)!        urhsl(i,k) = 0.5*hsl(i,k) + href(k,k)*tl(i,k) + bps(k)*dpsl(i)        vrhsl(i,k) = 0.5*hsm(i,k) + href(k,k)*tm(i,k) + bps(k)*dpsm(i)        usldm(i,k) = -phisl(i) + v3(i,k)*coriol*coslat -href(k,k)*tl(i,k) - bps(k)*dpsl(i)        vsldm(i,k) = -phism(i) - u3(i,k)*coriol*coslat -href(k,k)*tm(i,k) - bps(k)*dpsm(i)     end do!! Add pressure gradient terms to momentum tendencies!     if (k.ge.nprlev) then        do i=1,nlon           tmp        = rtv(i,k)*hybm(k)*rpmid(i,k)           dpslon(i,k) = rcoslat*tmp*pspsl(i)           dpslat(i,k) = rcoslat*tmp*pspsm(i)           urhsl(i,k) = urhsl(i,k) - dpslon(i,k)*coslat           vrhsl(i,k) = vrhsl(i,k) - dpslat(i,k)*coslat        end do     else        do i = 1,nlon           dpslon(i,k) = 0.           dpslat(i,k) = 0.        end do     end if  end do!! Interior levels of the hydrostatic term!  do k=1,plev-1     do l=k+1,plev        do i=1,nlon           urhsl(i,k) = urhsl(i,k) + href(l,k)*tl(i,l) + hsl(i,l)           vrhsl(i,k) = vrhsl(i,k) + href(l,k)*tm(i,l) + hsm(i,l)!           usldm(i,k) = usldm(i,k) - href(l,k)*tl(i,l)           vsldm(i,k) = vsldm(i,k) - href(l,k)*tm(i,l)        end do     end do  end do  if(is_first_step()) then     do k = 1,plevp        do i = 1,nlon           etadotm1(i,k) = etadot(i,k)        end do     end do  endif!! The modified etadotm1 will be used later for trajectory calculation in SCANSLT!  do k = 1,plevp     do i = 1,nlon        etadotm1(i,k) = 2.*etadot(i,k) - etadotm1(i,k)     end do  end do!! Compute vertical advection portion of Hortal Temperature correction!  npr = nprlev - 1  if(npr .lt. 1) npr = 1  do k = npr,plev-1     tmpk   = 0.5*hdel(k-1)/detai(k)     tmpkp1 = 0.5*hdel(k  )/detai(k)     do i = 1,nlon        trhsl(i,k) = trhsl(i,k) - (etadot(i,k  )*tmpk + etadot(i,k+1)*tmpkp1)*lpsstar(i)     end do  end do!! ... bottom level!  tmpk = 0.5*hdel(plev-1)/detai(plev)  do i = 1,nlon     trhsl(i,plev) = trhsl(i,plev) - etadot(i,plev)*tmpk*lpsstar(i)  end do!! Compute (by extrapolation) RHS terms for time n + 1/2!  if(is_first_step()) then     do k = 1,plev        do i = 1,nlon           urhs  (i,k) = urhsl (i,k)           vrhs  (i,k) = vrhsl (i,k)           trhs  (i,k) = trhsl (i,k)           prhs  (i,k) = prhsl (i,k)        end do     end do  endif!  do k = 1,plev     do i = 1,nlon        tmp         =       urhsl (i,k)        urhsl (i,k) = facm1*urhsl (i,k) + facm2*urhs  (i,k)        urhs  (i,k) = tmp        tmp         =       vrhsl (i,k)        vrhsl (i,k) = facm1*vrhsl (i,k) + facm2*vrhs  (i,k)        vrhs  (i,k) = tmp        tmp         =       trhsl (i,k)        trhsl (i,k) = facm1*trhsl (i,k) + facm2*trhs  (i,k)        trhs  (i,k) = tmp        tmp         =       prhsl (i,k)        prhsl (i,k) = facm1*prhsl (i,k) + facm2*prhs  (i,k)        prhs  (i,k) = tmp     end do  end do!! Prepare SLD arrays for interpolation by SCANSLT!! Append appropriate coefficients (including decentering epsilon values -- See! Notes 13,14,15)!  dt   = 0.5*ztodt  tmp1 = 0.5*ztodt/coslat  do k = 1,plev     tmp = cappa*t0(k)*hypi(plevp)/hypm(k)     do i = 1,nlon        urhsl (i,k) = tmp1*urhsl (i,k)        vrhsl (i,k) = tmp1*vrhsl (i,k)        trhsl (i,k) = dt  *trhsl (i,k)        tsld0a(i,k) = dt  *tsld0a(i,k)        prhsl (i,k) = dt  *prhsl (i,k)        t2    (i,k) = ztodt*t2   (i,k)        fu    (i,k) = ztodt*fu   (i,k)        fv    (i,k) = ztodt*fv   (i,k)!! Combine terms.! (Time split the midpoint results between arrival and departure!  points)!        tmp2          = lpsstar(i)*hortalc(k)        u3m1    (i,k) = (u3  (i,k)*coslat + onemeps*usldm(i,k)*dt)/coslat + &                                            onemeps*urhsl(i,k) + &                                            onemeps*fu   (i,k)*0.5        v3m1    (i,k) = (v3  (i,k)*coslat + onemeps*vsldm(i,k)*dt)/coslat + &                                            onemeps*vrhsl(i,k) + &                                            onemeps*fv   (i,k)*0.5#ifdef HADVTEST        t3m1    (i,k) = t3  (i,k)#else        t3m1    (i,k) = t3  (i,k) +         onemeps*tsldm(i,k)*dt + &                                            onemeps*trhsl(i,k) - tmp2 + &                                            onemeps*t2(i,k)*0.5#endif        prhssld (i,k) = (psldm(i,k)*dt + prhsl(i,k))*onemeps        tarrsld (i,k) = (trhsl(i,k) + &                         hybm(k)*tmp*tsld0a(i,k))*onepeps + tmp2 + onepeps*t2(i,k)*0.5        parrsld (i,k) = (prhsl(i,k) - hybd(k)*tsld0a(i,k))*onepeps        lnpssld (i,k) = logps(i) - lpsstar(i) - tsld0a(i,k)*onemeps        fu      (i,k) = urhsl(i,k)*coslat*onepeps + fu(i,k)*coslat*onepeps*0.5        fv      (i,k) = vrhsl(i,k)*coslat*onepeps + fv(i,k)*coslat*onepeps*0.5     end do  end do#ifdef VADVTEST  do k=2,plev     do i=1,nlon        etadot(i,k) = -0.5/(12.*86400.)     end do  end do#endif!  returnend subroutine grmult