122     &                              LBi, UBi, LBj, UBj,                 &
 
  123     &                              IminS, ImaxS, JminS, JmaxS,         &
 
  126     &                              pmask, rmask, umask, vmask,         &
 
  129     &                              pmask_wet, rmask_wet,               &
 
  130     &                              umask_wet, vmask_wet,               &
 
  132     &                              om_p, om_r, om_u, om_v,             &
 
  133     &                              on_p, on_r, on_u, on_v,             &
 
  139     &                              visc2_p, visc2_r,                   &
 
  142     &                              DiaRUfrc, DiaRVfrc,                 &
 
  143     &                              DiaU3wrk, DiaV3wrk,                 &
 
  154      integer, 
intent(in) :: ng, tile
 
  155      integer, 
intent(in) :: LBi, UBi, LBj, UBj
 
  156      integer, 
intent(in) :: IminS, ImaxS, JminS, JmaxS
 
  157      integer, 
intent(in) :: nrhs, nnew
 
  161      real(r8), 
intent(in) :: pmask(LBi:,LBj:)
 
  162      real(r8), 
intent(in) :: rmask(LBi:,LBj:)
 
  163      real(r8), 
intent(in) :: umask(LBi:,LBj:)
 
  164      real(r8), 
intent(in) :: vmask(LBi:,LBj:)
 
  167      real(r8), 
intent(in) :: pmask_wet(LBi:,LBj:)
 
  168      real(r8), 
intent(in) :: rmask_wet(LBi:,LBj:)
 
  169      real(r8), 
intent(in) :: umask_wet(LBi:,LBj:)
 
  170      real(r8), 
intent(in) :: vmask_wet(LBi:,LBj:)
 
  172      real(r8), 
intent(in) :: om_p(LBi:,LBj:)
 
  173      real(r8), 
intent(in) :: om_r(LBi:,LBj:)
 
  174      real(r8), 
intent(in) :: om_u(LBi:,LBj:)
 
  175      real(r8), 
intent(in) :: om_v(LBi:,LBj:)
 
  176      real(r8), 
intent(in) :: on_p(LBi:,LBj:)
 
  177      real(r8), 
intent(in) :: on_r(LBi:,LBj:)
 
  178      real(r8), 
intent(in) :: on_u(LBi:,LBj:)
 
  179      real(r8), 
intent(in) :: on_v(LBi:,LBj:)
 
  180      real(r8), 
intent(in) :: pm(LBi:,LBj:)
 
  181      real(r8), 
intent(in) :: pn(LBi:,LBj:)
 
  182      real(r8), 
intent(in) :: Hz(LBi:,LBj:,:)
 
  183      real(r8), 
intent(in) :: z_r(LBi:,LBj:,:)
 
  185      real(r8), 
intent(in) :: visc3d_r(LBi:,LBj:,:)
 
  187      real(r8), 
intent(in) :: visc2_p(LBi:,LBj:)
 
  188      real(r8), 
intent(in) :: visc2_r(LBi:,LBj:)
 
  190# ifdef DIAGNOSTICS_UV 
  191      real(r8), 
intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
 
  192      real(r8), 
intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
 
  193      real(r8), 
intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
 
  194      real(r8), 
intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
 
  196      real(r8), 
intent(inout) :: rufrc(LBi:,LBj:)
 
  197      real(r8), 
intent(inout) :: rvfrc(LBi:,LBj:)
 
  198      real(r8), 
intent(inout) :: u(LBi:,LBj:,:,:)
 
  199      real(r8), 
intent(inout) :: v(LBi:,LBj:,:,:)
 
  202      real(r8), 
intent(in) :: pmask(LBi:UBi,LBj:UBj)
 
  203      real(r8), 
intent(in) :: rmask(LBi:UBi,LBj:UBj)
 
  204      real(r8), 
intent(in) :: umask(LBi:UBi,LBj:UBj)
 
  205      real(r8), 
intent(in) :: vmask(LBi:UBi,LBj:UBj)
 
  208      real(r8), 
