127     &                              LBi, UBi, LBj, UBj,                 &
 
  128     &                              IminS, ImaxS, JminS, JmaxS,         &
 
  131     &                              pmask, rmask, umask, vmask,         &
 
  134     &                              pmask_wet, rmask_wet,               &
 
  135     &                              umask_wet, vmask_wet,               &
 
  137     &                              om_p, om_r, om_u, om_v,             &
 
  138     &                              on_p, on_r, on_u, on_v,             &
 
  142# ifdef UV_U3ADV_SPLIT
 
  143     &                              Uvis3d_r, Vvis3d_r,                 &
 
  148     &                              visc4_p, visc4_r,                   &
 
  151     &                              DiaRUfrc, DiaRVfrc,                 &
 
  152     &                              DiaU3wrk, DiaV3wrk,                 &
 
  164      integer, 
intent(in) :: ng, tile
 
  165      integer, 
intent(in) :: LBi, UBi, LBj, UBj
 
  166      integer, 
intent(in) :: IminS, ImaxS, JminS, JmaxS
 
  167      integer, 
intent(in) :: nrhs, nnew
 
  171      real(r8), 
intent(in) :: pmask(LBi:,LBj:)
 
  172      real(r8), 
intent(in) :: rmask(LBi:,LBj:)
 
  173      real(r8), 
intent(in) :: umask(LBi:,LBj:)
 
  174      real(r8), 
intent(in) :: vmask(LBi:,LBj:)
 
  177      real(r8), 
intent(in) :: pmask_wet(LBi:,LBj:)
 
  178      real(r8), 
intent(in) :: rmask_wet(LBi:,LBj:)
 
  179      real(r8), 
intent(in) :: umask_wet(LBi:,LBj:)
 
  180      real(r8), 
intent(in) :: vmask_wet(LBi:,LBj:)
 
  182      real(r8), 
intent(in) :: om_p(LBi:,LBj:)
 
  183      real(r8), 
intent(in) :: om_r(LBi:,LBj:)
 
  184      real(r8), 
intent(in) :: om_u(LBi:,LBj:)
 
  185      real(r8), 
intent(in) :: om_v(LBi:,LBj:)
 
  186      real(r8), 
intent(in) :: on_p(LBi:,LBj:)
 
  187      real(r8), 
intent(in) :: on_r(LBi:,LBj:)
 
  188      real(r8), 
intent(in) :: on_u(LBi:,LBj:)
 
  189      real(r8), 
intent(in) :: on_v(LBi:,LBj:)
 
  190      real(r8), 
intent(in) :: pm(LBi:,LBj:)
 
  191      real(r8), 
intent(in) :: pn(LBi:,LBj:)
 
  192      real(r8), 
intent(in) :: Hz(LBi:,LBj:,:)
 
  193      real(r8), 
intent(in) :: z_r(LBi:,LBj:,:)
 
  195#  ifdef UV_U3ADV_SPLIT 
  196      real(r8), 
intent(in) :: Uvis3d_r(LBi:,LBj:,:)
 
  197      real(r8), 
intent(in) :: Vvis3d_r(LBi:,LBj:,:)
 
  199      real(r8), 
intent(in) :: visc3d_r(LBi:,LBj:,:)
 
  202      real(r8), 
intent(in) :: visc4_p(LBi:,LBj:)
 
  203      real(r8), 
intent(in) :: visc4_r(LBi:,LBj:)
 
  205# ifdef DIAGNOSTICS_UV 
  206      real(r8), 
intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
 
  207      real(r8), 
intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
 
  208      real(r8), 
intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
 
  209      real(r8), 
intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
 
  211      real(r8), 
intent(inout) :: rufrc(LBi:,LBj:)
 
  212      real(r8), 
intent(inout) :: rvfrc(LBi:,LBj:)
 
  213      real(r8), 
intent(inout) :: u(LBi:,LBj:,:,:)
 
  214      real(r8), 
intent(inout) :: v(LBi:,LBj:,:,:)
 
  217      real(r8), 
intent(in) :: pmask(LBi:UBi,LBj:UBj)
 
  218      real(r8), 
intent(in) :: rmask(LBi:UBi,LBj:UBj)
 
  219      real(r8), 
intent(in) :: umask(LBi:UBi,LBj:UBj)
 
  220      real(r8), 
intent(in) :: vmask(LBi:UBi,LBj:UBj)
 
  223      real(r8), 
intent(in) :: pmask_wet(LBi:UBi,LBj:UBj)
 
  224      real(r8), 
intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
 
  225      real(r8), 
intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
 
  226      real(r8), 
intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
 
  228      real(r8), 
intent(in) :: om_p(LBi:UBi,LBj:UBj)
 
  229      real(r8), 
intent(in) :: om_r(LBi:UBi,LBj:UBj)
 
  230      real(r8), 
intent(in) :: om_u(LBi:UBi,LBj:UBj)
 
  231      real(r8), 
intent(in) :: om_v(LBi:UBi,LBj:UBj)
 
  232      real(r8), 
intent(in) :: on_p(LBi:UBi,LBj:UBj)
 
  233      real(r8), 
intent(in) :: on_r(LBi:UBi,LBj:UBj)
 
  234      real(r8), 
