ROMS_r215: /home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms

00001 !
00002 ! $Id: v3dmix.F,v 1.3 2005/09/22 14:51:39 pmarches Exp $
00003 !
00004 #include "cppdefs.h"
00005 #ifdef SOLVE3D
00006 
00007       subroutine v3dmix (tile)
00008       implicit none
00009       integer tile, trd, omp_get_thread_num
00010 # include "param.h"
00011 # include "private_scratch.h"
00012 # include "compute_tile_bounds.h"
00013       trd=omp_get_thread_num()
00014       call v3dmix_tile (Istr,Iend,Jstr,Jend,  A3d(1,1,trd),
00015      &                                                  A3d(1,2,trd),
00016      &                      A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd),
00017      &                      A2d(1,4,trd), A2d(1,5,trd), A2d(1,6,trd))
00018       return
00019       end
00020 
00021       subroutine v3dmix_tile (Istr,Iend,Jstr,Jend, LapV,VFs, VFz,
00022      &                                     VFx,VFe,wrk,Slope,dVdz)
00023 !
00024 !--------------------------------------------------------------------
00025 !  Compute lateral mixing terms for ETA-direction momentum equation.
00026 !--------------------------------------------------------------------
00027 !
00028       implicit none
00029 # include "param.h"
00030       integer Istr,Iend,Jstr,Jend, i,j,k,indx
00031       real LapV(PRIVATE_2D_SCRATCH_ARRAY,0:N),
00032      &      VFs(PRIVATE_2D_SCRATCH_ARRAY,0:N),
00033      &      VFz(PRIVATE_1D_SCRATCH_ARRAY,0:N),
00034      &          VFx(PRIVATE_2D_SCRATCH_ARRAY),
00035      &          VFe(PRIVATE_2D_SCRATCH_ARRAY),
00036      &          wrk(PRIVATE_2D_SCRATCH_ARRAY),
00037      &        Slope(PRIVATE_2D_SCRATCH_ARRAY),
00038      &         dVdz(PRIVATE_2D_SCRATCH_ARRAY), cff
00039 # include "grid.h"
00040 # include "ocean3d.h"
00041 # include "coupling.h"
00042 # include "mixing.h"
00043 # include "scalars.h"
00044 !
00045 #ifdef MPI
00046 #define LOCALLM Lmmpi
00047 #define LOCALMM Mmmpi
00048 #else
00049 #define LOCALLM Lm
00050 #define LOCALMM Mm
00051 #endif
00052 # include "compute_auxiliary_bounds.h"
00053 !
00054       indx=3-nstp
00055 # define nnew illegal
00056 
00057 # ifdef MIX_S_UV
00058 #  ifdef UV_VIS2
00059 !
00060 ! Compute horizontal Laplacian viscosity along constant S-surfaces.
00061 !------------------------------------------------------------------
00062 ! Compute horizontal viscous V-component momentum fluxes in XI- and
00063 ! ETA-directions, then apply them to compute harmonic viscous term.
00064 !
00065       do k=1,N
00066         do j=JstrV,Jend
00067           do i=Istr,Iend+1
00068             VFx(i,j)=
00069 #   ifdef SMAGORINSKY
00070      &               (visc2_p(i,j)+visc3d_p(i,j,k))
00071 #   else
00072      &                visc2_p(i,j)
00073 #   endif
00074      &                       *pmon_p(i,j)*0.25*( Hz(i,j,k)
00075      &                  +Hz(i-1,j,k)+Hz(i,j-1,k)+Hz(i-1,j-1,k))
00076      &                         *(v(i,j,k,nrhs)-v(i-1,j,k,nrhs))
00077 #   ifdef MASKING
00078      &                                              *pmask(i,j)
00079 #   endif
00080           enddo
00081         enddo
00082         do j=JstrV-1,Jend
00083           do i=Istr,Iend
00084             VFe(i,j)=
00085 #   ifdef SMAGORINSKY
00086      &               (visc2_r(i,j)+visc3d_r(i,j,k))
00087 #   else
00088      &                visc2_r(i,j)
00089 #   endif
00090      &                       *pnom_r(i,j)*Hz(i,j,k)
00091      &                         *(v(i,j+1,k,nrhs)-v(i,j,k,nrhs))
00092           enddo
00093         enddo
00094 
00095         do j=JstrV,Jend
00096           do i=Istr,Iend
00097             cff=0.25*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
00098      &             *(VFx(i+1,j)-VFx(i,j)+VFe(i,j)-VFe(i,j-1))
00099             rvfrc(i,j)=rvfrc(i,j)+cff
00100             v(i,j,k,indx)=v(i,j,k,indx)+cff*dt
00101           enddo
00102         enddo
00103       enddo
00104 #  endif /* UV_VIS2 */
00105 #  ifdef UV_VIS4
00106 !
00107 ! Compute horizontal biharmonic viscosity along constant S-surfaces.
00108 !-------------------------------------------------------------------
00109 ! The biharmonic operator is computed by applying the Laplacian
00110 ! operator twice.
00111 !
00112 ! Compute horizontal v-momentum flux [m^2/s] in the XI-direction.
00113 !
00114       do k=1,N
00115         do j=JstrV-1,Jend+1
00116           do i=max(Istr-1,1),min(Iend+2,LOCALLM+1)
00117             VFx(i,j)=pmon_p(i,j)*0.25*( Hz(i,j  ,k)+Hz(i-1,j  ,k)
00118      &                                 +Hz(i,j-1,k)+Hz(i-1,j-1,k))
00119      &                            *(v(i,j,k,nrhs)-v(i-1,j,k,nrhs))
00120 #   ifdef MASKING
00121      &                                                 *pmask(i,j)
00122 #   endif
00123           enddo
00124         enddo
00125 #   ifdef EW_PERIODIC
00126         if (WESTERN_EDGE) then
00127           do j=JstrV-1,Jend+1
00128             VFx(0,j)=pmon_p(Lm2,j)*(v(Lm2,j,k,nrhs)-v(L-3,j,k,nrhs))
00129      &                           *0.25*(Hz(Lm2,j  ,k)+Hz(L-3,j  ,k)+
00130      &                                  Hz(Lm2,j-1,k)+Hz(L-3,j-1,k))
00131 #    ifdef MASKING
00132             VFx(0,j)=VFx(0,j)*pmask(Lm2,j)
00133 #    endif /* MASKING */
00134           enddo
00135         endif
00136         if (EASTERN_EDGE) then
00137           do j=JstrV-1,Jend+1
00138             VFx(Lp,j)=pmon_p(3,j)*(v(3,j,k,nrhs)-v(2,j,k,nrhs))
00139      &                               *0.25*(Hz(3,j  ,k)+Hz(2,j  ,k)+
00140      &                                      Hz(3,j-1,k)+Hz(2,j-1,k))
00141 #    ifdef MASKING
00142             VFx(Lp,j)=VFx(Lp,j)*pmask(3,j)
00143 #    endif /* MASKING */
00144           enddo
00145         endif
00146 #   endif /* EW_PERIODIC */
00147 !
00148 !  Compute horizontal v-momentum flux [m^2/s] in the ETA-direction.
00149 !
00150         do j=max(JstrV-2,1),min(Jend+1,LOCALMM)
00151           do i=Istr-1,Iend+1
00152             VFe(i,j)=pnom_r(i,j)*Hz(i,j,k)*( v(i,j+1,k,nrhs)-
00153      &                                          v(i,j,k,nrhs))
00154           enddo
00155         enddo
00156 #   ifdef NS_PERIODIC
00157         if (SOUTHERN_EDGE) then
00158           do i=Istr-1,Iend+1
00159             VFe(i,0)=pnom_r(i,Mm2)*Hz(i,Mm2,k)*( v(i,Mm,k,nrhs)-
00160      &                                           v(i,Mm2,k,nrhs))
00161           enddo
00162         endif
00163         if (NORTHERN_EDGE) then
00164           do i=Istr-1,Iend+1
00165             VFe(i,M)=pnom_r(i,2)*Hz(i,2,k)*( v(i,3,k,nrhs)-
00166      &                                          v(i,2,k,nrhs))
00167           enddo
00168         endif
00169 #   endif /* NS_PERIODIC */
00170 !
00171 !  Compute first Laplacian, without mixing coefficient.
00172 !  Multiply by the metrics of the second Laplacian, [1/m/s].
00173 !  Save the result into work array "wrk".
00174 !
00175         do j=JstrV-1,Jend+1
00176           do i=Istr-1,Iend+1
00177             wrk(i,j)=(VFx(i+1,j)-VFx(i,j)+VFe(i,j)-VFe(i,j-1))*
00178      &               0.5*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))/
00179      &               (Hz(i,j,k)+Hz(i,j-1,k))
00180           enddo
00181         enddo
00182 !
00183 !  In cases other than periodicity, apply gradient or closed [free-
00184 !  slip (gamma2=1) or no-slip (gamma2=-1)] boundary conditions.
00185 !
00186 #   ifndef EW_PERIODIC
00187         if (WESTERN_EDGE) then
00188           do j=JstrV,Jend
00189 #    ifdef WESTERN_WALL
00190             wrk(Istr-1,j)=gamma2*wrk(Istr,j)
00191 #    else
00192             wrk(Istr-1,j)=0.
00193 #    endif
00194           enddo
00195         endif
00196         if (EASTERN_EDGE) then
00197           do j=JstrV,Jend
00198 #    ifdef EASTERN_WALL
00199             wrk(Iend+1,j)=gamma2*wrk(Iend,j)
00200 #    else
00201             wrk(Iend+1,j)=0.
00202 #    endif
00203           enddo
00204         endif
00205 #   endif /* !EW_PERIODIC */
00206 #   ifndef NS_PERIODIC
00207         if (SOUTHERN_EDGE) then
00208           do i=Istr-1,Iend+1
00209 #    ifdef SOUTHERN_WALL
00210             wrk(i,JstrV-1)=0.
00211 #    else
00212             wrk(i,JstrV-1)=wrk(i,JstrV)
00213 #    endif
00214           enddo
00215         endif
00216         if (NORTHERN_EDGE) then
00217           do i=Istr-1,Iend+1
00218 #    ifdef NORTHERN_WALL
00219             wrk(i,Jend+1)=0.
00220 #    else
00221             wrk(i,Jend+1)=wrk(i,Jend)
00222 #    endif
00223           enddo
00224         endif
00225 #   endif /* !NS_PERIODIC */
00226 !
00227 !  Compute d(LapV)/d(xi) and d(LapV)/d(eta) terms, then
00228 !  compute and add in biharmonic mixing [m^4/s^2] to r.h.s.
00229 !
00230         do j=JstrV,Jend
00231           do i=Istr,Iend+1
00232             VFx(i,j)=visc4_p(i,j)*pmon_p(i,j)
00233      &                 *0.25*(Hz(i,j  ,k)+Hz(i-1,j  ,k)+
00234      &                        Hz(i,j-1,k)+Hz(i-1,j-1,k))
00235      &                      *(wrk(i,j)-wrk(i-1,j))
00236 #   ifdef MASKING
00237             VFx(i,j)=VFx(i,j)*pmask(i,j)
00238 #   endif /* MASKING */
00239           enddo
00240         enddo
00241         do j=JstrV-1,Jend
00242           do i=Istr,Iend
00243             VFe(i,j)=visc4_r(i,j)*pnom_r(i,j)*Hz(i,j,k)
00244      &                           *(wrk(i,j+1)-wrk(i,j))
00245           enddo
00246         enddo
00247         cff=dt*0.25D0
00248         do j=JstrV,Jend
00249           do i=Istr,Iend
00250             v(i,j,k,indx)=v(i,j,k,indx)-cff*(pm(i,j)+pm(i,j-1))
00251      &                                     *(pn(i,j)+pn(i,j-1))
00252      &                                *( VFx(i+1,j)-VFx(i,j  )
00253      &                                  +VFe(i  ,j)-VFe(i,j-1))
00254           enddo
00255         enddo
00256       enddo
00257 #  endif /* UV_VIS4 */
00258 # endif /* MIX_S_UV */
00259 # ifdef MIX_GP_UV
00260 #  ifdef UV_VIS2
00261 !
00262 !  Add in horizontal Laplacian viscosity along rotated surfaces.
00263 !-----------------------------------------------------------------
00264 !  Compute d(v)/d(z) at horizontal V-points and vertical W-points.
00265 !
00266       do j=JstrV-1,Jend+1
00267         do k=1,N-1
00268           do i=Istr-1,Iend+1
00269             VFs(i,j,k)=(v(i,j,k+1,nrhs)-v(i,j,k,nrhs))/
00270      &                    (0.5*(z_r(i,j-1,k+1)-z_r(i,j-1,k)+
00271      &                          z_r(i,j  ,k+1)-z_r(i,j  ,k)))
00272           enddo
00273         enddo
00274         do i=Istr-1,Iend+1
00275           VFs(i,j,0)=0.
00276           VFs(i,j,N)=0.
00277         enddo
00278       enddo
00279 !-
00280       do k=N,1,-1
00281 !
00282 !  Compute contribution of the to XI-component viscous V-fluxes.
00283 !-------------------------------------------------------------------
00284 !  Average d(v)/d(z) to vertical RHO-points [1/s]. Notice that
00285 !  vertical K-loop is reversed to allow recursive use of 3D work
00286 !  array "VFs".
00287 !
00288         do j=JstrV-1,Jend+1
00289           do i=Istr-1,Iend+1
00290             dVdz(i,j)=0.5*(VFs(i,j,k)+VFs(i,j,k-1))
00291           enddo
00292         enddo
00293 !
00294 !  Compute slope (nondimensional) due mixing along geopotential
00295 !  surfaces; if applicable, add in slope due to mixing along
00296 !  isopycnal surfaces and apply land/sea mask.
00297 !
00298         do j=JstrV-1,Jend
00299           do i=Istr,Iend+1
00300             Slope(i,j)=0.5*(pm(i,j)+pm(i-1,j))
00301      &                    *(z_r(i,j,k)-z_r(i-1,j,k))
00302 #   ifdef MIX_EN_UV
00303             Slope(i,j)=Slope(i,j)+rhosx(i,j,k)
00304 #   endif
00305 #   ifdef MASKING
00306             Slope(i,j)=Slope(i,j)*umask(i,j)
00307 #   endif
00308           enddo
00309         enddo
00310 !
00311 !  Compute flux along along S-surfaces and subtract slope due to
00312 !  rotated surfaces, [m^2/s^2]. Save XI-component of the vertical
00313 !  flux into work array "wrk".
00314 !
00315         do j=JstrV,Jend
00316           do i=Istr,Iend+1
00317             cff=0.25*(pm(i,j)+pm(i,j-1)+pm(i-1,j)+pm(i-1,j-1))
00318      &                         *(v(i,j,k,nrhs)-v(i-1,j,k,nrhs))
00319 #   ifdef MASKING
00320      &                                               *pmask(i,j)
00321 #   endif
00322             cff=
00323 #   ifdef SMAGORINSKY
00324      &          (visc2_p(i,j)+visc3d_p(i,j,k))
00325 #   else
00326      &           visc2_p(i,j)
00327 #   endif
00328      &                     *(cff-0.25*(Slope(i,j)+Slope(i,j-1))
00329      &                                 *(dVdz(i,j)+dVdz(i-1,j)))
00330             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j-1))
00331             VFx(i,j)=cff*(Hz(i,j  ,k)+Hz(i-1,j  ,k)+
00332      &                    Hz(i,j-1,k)+Hz(i-1,j-1,k))/
00333      &                   (pn(i  ,j)+pn(i  ,j-1)+
00334      &                    pn(i-1,j)+pn(i-1,j-1))
00335           enddo
00336 !
00337 !  Add in divergence contribution of the XI-component V-fluxes.
00338 !  Also average XI-component of the vertical flux to horizontal
00339 !  V-points and load it back into work array "VFs".
00340 !
00341           cff=0.25*dt
00342           do i=Istr,Iend
00343             v(i,j,k,indx)=v(i,j,k,indx)+cff*(pm(i,j)+pm(i,j-1))
00344      &                                     *(pn(i,j)+pn(i,j-1))
00345      &                                   *(VFx(i+1,j)-VFx(i,j))
00346 
00347             VFs(i,j,k)=0.5*(wrk(i,j)+wrk(i+1,j))
00348           enddo
00349         enddo
00350 !
00351 !  Compute contribution of the to ETA-component viscous V-fluxes.
00352 !----------------------------------------------------------------
00353 !  Compute slope (nondimensional) due mixing along geopotential
00354 !  surfaces, if so prescribed, add in slope due to mixing along
00355 !  isopycnal surfaces and apply land/sea mask.
00356 !
00357         do j=JstrV-1,Jend+1
00358           do i=Istr,Iend
00359             Slope(i,j)=0.5*(pn(i,j)+pn(i,j-1))*
00360      &                    (z_r(i,j,k)-z_r(i,j-1,k))
00361 #   ifdef MIX_EN_UV
00362             Slope(i,j)=Slope(i,j)+rhose(i,j,k)
00363 #   endif
00364 #   ifdef MASKING
00365             Slope(i,j)=Slope(i,j)*vmask(i,j)
00366 #   endif
00367           enddo
00368         enddo
00369 !
00370 !  Compute flux along along S-surfaces and subtract slope due
00371 !  to rotated surfaces, [m^2/s^2].
00372 !
00373         do j=JstrV-1,Jend
00374           do i=Istr,Iend
00375             cff=pn(i,j)*(v(i,j+1,k,nrhs)-v(i,j,k,nrhs))
00376             cff=
00377 #   ifdef SMAGORINSKY
00378      &          (visc2_r(i,j)+visc3d_r(i,j,k))
00379 #   else
00380      &           visc2_r(i,j)
00381 #   endif
00382      &                     *(cff-0.25*(Slope(i,j)+Slope(i,j+1))
00383      &                                 *(dVdz(i,j)+dVdz(i,j+1)))
00384             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j+1))
00385             VFe(i,j)=cff*Hz(i,j,k)/pm(i,j)
00386           enddo
00387         enddo
00388 !
00389 !  Add in divergence contribution of the ETA-component V-fluxes.
00390 !  Also average ETA-component of the vertical flux to horizontal
00391 !  V-points and add it into work array "VFs".
00392 !
00393         cff=dt*0.25D0
00394         do j=JstrV,Jend
00395           do i=Istr,Iend
00396             v(i,j,k,indx)=v(i,j,k,indx)+cff*(pm(i,j)+pm(i,j-1))
00397      &                                     *(pn(i,j)+pn(i,j-1))
00398      &                                   *(VFe(i,j)-VFe(i,j-1))
00399 
00400             VFs(i,j,k)=VFs(i,j,k)+0.5*(wrk(i,j-1)+wrk(i,j))
00401           enddo
00402         enddo
00403       enddo
00404 !
00405 !  Average vertical flux to vertical W-points.
00406 !  Apply no flux boundary conditions.
00407 !
00408       do j=JstrV,Jend
00409         do k=1,N-1
00410           do i=Istr,Iend
00411             VFz(i,k)=0.5*(VFs(i,j,k)+VFs(i,j,k+1))
00412           enddo
00413         enddo
00414         do i=Istr,Iend
00415           VFz(i,0)=0.
00416           VFz(i,N)=0.
00417         enddo
00418 !
00419 !  Add in the vertical divergence of the flux [m^4/s^2] due
00420 !  to sloping S-surfaces.
00421 !
00422         do k=1,N
00423           do i=Istr,Iend
00424             v(i,j,k,indx)=v(i,j,k,indx)-dt*(VFz(i,k)-VFz(i,k-1))
00425           enddo
00426         enddo
00427       enddo
00428 #  endif /* UV_VIS2 */
00429 #  ifdef UV_VIS4
00430 !
00431 !  Add in horizontal biharmonic viscosity along rotated surfaces.
00432 !  The biharmonic operator is computed by applying the rotated
00433 !  Laplacian operator twice.
00434 !--------------------------------------------------------------------
00435 !  Evaluate the first rotated Laplacian operator.
00436 !--------------------------------------------------------------------
00437 !  Compute d(v)/d(z) at horizontal V-points and vertical W-points.
00438 !
00439       do k=1,N-1
00440         do j=max(JstrV-2,1),min(Jend+2,LOCALMM+1)
00441           do i=max(Istr-2,0),min(Iend+2,LOCALLM+1)
00442             VFs(i,j,k)=(v(i,j,k+1,nrhs)-v(i,j,k,nrhs))/
00443      &                 (0.5*(z_r(i,j-1,k+1)-z_r(i,j-1,k)+
00444      &                      z_r(i,j  ,k+1)-z_r(i,j  ,k)))
00445           enddo
00446         enddo
00447 #   ifdef EW_PERIODIC
00448         if (WESTERN_EDGE) then
00449           do j=JstrV-1,Jend+1
00450             VFs(-1,j,k)=(v(L-3,j,k+1,nrhs)-v(L-3,j,k,nrhs))/
00451      &                  (0.5*(z_r(L-3,j-1,k+1)-z_r(L-3,j-1,k)+
00452      &                       z_r(L-3,j  ,k+1)-z_r(L-3,j  ,k)))
00453           enddo
00454         endif
00455         if (EASTERN_EDGE) then
00456           do j=JstrV-1,Jend+1
00457             VFs(Lp,j,k)=(v(3,j,k+1,nrhs)-v(3,j,k,nrhs))/
00458      &                  (0.5*(z_r(3,j-1,k+1)-z_r(3,j-1,k)+
00459      &                       z_r(3,j  ,k+1)-z_r(3,j  ,k)))
00460           enddo
00461         endif
00462 #   endif /* EW_PERIODIC */
00463 #   ifdef NS_PERIODIC
00464         if (SOUTHERN_EDGE) then
00465           do i=Istr-1,Iend+1
00466             VFs(i,0,k)=(v(i,Mm2,k+1,nrhs)-v(i,Mm2,k,nrhs))/
00467      &                     (0.5*(z_r(i,M-3,k+1)-z_r(i,M-3,k)+
00468      &                          z_r(i,Mm2,k+1)-z_r(i,Mm2,k)))
00469           enddo
00470         endif
00471         if (NORTHERN_EDGE) then
00472           do i=Istr-1,Iend+1
00473             VFs(i,Mp,k)=(v(i,3,k+1,nrhs)-v(i,3,k,nrhs))/
00474      &                      (0.5*(z_r(i,2,k+1)-z_r(i,2,k)+
00475      &                           z_r(i,3,k+1)-z_r(i,3,k)))
00476           enddo
00477         endif
00478 #   endif /* NS_PERIODIC */
00479         do j=JstrV-2,Jend+2
00480           do i=Istr-2,Iend+2
00481             VFs(i,j,0)=0.
00482             VFs(i,j,N)=0.
00483           enddo
00484         enddo
00485       enddo
00486       do k=N,1,-1
00487 !
00488 !  Average d(v)/d(z) to vertical RHO-points [1/s]. Notice that the
00489 !  vertical k-loop is reversed to allow recursive use of 3D work
00490 !  array "VFs".
00491 !
00492         do j=JstrV-2,Jend+2
00493           do i=Istr-2,Iend+2
00494             dVdz(i,j)=0.5*(VFs(i,j,k)+VFs(i,j,k-1))
00495           enddo
00496         enddo
00497 !
00498 !  First Laplacian: compute lateral mixing contribution to horizontal
00499 !                   flux in the XI-direction.
00500 !--------------------------------------------------------------------
00501 !  Compute XI-component of the slope of the S-coordinate isosurfaces  
00502 !  with respect to geopotential surfaces; if applicable, add in slope 
00503 !  due to mixing along isopycnal surfaces and apply land/sea mask.   
00504 !
00505         do j=JstrV-2,Jend+1
00506           do i=max(Istr-1,1),min(Iend+2,LOCALLM+1)
00507             Slope(i,j)=0.5*(pm(i,j)+pm(i-1,j))*
00508      &                    (z_r(i,j,k)-z_r(i-1,j,k))
00509 #   ifdef MIX_EN_UV
00510             Slope(i,j)=Slope(i,j)+rhosx(i,j,k)
00511 #   endif 
00512 #   ifdef MASKING
00513             Slope(i,j)=Slope(i,j)*umask(i,j)
00514 #   endif
00515           enddo
00516         enddo
00517 #   ifdef EW_PERIODIC
00518         if (WESTERN_EDGE) then
00519           do j=JstrV-2,Jend+1
00520             Slope(0,j)=0.5*(pm(Lm2,j)+pm(L-3,j))*
00521      &                    (z_r(Lm2,j,k)-z_r(L-3,j,k))
00522 #   ifdef MIX_EN_UV
00523             Slope(0,j)=Slope(0,j)+rhosx(Lm2,j,k)
00524 #   endif
00525 #   ifdef MASKING
00526             Slope(0,j)=Slope(0,j)*umask(Lm2,j)
00527 #   endif
00528           enddo
00529         endif
00530         if (EASTERN_EDGE) then
00531           do j=JstrV-2,Jend+1
00532             Slope(Lp,j)=0.5*(pm(3,j)+pm(2,j))*
00533      &                     (z_r(3,j,k)-z_r(2,j,k))
00534 #   ifdef MIX_EN_UV
00535             Slope(Lp,j)=Slope(Lp,j)+rhosx(3,j,k)
00536 #   endif
00537 #   ifdef MASKING
00538             Slope(Lp,j)=Slope(Lp,j)*umask(3,j)
00539 #   endif /* MASKING */
00540           enddo
00541         endif
00542 #   endif /* EW_PERIODIC */
00543 !
00544 !  Compute XI-component V viscous flux (without mixing coefficient)
00545 !  along S-surfaces and subtract slope due to rotated surfaces,
00546 !  [m^2/s]. Save XI-component of the vertical flux into array "wrk".
00547 !
00548 
00549         do j=JstrV-1,Jend+1
00550           do i=max(Istr-1,1),min(Iend+2,LOCALLM+1)
00551             VFx(i,j)=0.25*(v(i,j,k,nrhs)-v(i-1,j,k,nrhs))*
00552      &                   (pm(i,j)+pm(i,j-1)+pm(i-1,j)+pm(i-1,j-1))
00553 #   ifdef MASKING
00554             VFx(i,j)=VFx(i,j)*pmask(i,j)
00555 #   endif
00556           enddo
00557         enddo
00558 #   ifdef EW_PERIODIC
00559         if (WESTERN_EDGE) then
00560           do j=JstrV-1,Jend+1
00561             VFx(0,j)=0.25*(v(Lm2,j,k,nrhs)-v(L-3,j,k,nrhs))*
00562      &               (pm(Lm2,j)+pm(Lm2,j-1)+pm(L-3,j)+pm(L-3,j-1))
00563 #   ifdef MASKING
00564             VFx(0,j)=VFx(0,j)*pmask(Lm2,j)
00565 #   endif
00566           enddo
00567         endif
00568         if (EASTERN_EDGE) then
00569           do j=JstrV-1,Jend+1
00570             VFx(Lp,j)=0.25*(v(3,j,k,nrhs)-v(2,j,k,nrhs))*
00571      &                    (pm(3,j)+pm(3,j-1)+pm(2,j)+pm(2,j-1))
00572 #   ifdef MASKING
00573             VFx(Lp,j)=VFx(Lp,j)*pmask(3,j)
00574 #   endif
00575           enddo
00576         endif
00577 #   endif
00578         do j=JstrV-1,Jend+1
00579           do i=Istr-1,Iend+2
00580             cff=VFx(i,j)-0.25*(Slope(i,j)+Slope(i,j-1))*
00581      &                          (dVdz(i,j)+dVdz(i-1,j))
00582             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j-1))
00583             VFx(i,j)=cff*(Hz(i,j  ,k)+Hz(i-1,j  ,k)+
00584      &                    Hz(i,j-1,k)+Hz(i-1,j-1,k))/
00585      &                   (pn(i  ,j)+pn(i  ,j-1)+
00586      &                    pn(i-1,j)+pn(i-1,j-1))
00587           enddo
00588 !
00589 !  Add in divergence contribution of the XI-component V-fluxes.
00590 !  Also average XI-component of the vertical flux to horizontal
00591 !  V-points and load it back into work array "VFs".
00592 !
00593           do i=Istr-1,Iend+1
00594             LapV(i,j,k)=VFx(i+1,j)-VFx(i,j)
00595             VFs(i,j,k)=0.5*(wrk(i,j)+wrk(i+1,j))
00596           enddo
00597         enddo
00598 !
00599 !  First Laplacian: compute lateral mixing contribution to horizontal
00600 !                   flux in the ETA-direction.
00601 !--------------------------------------------------------------------
00602 !  Compute ETA-component of the slope of the S-coordinate isosurfaces
00603 !  with respect to geopotential surfaces; if applicable, add in slope
00604 !  due to mixing along isopycnal surfaces and apply land/sea mask.
00605 !
00606         do j=max(JstrV-2,1),min(Jend+2,LOCALMM+1)
00607           do i=Istr-1,Iend+1
00608             Slope(i,j)=0.5*(pn(i,j)+pn(i,j-1))*
00609      &                    (z_r(i,j,k)-z_r(i,j-1,k))
00610 #   ifdef MIX_EN_UV
00611             Slope(i,j)=Slope(i,j)+rhose(i,j,k)
00612 #   endif
00613 #   ifdef MASKING
00614             Slope(i,j)=Slope(i,j)*vmask(i,j)
00615 #   endif
00616           enddo
00617         enddo
00618 #   ifdef NS_PERIODIC
00619         if (SOUTHERN_EDGE) then
00620           do i=Istr-1,Iend+1
00621             Slope(i,0)=0.5*(pn(i,Mm2)+pn(i,M-3))*
00622      &                 (z_r(i,Mm2,k)-z_r(i,M-3,k))
00623 #    ifdef MIX_EN_UV
00624             Slope(i,0)=Slope(i,0)+rhose(i,Mm2,k)
00625 #    endif
00626 #    ifdef MASKING
00627             Slope(i,0)=Slope(i,0)*vmask(i,Mm2)
00628 #    endif
00629           enddo
00630         endif
00631         if (NORTHERN_EDGE) then
00632           do i=Istr-1,Iend+1
00633             Slope(i,Mp)=0.5*(pn(i,3)+pn(i,2))*(z_r(i,3,k)-z_r(i,2,k))
00634 #    ifdef MIX_EN_UV
00635             Slope(i,Mp)=Slope(i,Mp)+rhose(i,3,k)
00636 #    endif
00637 #    ifdef MASKING
00638             Slope(i,Mp)=Slope(i,Mp)*vmask(i,3)
00639 #    endif
00640           enddo
00641         endif
00642 #   endif /* NS_PERIODIC */
00643 !
00644 !  Compute ETA-component horizontal (along S-surfaces) viscous
00645 !  V-momentum fluxes [m^4/s^2].
00646 !
00647         do j=max(JstrV-2,1),min(Jend+1,LOCALMM)
00648           do i=Istr-1,Iend+1
00649             VFe(i,j)=pn(i,j)*(v(i,j+1,k,nrhs)-v(i,j,k,nrhs))
00650           enddo
00651         enddo
00652 #   ifdef NS_PERIODIC
00653         if (SOUTHERN_EDGE) then
00654           do i=Istr-1,Iend+1
00655             VFe(i,0)=pn(i,Mm2)*(v(i,Mm,k,nrhs)-v(i,Mm2,k,nrhs))
00656           enddo
00657         endif
00658         if (NORTHERN_EDGE) then
00659           do i=Istr-1,Iend+1
00660             VFe(i,M)=pn(i,2)*(v(i,3,k,nrhs)-v(i,2,k,nrhs))
00661           enddo
00662         endif
00663 #   endif /* NS_PERIODIC */
00664         do j=JstrV-2,Jend+1
00665           do i=Istr-1,Iend+1
00666             cff=VFe(i,j)-0.5*(Slope(i,j)+Slope(i,j+1))*
00667      &                          (dVdz(i,j)+dVdz(i,j+1))
00668             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j+1))
00669             VFe(i,j)=cff*Hz(i,j,k)/pm(i,j)
00670           enddo
00671         enddo
00672 !
00673 !  Add in divergence contribution of the ETA-component V-fluxes.
00674 !  Also average ETA-component of the vertical flux to horizontal
00675 !  V-points and add it in work array "VFs".
00676 !
00677         do j=JstrV-1,Jend+1
00678           do i=Istr-1,Iend+1
00679             LapV(i,j,k)=LapV(i,j,k)+VFe(i,j)-VFe(i,j-1)
00680             VFs(i,j,k)=VFs(i,j,k)+0.5*(wrk(i,j-1)+wrk(i,j))
00681           enddo
00682         enddo
00683       enddo
00684 !
00685 !  Average vertical flux to vertical W-points.
00686 !  Apply no flux top/bottom boundary conditions.
00687 !
00688       do j=JstrV-1,Jend+1
00689         do k=1,N-1
00690           do i=Istr-1,Iend+1
00691             VFz(i,k)=0.5*(VFs(i,j,k)+VFs(i,j,k+1))
00692           enddo
00693         enddo
00694         do i=Istr-1,Iend+1
00695           VFz(i,0)=0.
00696           VFz(i,N)=0.
00697         enddo
00698 !
00699 !  Add in the vertical divergence of the flux due to sloping
00700 !  S-surfaces in the ETA-direction. Multiply the first Laplacian
00701 !  by the metrics of the second Laplacian, [1/m/s].
00702 !
00703         do k=1,N
00704           do i=Istr-1,Iend+1
00705             LapV(i,j,k)=( LapV(i,j,k)*0.25*(pn(i,j)+pn(i,j-1))*
00706      &                    (pm(i,j)+pm(i,j-1))-VFz(i,k)+VFz(i,k-1)
00707      &                            )/(0.5*(Hz(i,j,k)+Hz(i,j-1,k)))
00708           enddo
00709         enddo
00710       enddo
00711 !
00712 !   In cases other than periodicity, apply gradient or closed
00713 !  [free-slip (gamma2=1) or no-slip (gamma2=-1)] boundary conditions.
00714 !
00715 #   ifndef EW_PERIODIC
00716       if (WESTERN_EDGE) then
00717         do k=1,N
00718           do j=JstrV-1,Jend+1
00719 #    ifdef WESTERN_WALL
00720             LapV(Istr-1,j,k)=gamma2*LapV(Istr,j,k)
00721 #    else
00722             LapV(Istr-1,j,k)=0.
00723 #    endif
00724           enddo
00725         enddo
00726       endif
00727       if (EASTERN_EDGE) then
00728         do k=1,N
00729           do j=JstrV-1,Jend+1
00730 #    ifdef EASTERN_WALL
00731             LapV(Iend+1,j,k)=gamma2*LapV(Iend,j,k)
00732 #    else
00733             LapV(Iend+1,j,k)=0.
00734 #    endif
00735           enddo
00736         enddo
00737       endif
00738 #   endif /* !EW_PERIODIC */
00739 #   ifndef NS_PERIODIC
00740       if (SOUTHERN_EDGE) then
00741         do k=1,N
00742           do i=Istr-1,Iend+1
00743 #    ifdef SOUTHERN_WALL
00744             LapV(i,JstrV-1,k)=0.
00745 #    else
00746             LapV(i,JstrV-1,k)=LapV(i,JstrV,k)
00747 #    endif
00748           enddo
00749         enddo
00750       endif
00751       if (NORTHERN_EDGE) then
00752         do k=1,N
00753           do i=Istr-1,Iend+1
00754 #    ifdef NORTHERN_WALL
00755             LapV(i,Jend+1,k)=0.
00756 #    else
00757             LapV(i,Jend+1,k)=LapV(i,Jend,k)
00758 #    endif
00759           enddo
00760         enddo
00761       endif
00762 #   endif /* !NS_PERIODIC */
00763 !
00764 !  Evaluate the second rotated Laplacian operator.
00765 !--------------------------------------------------------------------
00766 !  Compute VFs=d(LAP(v))/d(z) at horizontal V-points and vertical
00767 !  W-points; set verical boundary conditions for it.
00768 !
00769       do j=JstrV-1,Jend+1
00770         do k=1,N-1
00771           do i=Istr-1,Iend+1
00772             VFs(i,j,k)=(LapV(i,j,k+1)-LapV(i,j,k))/
00773      &                 (0.5*((z_r(i,j-1,k+1)-z_r(i,j-1,k))+
00774      &                      (z_r(i,j  ,k+1)-z_r(i,j  ,k))))
00775      &                     
00776           enddo
00777         enddo
00778         do i=Istr-1,Iend+1
00779           VFs(i,j,0)=0.
00780           VFs(i,j,N)=0.
00781         enddo
00782       enddo
00783       do k=N,1,-1
00784 !
00785 !  Compute contribution of the lateral XI-component horizontal
00786 !  viscous LapU-fluxes.
00787 !--------------------------------------------------------------------
00788 !  Average d(LapV)/d(z) to vertical RHO-points [1/m^2/s]. Notice
00789 !  that the vertical k-loop is reversed to allow recursive reuse
00790 !  of 3D work array "VFs".
00791 !
00792         do j=JstrV-1,Jend+1
00793           do i=Istr-1,Iend+1
00794             dVdz(i,j)=0.5*(VFs(i,j,k)+VFs(i,j,k-1))
00795           enddo
00796         enddo
00797 !
00798 !  Compute slope (nondimensional) due mixing along geopotential
00799 !  surfaces; if applicable, add in slope (nondimensional) due
00800 !  to mixing along isopycnal surfaces and apply land/sea mask.
00801 !
00802         do j=JstrV-1,Jend
00803           do i=Istr,Iend+1
00804             Slope(i,j)=0.5*(pm(i,j)+pm(i-1,j))*
00805      &                    (z_r(i,j,k)-z_r(i-1,j,k))
00806 #   ifdef MIX_EN_UV
00807             Slope(i,j)=Slope(i,j)+rhosx(i,j,k)
00808 #   endif
00809 #   ifdef MASKING
00810             Slope(i,j)=Slope(i,j)*umask(i,j)
00811 #   endif
00812           enddo
00813         enddo
00814 !
00815 !  Compute XI-component viscous LapV-fluxes along S-surfaces
00816 !  and subtract the slope due to rotated surfaces, [m^2/s^2].
00817 !  Save XI-component of the vertical flux in work array "wrk".
00818 !
00819         do j=JstrV,Jend
00820           do i=Istr,Iend+1
00821             cff=0.25*(LapV(i,j,k)-LapV(i-1,j,k))*
00822      &              (pm(i,j)+pm(i,j-1)+pm(i-1,j)+pm(i-1,j-1))
00823 #   ifdef MASKING
00824             cff=cff*pmask(i,j)
00825 #   endif
00826             cff=-visc4_p(i,j)*(cff-0.25*(Slope(i,j)+Slope(i,j-1))
00827      &                                   *(dVdz(i,j)+dVdz(i-1,j)))
00828             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j-1))
00829             VFx(i,j)=cff*(Hz(i,j  ,k)+Hz(i-1,j  ,k)+
00830      &                    Hz(i,j-1,k)+Hz(i-1,j-1,k))/
00831      &                   (pn(i,j)+pn(i,j-1)+pn(i-1,j)+pn(i-1,j-1))
00832           enddo
00833 !
00834 !  Add in divergence contribution of the XI-component LapV-fluxes.
00835 !  Also average XI-component of the vertical flux to horizontal
00836 !  V-points and load it back into work array "VFs".
00837 !
00838           cff=dt*0.25D0
00839           do i=Istr,Iend
00840 c*          rv(i,j,k,nrhs)=rv(i,j,k,nrhs)+VFx(i+1,j)-VFx(i,j)
00841 
00842             v(i,j,k,indx)=v(i,j,k,indx)+cff*(pm(i,j)+pm(i,j-1))
00843      &                                     *(pn(i,j)+pn(i,j-1))
00844      &                                   *(VFx(i+1,j)-VFx(i,j))
00845 
00846             VFs(i,j,k)=0.5*(wrk(i,j)+wrk(i+1,j))
00847           enddo
00848         enddo
00849 !
00850 !  Compute contribution of the lateral ETA-component horizontal
00851 !  viscous LapV-fluxes.
00852 !--------------------------------------------------------------------
00853 !  Compute slope (nondimensional) due mixing along geopotential
00854 !  surfaces; if applicable, add in slope due to mixing along
00855 !  isopycnal surfaces; apply land/sea mask.
00856 !
00857         do j=JstrV-1,Jend+1
00858           do i=Istr,Iend
00859             Slope(i,j)=0.5*(pn(i,j)+pn(i,j-1))*
00860      &                    (z_r(i,j,k)-z_r(i,j-1,k))
00861 #   ifdef MIX_EN_UV
00862             Slope(i,j)=Slope(i,j)+rhose(i,j,k)
00863 #   endif
00864 #   ifdef MASKING
00865             Slope(i,j)=Slope(i,j)*vmask(i,j)
00866 #   endif
00867           enddo
00868         enddo
00869 !
00870 !  Compute ETA-component viscous LapV-fluxes along along S-surfaces
00871 !  and subtract the slope due to rotated surfaces, [m^2/s^2]; Save
00872 !  ETA-component of the vertical flux in work array "wrk".
00873 !
00874         do j=JstrV-1,Jend
00875           do i=Istr,Iend
00876             cff=pn(i,j)*(LapV(i,j+1,k)-LapV(i,j,k))
00877             cff=-visc4_r(i,j)*(cff-0.25*(Slope(i,j)+Slope(i,j+1))*
00878      &                             (dVdz(i,j)+dVdz(i,j+1)))
00879             wrk(i,j)=cff*0.5*(Slope(i,j)+Slope(i,j+1))
00880             VFe(i,j)=cff*Hz(i,j,k)/pm(i,j)
00881           enddo
00882         enddo
00883 !
00884 !  Add in divergence contribution of the ETA-component LapV-fluxes.
00885 !  Also average ETA-component of the vertical flux to horizontal
00886 !  V-points and add it in work array "VFs".
00887 !
00888         cff=dt*0.25D0
00889         do j=JstrV,Jend
00890           do i=Istr,Iend
00891             v(i,j,k,indx)=v(i,j,k,indx)+cff*(pm(i,j)+pm(i,j-1))
00892      &                                     *(pn(i,j)+pn(i,j-1))
00893      &                                   *(VFe(i,j)-VFe(i,j-1))
00894 
00895             VFs(i,j,k)=VFs(i,j,k)+0.5*(wrk(i,j-1)+wrk(i,j))
00896           enddo
00897         enddo
00898       enddo
00899 !
00900 !  Average vertical flux to vertical W-points.
00901 !  Apply no flux boundary conditions.
00902 !
00903       do j=JstrV,Jend
00904         do k=1,N-1
00905           do i=Istr,Iend
00906             VFz(i,k)=0.5*(VFs(i,j,k)+VFs(i,j,k+1))
00907           enddo
00908         enddo
00909         do i=Istr,Iend
00910           VFz(i,0)=0.
00911           VFz(i,N)=0.
00912         enddo
00913 !
00914 !  Add in divergence contribution of vertical component fluxes
00915 !  due to the sloping S-surfaces.
00916 !
00917         do k=1,N
00918           do i=Istr,Iend
00919             v(i,j,k,indx)=v(i,j,k,indx)-dt*(VFz(i,k)-VFz(i,k-1))
00920           enddo
00921         enddo
00922       enddo
00923 #  endif /* UV_VIS4 */
00924 # endif /* MIX_GP_UV */
00925 #else
00926       subroutine v3dmix_empty
00927 #endif /* SOLVE3D */
00928       return
00929       end
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/JUNK/v3dmix.F