intent(in) :: pmask_wet(LBi:UBi,LBj:UBj)
 
  209      real(r8), 
intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
 
  210      real(r8), 
intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
 
  211      real(r8), 
intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
 
  213      real(r8), 
intent(in) :: om_p(LBi:UBi,LBj:UBj)
 
  214      real(r8), 
intent(in) :: om_r(LBi:UBi,LBj:UBj)
 
  215      real(r8), 
intent(in) :: om_u(LBi:UBi,LBj:UBj)
 
  216      real(r8), 
intent(in) :: om_v(LBi:UBi,LBj:UBj)
 
  217      real(r8), 
intent(in) :: on_p(LBi:UBi,LBj:UBj)
 
  218      real(r8), 
intent(in) :: on_r(LBi:UBi,LBj:UBj)
 
  219      real(r8), 
intent(in) :: on_u(LBi:UBi,LBj:UBj)
 
  220      real(r8), 
intent(in) :: on_v(LBi:UBi,LBj:UBj)
 
  221      real(r8), 
intent(in) :: pm(LBi:UBi,LBj:UBj)
 
  222      real(r8), 
intent(in) :: pn(LBi:UBi,LBj:UBj)
 
  223      real(r8), 
intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
 
  224      real(r8), 
intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
 
  226      real(r8), 
intent(in) :: visc3d_r(LBi:UBi,LBj:UBj,N(ng))
 
  228      real(r8), 
intent(in) :: visc2_p(LBi:UBi,LBj:UBj)
 
  229      real(r8), 
intent(in) :: visc2_r(LBi:UBi,LBj:UBj)
 
  231# ifdef DIAGNOSTICS_UV 
  232      real(r8), 
intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
 
  233      real(r8), 
intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
 
  234      real(r8), 
intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
 
  235      real(r8), 
intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
 
  237      real(r8), 
intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
 
  238      real(r8), 
intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
 
  239      real(r8), 
intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
 
  240      real(r8), 
intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
 
  245      integer :: i, j, k, k1, k2
 
  247      real(r8) :: cff, fac1, fac2, pm_p, pn_p
 
  248      real(r8) :: cff1, cff2, cff3, cff4
 
  249      real(r8) :: cff5, cff6, cff7, cff8
 
  250      real(r8) :: dmUdz, dnUdz, dmVdz, dnVdz
 
  254      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
 
  255      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
 
  256      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
 
  257      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
 
  259      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFse
 
  260      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFsx
 
  261      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFse
 
  262      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFsx
 
  263      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmUde
 
  264      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmVde
 
  265      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnUdx
 
  266      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnVdx
 
  267      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dUdz
 
  268      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dVdz
 
  269      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
 
  270      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
 
  271      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
 
  272      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
 
  274#include "set_bounds.h" 
  293      k_loop : 