intent(in) :: on_u(LBi:UBi,LBj:UBj)
 
  235      real(r8), 
intent(in) :: on_v(LBi:UBi,LBj:UBj)
 
  236      real(r8), 
intent(in) :: pm(LBi:UBi,LBj:UBj)
 
  237      real(r8), 
intent(in) :: pn(LBi:UBi,LBj:UBj)
 
  238      real(r8), 
intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
 
  239      real(r8), 
intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
 
  241#  ifdef UV_U3ADV_SPLIT 
  242      real(r8), 
intent(in) :: Uvis3d_r(LBi:UBi,LBj:UBj,N(ng))
 
  243      real(r8), 
intent(in) :: Vvis3d_r(LBi:UBi,LBj:UBj,N(ng))
 
  245      real(r8), 
intent(in) :: visc3d_r(LBi:UBi,LBj:UBj,N(ng))
 
  248      real(r8), 
intent(in) :: visc4_p(LBi:UBi,LBj:UBj)
 
  249      real(r8), 
intent(in) :: visc4_r(LBi:UBi,LBj:UBj)
 
  251# ifdef DIAGNOSTICS_UV 
  252      real(r8), 
intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
 
  253      real(r8), 
intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
 
  254      real(r8), 
intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
 
  255      real(r8), 
intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
 
  257      real(r8), 
intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
 
  258      real(r8), 
intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
 
  259      real(r8), 
intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
 
  260      real(r8), 
intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
 
  265      integer :: i, j, k, k1, k2
 
  267      real(r8) :: cff, fac1, fac2, pm_p, pn_p
 
  268      real(r8) :: cff1, cff2, cff3, cff4
 
  269      real(r8) :: cff5, cff6, cff7, cff8
 
  270      real(r8) :: dmUdz, dnUdz, dmVdz, dnVdz
 
  272      real(r8) :: Uvis_p, Vvis_p, visc_p
 
  275      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapU
 
  276      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapV
 
  278      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
 
  279      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
 
  280      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
 
  281      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
 
  283      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFse
 
  284      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFsx
 
  285      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFse
 
  286      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFsx
 
  287      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmUde
 
  288      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmVde
 
  289      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnUdx
 
  290      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnVdx
 
  291      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dUdz
 
  292      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dVdz
 
  293      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
 
  294      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
 
  295      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
 
  296      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
 
  298#include "set_bounds.h" 
  323      k_loop1 : 
DO k=0,n(ng)
 
  332              cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
  337              cff=cff*umask_wet(i,j)
 
  339              ufx(i,j)=cff*(z_r(i  ,j,k+1)-                             &
 
  345              cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
  350              cff=cff*vmask_wet(i,j)
 
  352              vfe(i,j)=cff*(z_r(i,j  ,k+1)-                             &
 
  359              dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+                        &
 
  361              dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+                        &
 
  367              dzdx_r(i,j,k2)=0.5_r8*(ufx(i  ,j)+                        &
 
  369              dzde_r(i,j,k2)=0.5_r8*(vfe(i,j  )+                        &
 
  383              cff=cff*rmask_wet(i,j)
 
  385              dnudx(i,j,k2)=cff*((pn(i  ,j)+pn(i+1,j))*                 &
 
  386     &                           u(i+1,j,k+1,nrhs)-                     &
 
  387     &                           (pn(i-1,j)+pn(i  ,j))*                 &
 
  394              cff=0.125_r8*(pn(i-1,j  )+pn(i,j  )+                      &
 
  395     &                      pn(i-1,j-1)+pn(i,j-1))
 
  400              cff=cff*pmask_wet(i,j)
 
  402              dmude(i,j,k2)=cff*((pm(i-1,j  )+pm(i,j  ))*               &
 
  403     &                           u(i,j  ,k+1,nrhs)-                     &
 
  404     &                           (pm(i-1,j-1)+pm(i,j-1))*               &
 
  411              cff=0.125_r8*(pm(i-1,j  )+pm(i,j  )+                      &
 
  412     &                      pm(i-1,j-1)+pm(i,j-1))
 
  417              cff=cff*pmask_wet(i,j)
 
  419              dnvdx(i,j,k2)=cff*((pn(i  ,j-1)+pn(i  ,j))*               &
 
  420     &                           v(i  ,j,k+1,nrhs)-                     &
 
  421     &                           (pn(i-1,j-1)+pn(i-1,j))*               &
 
  433              cff=cff*rmask_wet(i,j)
 
  435              dmvde(i,j,k2)=cff*((pm(i,j  )+pm(i,j+1))*                 &
 
  436     &                           v(i,j+1,k+1,nrhs)-                     &
 
  437     &                           (pm(i,j-1)+pm(i,j  ))*                 &
 
  443        IF ((k.eq.0).or.(k.eq.n(ng))) 
THEN 
  470              cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)-                       &
 
  474              dudz(i,j,k2)=cff*(u(i,j,k+1,nrhs)-                        &
 
  481              cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)-                       &
 
  485              dvdz(i,j,k2)=cff*(v(i,j,k+1,nrhs)-                        &
 
  497              cff1=min(dzdx_r(i,j,k1),0.0_r8)
 
  498              cff2=max(dzdx_r(i,j,k1),0.0_r8)
 
  499              cff3=min(dzde_r(i,j,k1),0.0_r8)
 
  500              cff4=max(dzde_r(i,j,k1),0.0_r8)
 
  501              cff=on_r(i,j)*(dnudx(i,j,k1)-                             &
 
  503     &                       (cff1*(dudz(i  ,j,k1)+                     &
 
  505     &                        cff2*(dudz(i  ,j,k2)+                     &
 
  506     &                              dudz(i+1,j,k1))))-                  &
 
  507     &            om_r(i,j)*(dmvde(i,j,k1)-                             &
 
  509     &                       (cff3*(dvdz(i,j  ,k1)+                     &
 
  511     &                        cff4*(dvdz(i,j  ,k2)+                     &
 
  517              cff=cff*rmask_wet(i,j)
 
  520# ifdef UV_U3ADV_SPLIT 
  521              ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
 
  522              vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
 
  524              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
 
  525              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
 
  528              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
 
  529              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
 
  536              pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+                      &
 
  537     &                      pm(i  ,j-1)+pm(i  ,j))
 
  538              pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+                      &
 
  539     &                      pn(i  ,j-1)+pn(i  ,j))
 
  540              cff1=min(dzdx_p(i,j,k1),0.0_r8)
 
  541              cff2=max(dzdx_p(i,j,k1),0.0_r8)
 
  542              cff3=min(dzde_p(i,j,k1),0.0_r8)
 
  543              cff4=max(dzde_p(i,j,k1),0.0_r8)
 
  544              cff=on_p(i,j)*(dnvdx(i,j,k1)-                             &
 
  546     &                       (cff1*(dvdz(i-1,j,k1)+                     &
 
  548     &                        cff2*(dvdz(i-1,j,k2)+                     &
 
  549     &                              dvdz(i  ,j,k1))))+                  &
 
  550     &            om_p(i,j)*(dmude(i,j,k1)-                             &
 
  552     &                       (cff3*(dudz(i,j-1,k1)+                     &
 
  554     &                        cff4*(dudz(i,j-1,k2)+                     &
 
  560              cff=cff*pmask_wet(i,j)
 
  563# ifdef UV_U3ADV_SPLIT 
  565     &               (uvis3d_r(i-1,j-1,k)+uvis3d_r(i-1,j,k)+            &
 
  566     &                uvis3d_r(i  ,j-1,k)+uvis3d_r(i  ,j,k))
 
  568     &               (vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+            &
 
  569     &                vvis3d_r(i  ,j-1,k)+vvis3d_r(i  ,j,k))
 
  570              ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
 
  571              vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
 
  574     &               (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+            &
 
  575     &                visc3d_r(i  ,j-1,k)+visc3d_r(i  ,j,k))
 
  576              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
 
  577              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
 
  580              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
 
  581              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
 
  593# ifdef UV_U3ADV_SPLIT 
  595     &              (uvis3d_r(i-1,j,k  )+uvis3d_r(i,j,k  )+             &
 
  596     &               uvis3d_r(i-1,j,k+1)+uvis3d_r(i,j,k+1))
 
  599     &              (visc3d_r(i-1,j,k  )+visc3d_r(i,j,k  )+             &
 
  600     &               visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
 
  605                cff=0.25_r8*(visc4_r(i-1,j)+visc4_r(i,j))
 
  609                cff=0.5_r8*(pn(i-1,j)+pn(i,j))
 
  610                dnudz=cff*dudz(i,j,k2)
 
  611                dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
  615                cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
  616                dmudz=cff*dudz(i,j,k2)
 
  617                dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
  622                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
  623                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
  624                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
  625                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
  627     &                       (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+        &
 
  628     &                        cff2*(cff2*dnudz-dnudx(i  ,j,k2))+        &
 
  629     &                        cff3*(cff3*dnudz-dnudx(i-1,j,k2))+        &
 
  630     &                        cff4*(cff4*dnudz-dnudx(i  ,j,k1)))
 
  632                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
  633                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
  634                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
  635                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
  637     &                       (cff1*(cff1*dmudz-dmude(i,j  ,k1))+        &
 
  638     &                        cff2*(cff2*dmudz-dmude(i,j+1,k2))+        &
 
  639     &                        cff3*(cff3*dmudz-dmude(i,j  ,k2))+        &
 
  640     &                        cff4*(cff4*dmudz-dmude(i,j+1,k1)))
 
  642                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
  643                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
  644                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
  645                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
  646                cff5=min(dzdx_p(i,j  ,k1),0.0_r8)
 
  647                cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
 
  648                cff7=max(dzdx_p(i,j  ,k2),0.0_r8)
 
  649                cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
 
  650                ufsx(i,j,k2)=ufsx(i,j,k2)+                              &
 
  652     &                       (cff1*(cff5*dnvdz-dnvdx(i,j  ,k1))+        &
 
  653     &                        cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+        &
 
  654     &                        cff3*(cff7*dnvdz-dnvdx(i,j  ,k2))+        &
 
  655     &                        cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
 
  657                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
  658                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
  659                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
  660                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
  661                cff5=min(dzde_r(i-1,j,k1),0.0_r8)
 
  662                cff6=min(dzde_r(i  ,j,k2),0.0_r8)
 
  663                cff7=max(dzde_r(i-1,j,k2),0.0_r8)
 
  664                cff8=max(dzde_r(i  ,j,k1),0.0_r8)
 
  665                ufse(i,j,k2)=ufse(i,j,k2)-                              &
 
  667     &                       (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+        &
 
  668     &                        cff2*(cff6*dmvdz-dmvde(i  ,j,k2))+        &
 
  669     &                        cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+        &
 
  670     &                        cff4*(cff8*dmvdz-dmvde(i  ,j,k1)))
 
  677# ifdef UV_U3ADV_SPLIT 
  679     &              (vvis3d_r(i,j-1,k  )+vvis3d_r(i,j,k  )+             &
 
  680     &               vvis3d_r(i,j-1,k+1)+vvis3d_r(i,j,k+1))
 
  683     &              (visc3d_r(i,j-1,k  )+visc3d_r(i,j,k  )+             &
 
  684     &               visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
 
  689                cff=0.25_r8*(visc4_r(i,j-1)+visc4_r(i,j))
 
  693                cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
  694                dnudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
  698                dnvdz=cff*dvdz(i,j,k2)
 
  699                cff=0.5_r8*(pm(i,j-1)+pm(i,j))
 
  700                dmudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
  704                dmvdz=cff*dvdz(i,j,k2)
 
  706                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
  707                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
  708                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
  709                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
  711     &                       (cff1*(cff1*dnvdz-dnvdx(i  ,j,k1))+        &
 
  712     &                        cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+        &
 
  713     &                        cff3*(cff3*dnvdz-dnvdx(i  ,j,k2))+        &
 
  714     &                        cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
 
  716                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
  717                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
  718                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
  719                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
  721     &                       (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+        &
 
  722     &                        cff2*(cff2*dmvdz-dmvde(i,j  ,k2))+        &
 
  723     &                        cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+        &
 
  724     &                        cff4*(cff4*dmvdz-dmvde(i,j  ,k1)))
 
  726                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
  727                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
  728                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
  729                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
  730                cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
 
  731                cff6=min(dzdx_r(i,j  ,k2),0.0_r8)
 
  732                cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
 
  733                cff8=max(dzdx_r(i,j  ,k1),0.0_r8)
 
  734                vfsx(i,j,k2)=vfsx(i,j,k2)-                              &
 
  736     &                       (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+        &
 
  737     &                        cff2*(cff6*dnudz-dnudx(i,j  ,k2))+        &
 
  738     &                        cff3*(cff7*dnudz-dnudx(i,j-1,k2))+        &
 
  739     &                        cff4*(cff8*dnudz-dnudx(i,j  ,k1)))
 
  741                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
  742                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
  743                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
  744                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
  745                cff5=min(dzde_p(i  ,j,k1),0.0_r8)
 
  746                cff6=min(dzde_p(i+1,j,k2),0.0_r8)
 
  747                cff7=max(dzde_p(i  ,j,k2),0.0_r8)
 
  748                cff8=max(dzde_p(i+1,j,k1),0.0_r8)
 
  749                vfse(i,j,k2)=vfse(i,j,k2)+                              &
 
  751     &                       (cff1*(cff5*dmudz-dmude(i  ,j,k1))+        &
 
  752     &                        cff2*(cff6*dmudz-dmude(i+1,j,k2))+        &
 
  753     &                        cff3*(cff7*dmudz-dmude(i  ,j,k2))+        &
 
  754     &                        cff4*(cff8*dmudz-dmude(i+1,j,k1)))
 
  763              cff=0.125_r8*(pm(i-1,j)+pm(i,j))*                         &
 
  764     &                     (pn(i-1,j)+pn(i,j))
 
  765              cff1=1.0_r8/(0.5_r8*(hz(i-1,j,k)+hz(i,j,k)))
 
  766              lapu(i,j,k)=cff*((pn(i-1,j)+pn(i,j))*                     &
 
  767                               (ufx(i,j)-ufx(i-1,j))+                   &
 
  768     &                         (pm(i-1,j)+pm(i,j))*                     &
 
  769     &                         (ufe(i,j+1)-ufe(i,j)))+                  &
 
  770     &                    cff1*((ufsx(i,j,k2)+ufse(i,j,k2))-            &
 
  771     &                          (ufsx(i,j,k1)+ufse(i,j,k1)))
 
  773              lapu(i,j,k)=lapu(i,j,k)*umask(i,j)
 
  776              lapu(i,j,k)=lapu(i,j,k)*umask_wet(i,j)
 
  783              cff=0.125_r8*(pm(i,j)+pm(i,j-1))*                         &
 
  784     &                     (pn(i,j)+pn(i,j-1))
 
  785              cff1=1.0_r8/(0.5_r8*(hz(i,j-1,k)+hz(i,j,k)))
 
  786              lapv(i,j,k)=cff*((pn(i,j-1)+pn(i,j))*                     &
 
  787     &                         (vfx(i+1,j)-vfx(i,j))-                   &
 
  788     &                         (pm(i,j-1)+pm(i,j))*                     &
 
  789     &                         (vfe(i,j)-vfe(i,j-1)))+                  &
 
  790     &                    cff1*((vfsx(i,j,k2)+vfse(i,j,k2))-            &
 
  791     &                          (vfsx(i,j,k1)+vfse(i,j,k1)))
 
  793              lapv(i,j,k)=lapv(i,j,k)*vmask(i,j)
 
  796              lapv(i,j,k)=lapv(i,j,k)*vmask_wet(i,j)
 
  807        IF (
domain(ng)%Western_Edge(tile)) 
THEN 
  811                lapu(istru-1,j,k)=0.0_r8
 
  817                lapu(istru-1,j,k)=lapu(istru,j,k)
 
  824                lapv(istr-1,j,k)=
gamma2(ng)*lapv(istr,j,k)
 
  830                lapv(istr-1,j,k)=0.0_r8
 
  838        IF (
domain(ng)%Eastern_Edge(tile)) 
THEN 
  842                lapu(iend+1,j,k)=0.0_r8
 
  848                lapu(iend+1,j,k)=lapu(iend,j,k)
 
  855                lapv(iend+1,j,k)=
gamma2(ng)*lapv(iend,j,k)
 
  861                lapv(iend+1,j,k)=0.0_r8
 
  869        IF (
domain(ng)%Southern_Edge(tile)) 
THEN 
  873                lapu(i,jstr-1,k)=
gamma2(ng)*lapu(i,jstr,k)
 
  879                lapu(i,jstr-1,k)=0.0_r8
 
  886                lapv(i,jstrv-1,k)=0.0_r8
 
  892                lapv(i,jstrv-1,k)=lapv(i,jstrv,k)
 
  900        IF (
domain(ng)%Northern_Edge(tile)) 
THEN 
  904                lapu(i,jend+1,k)=
gamma2(ng)*lapu(i,jend,k)
 
  910                lapu(i,jend+1,k)=0.0_r8
 
  917                lapv(i,jend+1,k)=0.0_r8
 
  923                lapv(i,jend+1,k)=lapv(i,jend,k)
 
  932        IF (
domain(ng)%SouthWest_Corner(tile)) 
THEN 
  934            lapu(istr  ,jstr-1,k)=0.5_r8*(lapu(istr+1,jstr-1,k)+        &
 
  935     &                                    lapu(istr  ,jstr  ,k))
 
  936            lapv(istr-1,jstr  ,k)=0.5_r8*(lapv(istr-1,jstr+1,k)+        &
 
  937     &                                    lapv(istr  ,jstr  ,k))
 
  944        IF (
domain(ng)%SouthEast_Corner(tile)) 
THEN 
  946            lapu(iend+1,jstr-1,k)=0.5_r8*(lapu(iend  ,jstr-1,k)+        &
 
  947     &                                    lapu(iend+1,jstr  ,k))
 
  948            lapv(iend+1,jstr  ,k)=0.5_r8*(lapv(iend  ,jstr  ,k)+        &
 
  949     &                                    lapv(iend+1,jstr+1,k))
 
  956        IF (
domain(ng)%NorthWest_Corner(tile)) 
THEN 
  958            lapu(istr  ,jend+1,k)=0.5_r8*(lapu(istr+1,jend+1,k)+        &
 
  959     &                                    lapu(istr  ,jend  ,k))
 
  960            lapv(istr-1,jend+1,k)=0.5_r8*(lapv(istr  ,jend+1,k)+        &
 
  961     &                                    lapv(istr-1,jend  ,k))
 
  968        IF (
domain(ng)%NorthEast_Corner(tile)) 
THEN 
  970            lapu(iend+1,jend+1,k)=0.5_r8*(lapu(iend  ,jend+1,k)+        &
 
  971     &                                    lapu(iend+1,jend  ,k))
 
  972            lapv(iend+1,jend+1,k)=0.5_r8*(lapv(iend  ,jend+1,k)+        &
 
  973     &                                    lapv(iend+1,jend  ,k))
 
  982      k_loop2 : 
DO k=0,n(ng)
 
  991              cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
  996              cff=cff*umask_wet(i,j)
 
  998              ufx(i,j)=cff*(z_r(i  ,j,k+1)-                             &
 
 1004              cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
 1009              cff=cff*vmask_wet(i,j)
 
 1011              vfe(i,j)=cff*(z_r(i,j  ,k+1)-                             &
 
 1018              dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+                        &
 
 1020              dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+                        &
 
 1026              dzdx_r(i,j,k2)=0.5_r8*(ufx(i  ,j)+                        &
 
 1028              dzde_r(i,j,k2)=0.5_r8*(vfe(i,j  )+                        &
 
 1043              cff=cff*rmask_wet(i,j)
 
 1045              dnudx(i,j,k2)=cff*((pn(i  ,j)+pn(i+1,j))*                 &
 
 1046     &                           lapu(i+1,j,k+1)-                       &
 
 1047     &                           (pn(i-1,j)+pn(i  ,j))*                 &
 
 1054              cff=0.125_r8*(pn(i-1,j  )+pn(i,j  )+                      &
 
 1055     &                      pn(i-1,j-1)+pn(i,j-1))
 
 1060              cff=cff*pmask_wet(i,j)
 
 1062              dmude(i,j,k2)=cff*((pm(i-1,j  )+pm(i,j  ))*               &
 
 1064     &                           (pm(i-1,j-1)+pm(i,j-1))*               &
 
 1071              cff=0.125_r8*(pm(i-1,j  )+pm(i,j  )+                      &
 
 1072     &                      pm(i-1,j-1)+pm(i,j-1))
 
 1077              cff=cff*pmask_wet(i,j)
 
 1079              dnvdx(i,j,k2)=cff*((pn(i  ,j-1)+pn(i  ,j))*               &
 
 1081     &                           (pn(i-1,j-1)+pn(i-1,j))*               &
 
 1093              cff=cff*rmask_wet(i,j)
 
 1095              dmvde(i,j,k2)=cff*((pm(i,j  )+pm(i,j+1))*                 &
 
 1096     &                           lapv(i,j+1,k+1)-                       &
 
 1097     &                           (pm(i,j-1)+pm(i,j  ))*                 &
 
 1103        IF ((k.eq.0).or.(k.eq.n(ng))) 
THEN 
 1130              cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)-                       &
 
 1134              dudz(i,j,k2)=cff*(lapu(i,j,k+1)-                          &
 
 1140              cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)-                       &
 
 1144              dvdz(i,j,k2)=cff*(lapv(i,j,k+1)-                          &
 
 1156              cff1=min(dzdx_r(i,j,k1),0.0_r8)
 
 1157              cff2=max(dzdx_r(i,j,k1),0.0_r8)
 
 1158              cff3=min(dzde_r(i,j,k1),0.0_r8)
 
 1159              cff4=max(dzde_r(i,j,k1),0.0_r8)
 
 1161     &            (on_r(i,j)*(dnudx(i,j,k1)-                            &
 
 1163     &                        (cff1*(dudz(i  ,j,k1)+                    &
 
 1164     &                               dudz(i+1,j,k2))+                   &
 
 1165     &                         cff2*(dudz(i  ,j,k2)+                    &
 
 1166     &                               dudz(i+1,j,k1))))-                 &
 
 1167     &             om_r(i,j)*(dmvde(i,j,k1)-                            &
 
 1169     &                        (cff3*(dvdz(i,j  ,k1)+                    &
 
 1170     &                               dvdz(i,j+1,k2))+                   &
 
 1171     &                         cff4*(dvdz(i,j  ,k2)+                    &
 
 1172     &                               dvdz(i,j+1,k1)))))
 
 1177              cff=cff*rmask_wet(i,j)
 
 1180# ifdef UV_U3ADV_SPLIT 
 1181              ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
 
 1182              vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
 
 1184              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
 
 1185              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
 
 1188              ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
 
 1189              vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
 
 1196              pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+                      &
 
 1197     &                      pm(i  ,j-1)+pm(i  ,j))
 
 1198              pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+                      &
 
 1199     &                      pn(i  ,j-1)+pn(i  ,j))
 
 1200              cff1=min(dzdx_p(i,j,k1),0.0_r8)
 
 1201              cff2=max(dzdx_p(i,j,k1),0.0_r8)
 
 1202              cff3=min(dzde_p(i,j,k1),0.0_r8)
 
 1203              cff4=max(dzde_p(i,j,k1),0.0_r8)
 
 1205     &            (hz(i-1,j  ,k)+hz(i,j  ,k)+                           &
 
 1206     &             hz(i-1,j-1,k)+hz(i,j-1,k))*                          &
 
 1207     &             (on_p(i,j)*(dnvdx(i,j,k1)-                           &
 
 1209     &                         (cff1*(dvdz(i-1,j,k1)+                   &
 
 1211     &                          cff2*(dvdz(i-1,j,k2)+                   &
 
 1212     &                                dvdz(i  ,j,k1))))+                &
 
 1213     &              om_p(i,j)*(dmude(i,j,k1)-                           &
 
 1215     &                         (cff3*(dudz(i,j-1,k1)+                   &
 
 1217     &                          cff4*(dudz(i,j-1,k2)+                   &
 
 1223              cff=cff*pmask_wet(i,j)
 
 1226# ifdef UV_U3ADV_SPLIT 
 1228     &               (uvis3d_r(i-1,j-1,k)+uvis3d_r(i-1,j,k)+            &
 
 1229     &                uvis3d_r(i  ,j-1,k)+uvis3d_r(i  ,j,k))
 
 1231     &               (vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+            &
 
 1232     &                vvis3d_r(i  ,j-1,k)+vvis3d_r(i  ,j,k))
 
 1233              ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
 
 1234              vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
 
 1237     &               (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+            &
 
 1238     &                visc3d_r(i  ,j-1,k)+visc3d_r(i  ,j,k))
 
 1239              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
 
 1240              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
 
 1243              ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
 
 1244              vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
 
 1252          IF (k.lt.n(ng)) 
THEN 
 1256# ifdef UV_U3ADV_SPLIT 
 1258     &              (uvis3d_r(i-1,j,k  )+uvis3d_r(i,j,k  )+             &
 
 1259     &               uvis3d_r(i-1,j,k+1)+uvis3d_r(i,j,k+1))
 
 1262     &              (visc3d_r(i-1,j,k  )+visc3d_r(i,j,k  )+             &
 
 1263     &               visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
 
 1268                cff=0.25_r8*(visc4_r(i-1,j)+visc4_r(i,j))
 
 1272                cff=0.5_r8*(pn(i-1,j)+pn(i,j))
 
 1273                dnudz=cff*dudz(i,j,k2)
 
 1274                dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
 1275     &                             dvdz(i  ,j+1,k2)+                    &
 
 1276     &                             dvdz(i-1,j  ,k2)+                    &
 
 1278                cff=0.5_r8*(pm(i-1,j)+pm(i,j))
 
 1279                dmudz=cff*dudz(i,j,k2)
 
 1280                dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+                    &
 
 1281     &                             dvdz(i  ,j+1,k2)+                    &
 
 1282     &                             dvdz(i-1,j  ,k2)+                    &
 
 1285                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
 1286                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
 1287                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
 1288                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
 1289                ufsx(i,j,k2)=fac1*                                      &
 
 1290     &                       (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+        &
 
 1291     &                        cff2*(cff2*dnudz-dnudx(i  ,j,k2))+        &
 
 1292     &                        cff3*(cff3*dnudz-dnudx(i-1,j,k2))+        &
 
 1293     &                        cff4*(cff4*dnudz-dnudx(i  ,j,k1)))
 
 1295                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
 1296                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
 1297                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
 1298                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
 1299                ufse(i,j,k2)=fac2*                                      &
 
 1300     &                       (cff1*(cff1*dmudz-dmude(i,j  ,k1))+        &
 
 1301     &                        cff2*(cff2*dmudz-dmude(i,j+1,k2))+        &
 
 1302     &                        cff3*(cff3*dmudz-dmude(i,j  ,k2))+        &
 
 1303     &                        cff4*(cff4*dmudz-dmude(i,j+1,k1)))
 
 1305                cff1=min(dzde_p(i,j  ,k1),0.0_r8)
 
 1306                cff2=min(dzde_p(i,j+1,k2),0.0_r8)
 
 1307                cff3=max(dzde_p(i,j  ,k2),0.0_r8)
 
 1308                cff4=max(dzde_p(i,j+1,k1),0.0_r8)
 
 1309                cff5=min(dzdx_p(i,j  ,k1),0.0_r8)
 
 1310                cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
 
 1311                cff7=max(dzdx_p(i,j  ,k2),0.0_r8)
 
 1312                cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
 
 1313                ufsx(i,j,k2)=ufsx(i,j,k2)+                              &
 
 1315     &                       (cff1*(cff5*dnvdz-dnvdx(i,j  ,k1))+        &
 
 1316     &                        cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+        &
 
 1317     &                        cff3*(cff7*dnvdz-dnvdx(i,j  ,k2))+        &
 
 1318     &                        cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
 
 1320                cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
 
 1321                cff2=min(dzdx_r(i  ,j,k2),0.0_r8)
 
 1322                cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
 
 1323                cff4=max(dzdx_r(i  ,j,k1),0.0_r8)
 
 1324                cff5=min(dzde_r(i-1,j,k1),0.0_r8)
 
 1325                cff6=min(dzde_r(i  ,j,k2),0.0_r8)
 
 1326                cff7=max(dzde_r(i-1,j,k2),0.0_r8)
 
 1327                cff8=max(dzde_r(i  ,j,k1),0.0_r8)
 
 1328                ufse(i,j,k2)=ufse(i,j,k2)-                              &
 
 1330     &                       (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+        &
 
 1331     &                        cff2*(cff6*dmvdz-dmvde(i  ,j,k2))+        &
 
 1332     &                        cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+        &
 
 1333     &                        cff4*(cff8*dmvdz-dmvde(i  ,j,k1)))
 
 1340# ifdef UV_U3ADV_SPLIT 
 1342     &              (vvis3d_r(i,j-1,k  )+vvis3d_r(i,j,k  )+             &
 
 1343     &               vvis3d_r(i,j-1,k+1)+vvis3d_r(i,j,k+1))
 
 1346     &              (visc3d_r(i,j-1,k  )+visc3d_r(i,j,k  )+             &
 
 1347     &               visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
 
 1352                cff=0.25_r8*(visc4_r(i,j-1)+visc4_r(i,j))
 
 1356                cff=0.5_r8*(pn(i,j-1)+pn(i,j))
 
 1357                dnudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
 1358     &                             dudz(i+1,j  ,k2)+                    &
 
 1359     &                             dudz(i  ,j-1,k2)+                    &
 
 1361                dnvdz=cff*dvdz(i,j,k2)
 
 1362                cff=0.5_r8*(pm(i,j-1)+pm(i,j))
 
 1363                dmudz=cff*0.25_r8*(dudz(i  ,j  ,k2)+                    &
 
 1364     &                             dudz(i+1,j  ,k2)+                    &
 
 1365     &                             dudz(i  ,j-1,k2)+                    &
 
 1367                dmvdz=cff*dvdz(i,j,k2)
 
 1369                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
 1370                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
 1371                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
 1372                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
 1373                vfsx(i,j,k2)=fac1*                                      &
 
 1374     &                       (cff1*(cff1*dnvdz-dnvdx(i  ,j,k1))+        &
 
 1375     &                        cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+        &
 
 1376     &                        cff3*(cff3*dnvdz-dnvdx(i  ,j,k2))+        &
 
 1377     &                        cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
 
 1379                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
 1380                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
 1381                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
 1382                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
 1383                vfse(i,j,k2)=fac2*                                      &
 
 1384     &                       (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+        &
 
 1385     &                        cff2*(cff2*dmvdz-dmvde(i,j  ,k2))+        &
 
 1386     &                        cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+        &
 
 1387     &                        cff4*(cff4*dmvdz-dmvde(i,j  ,k1)))
 
 1389                cff1=min(dzde_r(i,j-1,k1),0.0_r8)
 
 1390                cff2=min(dzde_r(i,j  ,k2),0.0_r8)
 
 1391                cff3=max(dzde_r(i,j-1,k2),0.0_r8)
 
 1392                cff4=max(dzde_r(i,j  ,k1),0.0_r8)
 
 1393                cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
 
 1394                cff6=min(dzdx_r(i,j  ,k2),0.0_r8)
 
 1395                cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
 
 1396                cff8=max(dzdx_r(i,j  ,k1),0.0_r8)
 
 1397                vfsx(i,j,k2)=vfsx(i,j,k2)-                              &
 
 1399     &                       (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+        &
 
 1400     &                        cff2*(cff6*dnudz-dnudx(i,j  ,k2))+        &
 
 1401     &                        cff3*(cff7*dnudz-dnudx(i,j-1,k2))+        &
 
 1402     &                        cff4*(cff8*dnudz-dnudx(i,j  ,k1)))
 
 1404                cff1=min(dzdx_p(i  ,j,k1),0.0_r8)
 
 1405                cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
 
 1406                cff3=max(dzdx_p(i  ,j,k2),0.0_r8)
 
 1407                cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
 
 1408                cff5=min(dzde_p(i  ,j,k1),0.0_r8)
 
 1409                cff6=min(dzde_p(i+1,j,k2),0.0_r8)
 
 1410                cff7=max(dzde_p(i  ,j,k2),0.0_r8)
 
 1411                cff8=max(dzde_p(i+1,j,k1),0.0_r8)
 
 1412                vfse(i,j,k2)=vfse(i,j,k2)+                              &
 
 1414     &                       (cff1*(cff5*dmudz-dmude(i  ,j,k1))+        &
 
 1415     &                        cff2*(cff6*dmudz-dmude(i+1,j,k2))+        &
 
 1416     &                        cff3*(cff7*dmudz-dmude(i  ,j,k2))+        &
 
 1417     &                        cff4*(cff8*dmudz-dmude(i+1,j,k1)))
 
 1425#ifdef DIAGNOSTICS_UV 
 1432              cff=
dt(ng)*0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
 
 1433              cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j  )-ufx(i-1,j))
 
 1434              cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i  ,j))
 
 1435              cff3=ufsx(i,j,k2)-ufsx(i,j,k1)
 
 1436              cff4=ufse(i,j,k2)-ufse(i,j,k1)
 
 1437              cff5=cff*(cff1+cff2)
 
 1438              cff6=
dt(ng)*(cff3+cff4)
 
 1439              rufrc(i,j)=rufrc(i,j)-cff1-cff2-cff3-cff4
 
 1440              u(i,j,k,nnew)=u(i,j,k,nnew)-cff5-cff6
 
 1441#ifdef DIAGNOSTICS_UV 
 1442              diarufrc(i,j,3,
m2hvis)=diarufrc(i,j,3,
m2hvis)-cff1-cff2-  &
 
 1444              diarufrc(i,j,3,
m2xvis)=diarufrc(i,j,3,
m2xvis)-cff1-cff3
 
 1445              diarufrc(i,j,3,
m2yvis)=diarufrc(i,j,3,
m2yvis)-cff2-cff4
 
 1446              diau3wrk(i,j,k,
m3hvis)=-cff5-cff6
 
 1447              diau3wrk(i,j,k,
m3xvis)=-cff*cff1-
dt(ng)*cff3
 
 1448              diau3wrk(i,j,k,
m3yvis)=-cff*cff2-
dt(ng)*cff4
 
 1455              cff=
dt(ng)*0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
 
 1456              cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j  ))
 
 1457              cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i  ,j)-vfe(i,j-1))
 
 1458              cff3=vfsx(i,j,k2)-vfsx(i,j,k1)
 
 1459              cff4=vfse(i,j,k2)-vfse(i,j,k1)
 
 1460              cff5=cff*(cff1-cff2)
 
 1461              cff6=
dt(ng)*(cff3+cff4)
 
 1462              rvfrc(i,j)=rvfrc(i,j)-cff1+cff2-cff3-cff4
 
 1463              v(i,j,k,nnew)=v(i,j,k,nnew)-cff5-cff6
 
 1464#ifdef DIAGNOSTICS_UV 
 1465              diarvfrc(i,j,3,
m2hvis)=diarvfrc(i,j,3,
m2hvis)-cff1+cff2-  &
 
 1467              diarvfrc(i,j,3,
m2xvis)=diarvfrc(i,j,3,
m2xvis)-cff1-cff3
 
 1468              diarvfrc(i,j,3,
m2yvis)=diarvfrc(i,j,3,
m2yvis)+cff2-cff4
 
 1469              diav3wrk(i,j,k,
m3hvis)=-cff5-cff6
 
 1470              diav3wrk(i,j,k,
m3xvis)=-cff*cff1-
dt(ng)*cff3
 
 1471              diav3wrk(i,j,k,
m3yvis)= cff*cff2-
dt(ng)*cff4