DO k=0,n(ng)
 
  302              cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
  307              cff=cff*umask_wet(i,j)
 
  309              ufx(i,j)=cff*(z_r(i  ,j,k+1)-                             &
 
  315              cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
  320              cff=cff*vmask_wet(i,j)
 
  322              vfe(i,j)=cff*(z_r(i,j  ,k+1)-                             &
 
  329              dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+                        &
 
  331              dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+                        &
 
  337              dzdx_r(i,j,k2)=0.5_r8*(ufx(i  ,j)+                        &
 
  339              dzde_r(i,j,k2)=0.5_r8*(vfe(i,j  )+                        &
 
  353              cff=cff*rmask_wet(i,j)
 
  355              dnudx(i,j,k2)=cff*((pn(i  ,j)+pn(i+1,j))*                 &
 
  356     &                           u(i+1,j,k+1,nrhs)-                     &
 
  357     &                           (pn(i-1,j)+pn(i  ,j))*                 &
 
  364              cff=0.125_r8*(pn(i-1,j  )+pn(i,j  )+                      &
 
  365     &                      pn(i-1,j-1)+pn(i,j-1))
 
  370              cff=cff*pmask_wet(i,j)
 
  372              dmude(i,j,k2)=cff*((pm(i-1,j  )+pm(i,j  ))*               &
 
  373     &                           u(i,j  ,k+1,nrhs)-                     &
 
  374     &                           (pm(i-1,j-1)+pm(i,j-1))*               &
 
  381              cff=0.125_r8*(pm(i-1,j  )+pm(i,j  )+                      &
 
  382     &                      pm(i-1,j-1)+pm(i,j-1))
 
  387              cff=cff*pmask_wet(i,j)
 
  389              dnvdx(i,j,k2)=cff*((pn(i  ,j-1)+pn(i  ,j))*               &
 
  390     &                           v(i  ,j,k+1,nrhs)-                     &
 
  391     &                           (pn(i-1,j-1)+pn(i-1,j))*               &
 
  403              cff=cff*rmask_wet(i,j)
 
  405              dmvde(i,j,k2)=cff*((pm(i,j  )+pm(i,j+1))*                 &
 
  406     &                           v(i,j+1,k+1,nrhs)-                     &
 
  407     &                           (pm(i,j-1)+pm(i,j  ))*                 &
 
  413        IF ((k.eq.0).or.(k.eq.n(ng))) 
THEN 
  440              cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)-z_r(i-1,j,k)+          &
 
  441     &                            z_r(i  ,j,k+1)-z_r(i  ,j,k)))
 
  442              dudz(i,j,k2)=cff*(u(i,j,k+1,nrhs)-                        &
 
  449              cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)-z_r(i,j-1,k)+          &
 
  450     &                            z_r(i,j  ,k+1)-z_r(i,j  ,k)))
 
  451              dvdz(i,j,k2)=cff*(v(i,j,k+1,nrhs)-                        &
 
  463              cff1=min(dzdx_r(i,j,k1),0.0_r8)
 
  464              cff2=max(dzdx_r(i,j,k1),0.0_r8)
 
  465              cff3=min(dzde_r(i,j,k1),0.0_r8)
 
  466              cff4=max(dzde_r(i,j,k1),0.0_r8)
 
  468     &            (on_r(i,j)*(dnudx(i,j,k1)-                            &
 
  470     &                        (cff1*(dudz(i  ,j,k1)+                    &
 
  472     &                         cff2*(dudz(i  ,j,k2)+                    &
 
  473     &                               dudz(i+1,j,k1))))-                 &
 
  474     &             om_r(i,j)*(dmvde(i,j,k1)-                            &
 
  476     &                        (cff3*(dvdz(i,j  ,k1)+                    &
 
  478     &                         cff4*(dvdz(i,j  ,k2)+                    &
 
  484              cff=cff*rmask_wet(i,j)
 
  487              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
 
  488              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
 
  490              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc2_r(i,j)*cff
 
  491              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc2_r(i,j)*cff
 
  498              pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+                      &
 
  499     &                      pm(i  ,j-1)+pm(i  ,j))
 
  500              pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+                      &
 
  501     &                      pn(i  ,j-1)+pn(i  ,j))
 
  502              cff1=min(dzdx_p(i,j,k1),0.0_r8)
 
  503              cff2=max(dzdx_p(i,j,k1),0.0_r8)
 
  504              cff3=min(dzde_p(i,j,k1),0.0_r8)
 
  505              cff4=max(dzde_p(i,j,k1),0.0_r8)
 
  507     &            (hz(i-1,j  ,k)+hz(i,j  ,k)+                           &
 
  508     &             hz(i-1,j-1,k)+hz(i,j-1,k))*                          &
 
  509     &            (on_p(i,j)*(dnvdx(i,j,k1)-                            &
 
  511     &                        (cff1*(dvdz(i-1,j,k1)+                    &
 
  513     &                         cff2*(dvdz(i-1,j,k2)+                    &
 
  514     &                               dvdz(i  ,j,k1))))+                 &
 
  515     &             om_p(i,j)*(dmude(i,j,k1)-                            &
 
  517     &                        (cff3*(dudz(i,j-1,k1)+                    &
 
  519     &                         cff4*(dudz(i,j-1,k2)+                    &
 
  525              cff=cff*pmask_wet(i,j)
 
  529     &               (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+            &
 
  530     &                visc3d_r(i  ,j-1,k)+visc3d_r(i  ,j,k))
 
  531              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
 
  532              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
 
  534              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc2_p(i,j)*cff
 
  535              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc2_p(i,j)*cff
 
  548     &              (visc3d_r(i-1,j,k  )+visc3d_r(i,j,k  )+             &
 
  549     &               visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
 
  553                cff=0.25_r8*(visc2_r(i-1,j)+visc2_r(i,j))
 
  557                cff=0.5_r8*(pn(i-1,j)+pn(i,j))
 
  558                dnudz=cff*dudz(i,j,k2)
 
  559                dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
  563                cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
  564                dmudz=cff*dudz(i,j,k2)
 
  565                dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
  570                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
  571                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
  572                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
  573                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
  575     &                       (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+        &
 
  576     &                        cff2*(cff2*dnudz-dnudx(i  ,j,k2))+        &
 
  577     &                        cff3*(cff3*dnudz-dnudx(i-1,j,k2))+        &
 
  578     &                        cff4*(cff4*dnudz-dnudx(i  ,j,k1)))
 
  580                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
  581                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
  582                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
  583                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
  585     &                       (cff1*(cff1*dmudz-dmude(i,j  ,k1))+        &
 
  586     &                        cff2*(cff2*dmudz-dmude(i,j+1,k2))+        &
 
  587     &                        cff3*(cff3*dmudz-dmude(i,j  ,k2))+        &
 
  588     &                        cff4*(cff4*dmudz-dmude(i,j+1,k1)))
 
  590                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
  591                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
  592                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
  593                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
  594                cff5=min(dzdx_p(i,j  ,k1),0.0_r8)
 
  595                cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
 
  596                cff7=max(dzdx_p(i,j  ,k2),0.0_r8)
 
  597                cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
 
  598                ufsx(i,j,k2)=ufsx(i,j,k2)+                              &
 
  600     &                       (cff1*(cff5*dnvdz-dnvdx(i,j  ,k1))+        &
 
  601     &                        cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+        &
 
  602     &                        cff3*(cff7*dnvdz-dnvdx(i,j  ,k2))+        &
 
  603     &                        cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
 
  605                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
  606                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
  607                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
  608                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
  609                cff5=min(dzde_r(i-1,j,k1),0.0_r8)
 
  610                cff6=min(dzde_r(i  ,j,k2),0.0_r8)
 
  611                cff7=max(dzde_r(i-1,j,k2),0.0_r8)
 
  612                cff8=max(dzde_r(i  ,j,k1),0.0_r8)
 
  613                ufse(i,j,k2)=ufse(i,j,k2)-                              &
 
  615     &                       (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+        &
 
  616     &                        cff2*(cff6*dmvdz-dmvde(i  ,j,k2))+        &
 
  617     &                        cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+        &
 
  618     &                        cff4*(cff8*dmvdz-dmvde(i  ,j,k1)))
 
  626     &              (visc3d_r(i,j-1,k  )+visc3d_r(i,j,k  )+             &
 
  627     &               visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
 
  631                cff=0.25_r8*(visc2_r(i,j-1)+visc2_r(i,j))
 
  635                cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
  636                dnudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
  640                dnvdz=cff*dvdz(i,j,k2)
 
  641                cff=0.5_r8*(pm(i,j-1)+pm(i,j))
 
  642                dmudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
  646                dmvdz=cff*dvdz(i,j,k2)
 
  648                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
  649                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
  650                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
  651                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
  653     &                       (cff1*(cff1*dnvdz-dnvdx(i  ,j,k1))+        &
 
  654     &                        cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+        &
 
  655     &                        cff3*(cff3*dnvdz-dnvdx(i  ,j,k2))+        &
 
  656     &                        cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
 
  658                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
  659                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
  660                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
  661                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
  663     &                       (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+        &
 
  664     &                        cff2*(cff2*dmvdz-dmvde(i,j  ,k2))+        &
 
  665     &                        cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+        &
 
  666     &                        cff4*(cff4*dmvdz-dmvde(i,j  ,k1)))
 
  668                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
  669                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
  670                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
  671                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
  672                cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
 
  673                cff6=min(dzdx_r(i,j  ,k2),0.0_r8)
 
  674                cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
 
  675                cff8=max(dzdx_r(i,j  ,k1),0.0_r8)
 
  676                vfsx(i,j,k2)=vfsx(i,j,k2)-                              &
 
  678     &                       (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+        &
 
  679     &                        cff2*(cff6*dnudz-dnudx(i,j  ,k2))+        &
 
  680     &                        cff3*(cff7*dnudz-dnudx(i,j-1,k2))+        &
 
  681     &                        cff4*(cff8*dnudz-dnudx(i,j  ,k1)))
 
  683                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
  684                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
  685                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
  686                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
  687                cff5=min(dzde_p(i  ,j,k1),0.0_r8)
 
  688                cff6=min(dzde_p(i+1,j,k2),0.0_r8)
 
  689                cff7=max(dzde_p(i  ,j,k2),0.0_r8)
 
  690                cff8=max(dzde_p(i+1,j,k1),0.0_r8)
 
  691                vfse(i,j,k2)=vfse(i,j,k2)+                              &
 
  693     &                       (cff1*(cff5*dmudz-dmude(i  ,j,k1))+        &
 
  694     &                        cff2*(cff6*dmudz-dmude(i+1,j,k2))+        &
 
  695     &                        cff3*(cff7*dmudz-dmude(i  ,j,k2))+        &
 
  696     &                        cff4*(cff8*dmudz-dmude(i+1,j,k1)))
 
  711              cff=
dt(ng)*0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
 
  712              cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j  )-ufx(i-1,j))
 
  713              cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i  ,j))
 
  714              cff3=ufsx(i,j,k2)-ufsx(i,j,k1)
 
  715              cff4=ufse(i,j,k2)-ufse(i,j,k1)
 
  717              cff6=
dt(ng)*(cff3+cff4)
 
  718              rufrc(i,j)=rufrc(i,j)+cff1+cff2+cff3+cff4
 
  719              u(i,j,k,nnew)=u(i,j,k,nnew)+cff5+cff6
 
  721              diarufrc(i,j,3,
m2hvis)=diarufrc(i,j,3,
m2hvis)+cff1+cff2+  &
 
  725              diau3wrk(i,j,k,
m3hvis)=cff5+cff6
 
  726              diau3wrk(i,j,k,
m3xvis)=cff*cff1+
dt(ng)*cff3
 
  727              diau3wrk(i,j,k,
m3yvis)=cff*cff2+
dt(ng)*cff4
 
  734              cff=
dt(ng)*0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
 
  735              cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j  ))
 
  736              cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i  ,j)-vfe(i,j-1))
 
  737              cff3=vfsx(i,j,k2)-vfsx(i,j,k1)
 
  738              cff4=vfse(i,j,k2)-vfse(i,j,k1)
 
  740              cff6=
dt(ng)*(cff3+cff4)
 
  741              rvfrc(i,j)=rvfrc(i,j)+cff1-cff2+cff3+cff4
 
  742              v(i,j,k,nnew)=v(i,j,k,nnew)+cff5+cff6
 
  744              diarvfrc(i,j,3,
m2hvis)=diarvfrc(i,j,3,
m2hvis)+cff1-cff2+  &
 
  748              diav3wrk(i,j,k,
m3hvis)=cff5+cff6
 
  749              diav3wrk(i,j,k,
m3xvis)= cff*cff1+
dt(ng)*cff3
 
  750              diav3wrk(i,j,k,
m3yvis)=-cff*cff2+
dt(ng)*cff4