ROMS_r215: /home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/uv3dmix

00001 !
00002 ! $Id: uv3dmix_S.F,v 1.3 2005/10/10 13:40:19 pmarches Exp $
00003 !
00004 #include "cppdefs.h"
00005 # if defined SOLVE3D && defined SPONGE_VIS2 && !defined UV_VIS2
00006 
00007 # define MIX_S_UV_SPG
00008 # if defined M3CLIMATOLOGY && !defined AGRIF
00009 #  define CLIMAT_SPONGE
00010 # endif
00011 
00012 # ifdef MIX_S_UV_SPG
00013 !
00014 !====================================================================
00015 !
00016 ! Compute horizontal harmonic viscosity along S-surfaces
00017 ! in SPONGE layers
00018 !
00019 !====================================================================
00020 !
00021       subroutine uv3dmix_spg (tile)
00022       implicit none
00023       integer tile, trd, omp_get_thread_num
00024 # include "param.h"
00025 # include "private_scratch.h"
00026 # include "compute_tile_bounds.h"
00027       trd=omp_get_thread_num()
00028       call uv3dmix_spg_tile (Istr,Iend,Jstr,Jend, 
00029      &                        A2d(1,1,trd), A2d(1,2,trd),
00030      &                        A2d(1,3,trd), A2d(1,4,trd))
00031       return
00032       end
00033  
00034       subroutine uv3dmix_spg_tile (Istr,Iend,Jstr,Jend, UFx,UFe,VFx,VFe)
00035 !
00036 ! Computes harmonic mixing of momentum as horizontal divergence of the 
00037 ! stress tensor. Components of the stress tensor are:
00038 !                         du     dv
00039 !         s_xx = -s_yy = ---- - -----
00040 !                         dx     dy
00041 !
00042 !                         du     dv
00043 !         s_xy =  s_yx = ---- + ----
00044 !                         dy     dx
00045 !
00046 ! Reference:
00047 !
00048 ! Wajsowicz, R.C, 1993: A consistent formulation of the anisotropic
00049 !     stress tensor for use in models of the large-scale ocean
00050 !     circulation, JCP, 105, 333-338.
00051 !
00052 ! Sadourny, R. and K. Maynard, 1997: Formulations of lateral
00053 !     diffusion in geophysical fluid dynamics models, In "Numerical
00054 !     Methods of Atmospheric and Oceanic Modelling". Lin, Laprise,
00055 !     and Ritchie, Eds., NRC Research Press, 547-556.
00056 !
00057 ! Griffies, S.M. and R.W. Hallberg, 2000: Biharmonic friction with
00058 !     a Smagorinsky-like viscosity for use in large-scale eddy-
00059 !     permitting ocean models, Mon. Wea. Rev., 128, 8, 2935-2946.
00060 !
00061 
00062       implicit none
00063       integer Istr,Iend,Jstr,Jend, i,j,k, indx
00064       real UFe(PRIVATE_2D_SCRATCH_ARRAY),
00065      &     UFx(PRIVATE_2D_SCRATCH_ARRAY), cff,
00066      &     VFe(PRIVATE_2D_SCRATCH_ARRAY), cff1,
00067      &     VFx(PRIVATE_2D_SCRATCH_ARRAY)
00068 # include "param.h"
00069 # include "scalars.h"
00070 # include "grid.h"
00071 # include "ocean3d.h"
00072 # include "coupling.h"
00073 # include "mixing.h"
00074 # include "climat.h"
00075 # ifdef DIAGNOSTICS_UV
00076 #  include "diagnostics.h"
00077 # endif
00078 !
00079 # include "compute_auxiliary_bounds.h"
00080 !
00081       indx=3-nstp    !--> time index for target arrays;
00082 
00083 !
00084 !  Compute flux-components of the horizontal divergence of the stress
00085 !  tensor (m5/s2) in XI- and ETA-directions.
00086 !
00087       do k=1,N
00088         do j=JstrV-1,Jend
00089           do i=IstrU-1,Iend
00090             cff=Hz(i,j,k)*visc2_r(i,j)*
00091 #  ifdef CLIMAT_SPONGE
00092      &    ( on_r(i,j)*pm(i,j)*
00093      &        ( pn_u(i+1,j)*(u(i+1,j,k,nstp)-uclm(i+1,j,k))
00094      &         -pn_u(i  ,j)*(u(i  ,j,k,nstp)-uclm(i  ,j,k)) )
00095      &     -om_r(i,j)*pn(i,j)*
00096      &        ( pm_v(i,j+1)*(v(i,j+1,k,nstp)-vclm(i,j+1,k))
00097      &         -pm_v(i,j  )*(v(i,j  ,k,nstp)-vclm(i,j  ,k)) ) )
00098 #  else
00099      &    ( on_r(i,j)*pm(i,j)*( pn_u(i+1,j)*u(i+1,j,k,nstp)
00100      &                         -pn_u(i  ,j)*u(i  ,j,k,nstp) )
00101      &     -om_r(i,j)*pn(i,j)*( pm_v(i,j+1)*v(i,j+1,k,nstp)
00102      &                         -pm_v(i,j  )*v(i,j  ,k,nstp) ) )
00103 #  endif
00104             UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
00105             VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
00106           enddo
00107         enddo
00108         do j=Jstr,Jend+1
00109           do i=Istr,Iend+1
00110             cff=0.0625*visc2_p(i,j)*
00111      &      (Hz(i-1,j,k)+Hz(i,j,k)+Hz(i-1,j-1,k)+Hz(i,j-1,k))*
00112 #  ifdef CLIMAT_SPONGE
00113      &      (  (pm(i-1,j)+pm(i,j)+pm(i-1,j-1)+pm(i,j-1))*on_p(i,j)
00114      &              *( pn_v(i  ,j)*(v(i  ,j,k,nstp)-vclm(i  ,j,k))
00115      &                -pn_v(i-1,j)*(v(i-1,j,k,nstp)-vclm(i-1,j,k)) )
00116      &        +(pn(i-1,j)+pn(i,j)+pn(i-1,j-1)+pn(i,j-1))*om_p(i,j)
00117      &              *( pm_u(i,j  )*(u(i,j  ,k,nstp)-uclm(i  ,j,k))
00118      &                -pm_u(i,j-1)*(u(i,j-1,k,nstp)-uclm(i-1,j,k)) ) )
00119 #  else
00120      &      (  (pm(i-1,j)+pm(i,j)+pm(i-1,j-1)+pm(i,j-1))*on_p(i,j)
00121      &                     *( pn_v(i  ,j)*v(i  ,j,k,nstp)
00122      &                       -pn_v(i-1,j)*v(i-1,j,k,nstp) )
00123      &        +(pn(i-1,j)+pn(i,j)+pn(i-1,j-1)+pn(i,j-1))*om_p(i,j)
00124      &                     *( pm_u(i,j  )*u(i,j  ,k,nstp)
00125      &                       -pm_u(i,j-1)*u(i,j-1,k,nstp) ) )
00126 #  endif
00127 #  ifdef MASKING
00128      &                                              *pmask(i,j)
00129 #  endif
00130             UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
00131             VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
00132           enddo
00133         enddo
00134 !
00135 ! Apply viscous terms. Note that at this stage arrays u,v(...,3-nstp)
00136 ! contain Hz*U and Hz*V with units of [m2/s].   Also compute vertical
00137 ! integral of viscous terms and add it into coupling terms for the
00138 ! barotropic mode
00139 !
00140         do j=Jstr,Jend
00141           do i=IstrU,Iend
00142             cff=pn_u(i,j)*(UFx(i,j)-UFx(i-1,j))
00143      &         +pm_u(i,j)*(UFe(i,j+1)-UFe(i,j))
00144             cff1=pm_u(i,j)*pn_u(i,j)*cff
00145             rufrc(i,j)=rufrc(i,j) + cff
00146             u(i,j,k,indx)=u(i,j,k,indx) + dt*cff1
00147 #  ifdef DIAGNOSTICS_UV
00148             MHmix(i,j,k,1) = MHmix(i,j,k,1)+
00149      &                            2*cff1/(Hz(i-1,j,k)+Hz(i,j,k))
00150 #  endif
00151           enddo
00152         enddo
00153         do j=JstrV,Jend
00154           do i=Istr,Iend
00155             cff=pn_v(i,j)*(VFx(i+1,j)-VFx(i,j))
00156      &         -pm_v(i,j)*(VFe(i,j)-VFe(i,j-1))
00157             cff1=pm_v(i,j)*pn_v(i,j)*cff 
00158             rvfrc(i,j)=rvfrc(i,j) + cff
00159             v(i,j,k,indx)=v(i,j,k,indx) + dt*cff1
00160 #  ifdef DIAGNOSTICS_UV
00161             MHmix(i,j,k,2) = MHmix(i,j,k,2)+
00162      &                            2*cff1/(Hz(i,j,k)+Hz(i,j-1,k))
00163 #  endif
00164           enddo
00165         enddo
00166       enddo
00167       return
00168       end
00169 
00170 # else
00171 !
00172 !====================================================================
00173 !
00174 ! Compute horizontal harmonic viscosity along geopotential surfaces
00175 ! in SPONGE layers
00176 !
00177 !====================================================================
00178 !
00179       subroutine uv3dmix_spg (tile)
00180       implicit none
00181       integer tile, trd, omp_get_thread_num
00182 #include "param.h"
00183 #include "private_scratch.h"
00184 #include "compute_tile_bounds.h"
00185       trd=omp_get_thread_num() 
00186       call uv3dmix_spg_tile (Istr,Iend,Jstr,Jend,
00187      &    A2d(1, 1,trd), A2d(1, 2,trd), A2d(1, 3,trd), A2d(1,4,trd),
00188      &    A2d(1, 5,trd), A2d(1, 7,trd), A2d(1, 9,trd), A2d(1,11,trd),
00189      &    A2d(1,13,trd), A2d(1,15,trd), A2d(1,17,trd), A2d(1,19,trd),
00190      &    A2d(1,21,trd), A2d(1,23,trd), A2d(1,25,trd), A2d(1,27,trd))
00191       return
00192       end
00193  
00194       subroutine uv3dmix_spg_tile (Istr,Iend,Jstr,Jend, UFx, UFe,
00195      &                             VFx,VFe, UFs,VFs, dnUdx,  dmUde,
00196      &                             dUdz,    dnVdx,   dmVde,  dVdz,
00197      &                             dZdx_r,  dZdx_p,  dZde_r, dZde_p)
00198       implicit none
00199       integer Istr,Iend,Jstr,Jend, i,j,k, k1,k2, indx
00200       real UFe(PRIVATE_2D_SCRATCH_ARRAY),
00201      &     VFe(PRIVATE_2D_SCRATCH_ARRAY),     cff,
00202      &     UFx(PRIVATE_2D_SCRATCH_ARRAY),     cff1,
00203      &     VFx(PRIVATE_2D_SCRATCH_ARRAY),     cff2,
00204      &     UFs(PRIVATE_2D_SCRATCH_ARRAY,2),   cff3,
00205      &     VFs(PRIVATE_2D_SCRATCH_ARRAY,2),   cff4,
00206      &   dmUde(PRIVATE_2D_SCRATCH_ARRAY,2),   cff5,
00207      &   dmVde(PRIVATE_2D_SCRATCH_ARRAY,2),   cff6,
00208      &   dnUdx(PRIVATE_2D_SCRATCH_ARRAY,2),   cff7,
00209      &   dnVdx(PRIVATE_2D_SCRATCH_ARRAY,2),   cff8,
00210      &    dUdz(PRIVATE_2D_SCRATCH_ARRAY,2),
00211      &    dVdz(PRIVATE_2D_SCRATCH_ARRAY,2),   dmUdz,
00212      &  dZde_p(PRIVATE_2D_SCRATCH_ARRAY,2),   dnUdz,
00213      &  dZde_r(PRIVATE_2D_SCRATCH_ARRAY,2),   dmVdz,
00214      &  dZdx_p(PRIVATE_2D_SCRATCH_ARRAY,2),   dnVdz,
00215      &  dZdx_r(PRIVATE_2D_SCRATCH_ARRAY,2)
00216  
00217 # include "param.h"
00218 # include "scalars.h"
00219 # include "grid.h"
00220 # include "ocean3d.h"
00221 # include "coupling.h"
00222 # include "mixing.h"
00223 # ifdef DIAGNOSTICS_UV
00224 #  include "diagnostics.h"
00225 # endif
00226 !
00227 # include "compute_auxiliary_bounds.h"
00228 
00229       indx=3-nstp    !--> time index for target arrays;
00230 
00231       k2=1
00232       do k=0,N,+1  !--> irreversible
00233         k1=k2
00234         k2=3-k1
00235         if (k.lt.N) then
00236           do j=Jstr-1,Jend+1
00237             do i=IstrU-1,Iend+1
00238               UFx(i,j)=(z_r(i,j,k+1)-z_r(i-1,j,k+1))*pm_u(i,j)
00239 # ifdef MASKING
00240      &                                            *umask(i,j)
00241 # endif
00242             enddo
00243           enddo
00244           do j=JstrV-1,Jend+1
00245             do i=Istr-1,Iend+1
00246               VFe(i,j)=(z_r(i,j,k+1)-z_r(i,j-1,k+1))*pn_v(i,j)
00247 # ifdef MASKING
00248      &                                            *vmask(i,j)
00249 # endif
00250             enddo
00251           enddo
00252           do j=JstrV-1,Jend
00253             do i=IstrU-1,Iend
00254               dnUdx(i,j,k2)=pm(i,j)*
00255 #  ifdef CLIMAT_SPONGE
00256      &        ( pn_u(i+1,j)*(u(i+1,j,k+1,nstp)-uclm(i+1,j,k+1))
00257      &         -pn_u(i  ,j)*(u(i  ,j,k+1,nstp)-uclm(i  ,j,k+1)) )
00258 #  else
00259      &              (pn_u(i+1,j)*u(i+1,j,k+1,nstp)
00260      &              -pn_u(i  ,j)*u(i  ,j,k+1,nstp))
00261 #  endif
00262 # ifdef MASKING
00263      &                                            *rmask(i,j)
00264 # endif
00265               dmVde(i,j,k2)=pn(i,j)*
00266 #  ifdef CLIMAT_SPONGE
00267      &        ( pm_v(i,j+1)*(v(i,j+1,k+1,nstp)-vclm(i,j+1,k+1))
00268      &         -pm_v(i,j  )*(v(i,j  ,k+1,nstp)-vclm(i,j  ,k+1)) )
00269 #  else
00270      &              (pm_v(i,j+1)*v(i,j+1,k+1,nstp)
00271      &              -pm_v(i,j  )*v(i,j  ,k+1,nstp))
00272 #  endif
00273 # ifdef MASKING
00274      &                                            *rmask(i,j)
00275 # endif
00276               dZdx_r(i,j,k2)=0.5*(UFx(i,j)+UFx(i+1,j))
00277               dZde_r(i,j,k2)=0.5*(VFe(i,j)+VFe(i,j+1))
00278             enddo
00279           enddo
00280           do j=Jstr,Jend+1
00281             do i=Istr,Iend+1
00282               dmUde(i,j,k2)=0.25*(pn(i,j)+pn(i-1,j)+pn(i,j-1)
00283      &                                            +pn(i-1,j-1))
00284 #  ifdef CLIMAT_SPONGE
00285      &       *( pm_u(i,j  )*(u(i,j  ,k+1,nstp)-uclm(i+1,j,k+1))
00286      &         -pm_u(i,j-1)*(u(i,j-1,k+1,nstp)-uclm(i  ,j,k+1)) )
00287 #  else
00288      &             *(pm_u(i,j  )*u(i,j  ,k+1,nstp)
00289      &              -pm_u(i,j-1)*u(i,j-1,k+1,nstp))
00290 #  endif
00291 # ifdef MASKING
00292      &                                             *pmask(i,j)
00293 # endif
00294               dnVdx(i,j,k2)=0.25*(pm(i,j)+pm(i-1,j)+pm(i,j-1)
00295      &                                            +pm(i-1,j-1))
00296 #  ifdef CLIMAT_SPONGE
00297      &       *( pn_v(i  ,j)*(v(i  ,j,k+1,nstp)-vclm(i  ,j,k+1))
00298      &         -pn_v(i-1,j)*(v(i-1,j,k+1,nstp)-vclm(i-1,j,k+1)) )
00299 #  else
00300      &             *(pn_v(i  ,j)*v(i  ,j,k+1,nstp)
00301      &              -pn_v(i-1,j)*v(i-1,j,k+1,nstp))
00302 #  endif
00303 # ifdef MASKING
00304      &                                             *pmask(i,j)
00305 # endif
00306               dZde_p(i,j,k2)=0.5*(VFe(i-1,j)+VFe(i,j))
00307               dZdx_p(i,j,k2)=0.5*(UFx(i,j-1)+UFx(i,j))
00308             enddo
00309           enddo          !--> discard UFx,VFe, keep all others
00310         endif
00311  
00312         if (k.eq.0 .or. k.eq.N) then
00313           do j=Jstr-1,Jend+1
00314             do i=IstrU-1,Iend+1
00315               dUdz(i,j,k2)=0.
00316               UFs(i,j,k2)=0.
00317             enddo
00318           enddo
00319           do j=JstrV-1,Jend+1
00320             do i=Istr-1,Iend+1
00321               dVdz(i,j,k2)=0.
00322               VFs(i,j,k2)=0.
00323             enddo
00324           enddo
00325         else
00326           do j=Jstr-1,Jend+1
00327             do i=IstrU-1,Iend+1
00328 #  ifdef CLIMAT_SPONGE
00329               dUdz(i,j,k2)=2.*(u(i,j,k+1,nstp)-u(i,j,k,nstp)
00330      &                            -uclm(i,j,k+1)+uclm(i,j,k))
00331 #  else
00332               dUdz(i,j,k2)=2.*(u(i,j,k+1,nstp)-u(i,j,k,nstp))
00333 #  endif
00334      &                          /( z_r(i-1,j,k+1)-z_r(i-1,j,k)
00335      &                            +z_r(i  ,j,k+1)-z_r(i  ,j,k))
00336             enddo
00337           enddo
00338           do j=JstrV-1,Jend+1
00339             do i=Istr-1,Iend+1
00340 #  ifdef CLIMAT_SPONGE
00341               dVdz(i,j,k2)=2.*(v(i,j,k+1,nstp)-v(i,j,k,nstp)
00342      &                            -vclm(i,j,k+1)+vclm(i,j,k))
00343 #  else
00344               dVdz(i,j,k2)=2.*(v(i,j,k+1,nstp)-v(i,j,k,nstp))
00345 #  endif
00346      &                          /( z_r(i,j-1,k+1)-z_r(i,j-1,k)
00347      &                            +z_r(i,j  ,k+1)-z_r(i,j  ,k))
00348             enddo
00349           enddo
00350         endif
00351 !
00352 ! Compute components of the rotated viscous flux [m5/s2] along
00353 ! geopotential surfaces in the XI- and ETA-directions.
00354 !
00355         if (k.gt.0) then
00356           do j=JstrV-1,Jend
00357             do i=IstrU-1,Iend
00358               cff=visc2_r(i,j)*Hz(i,j,k)*(
00359      &         om_r(i,j)*( dnUdx(i,j,k1) - 0.5*pn(i,j)*(
00360      &           min(dZdx_r(i,j,k1),0.)*(dUdz(i,j,k1)+dUdz(i+1,j,k2))
00361      &          +max(dZdx_r(i,j,k1),0.)*(dUdz(i,j,k2)+dUdz(i+1,j,k1))
00362      &                                                             ))
00363      &        -om_r(i,j)*( dmVde(i,j,k1) - 0.5*pm(i,j)*(
00364      &           min(dZde_r(i,j,k1),0.)*(dVdz(i,j,k1)+dVdz(i,j+1,k2))
00365      &          +max(dZde_r(i,j,k1),0.)*(dVdz(i,j,k2)+dVdz(i,j+1,k1))
00366      &                                                            )))
00367 # ifdef MASKING
00368      &                                                   *rmask(i,j)
00369 # endif
00370               UFx(i,j)=om_r(i,j)*om_r(i,j)*cff
00371               VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
00372             enddo
00373           enddo
00374           do j=Jstr,Jend+1
00375             do i=Istr,Iend+1
00376               cff=0.25*visc2_p(i,j)
00377      &             *(Hz(i,j,k)+Hz(i-1,j,k)+Hz(i,j-1,k)+Hz(i-1,j-1,k))
00378      &       *( om_p(i,j)*( dnVdx(i,j,k1)-0.125*( pn(i,j)+pn(i-1,j)
00379      &                                        +pn(i,j-1)+pn(i-1,j-1))
00380      &        *( min(dZdx_p(i,j,k1),0.)*(dVdz(i-1,j,k1)+dVdz(i,j,k2))
00381      &          +max(dZdx_p(i,j,k1),0.)*(dVdz(i-1,j,k2)+dVdz(i,j,k1))
00382      &                                                             )) 
00383      &        + om_p(i,j)*( dmUde(i,j,k1)-0.125*( pm(i,j)+pm(i-1,j)
00384      &                                        +pm(i,j-1)+pm(i-1,j-1))
00385      &        *( min(dZde_p(i,j,k1),0.)*(dUdz(i,j-1,k1)+dUdz(i,j,k2))
00386      &          +max(dZde_p(i,j,k1),0.)*(dUdz(i,j-1,k2)+dUdz(i,j,k1))
00387      &                                                            )))
00388 # ifdef MASKING
00389      &                                                   *pmask(i,j)
00390 # endif
00391               UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
00392               VFx(i,j)=om_p(i,j)*om_p(i,j)*cff
00393             enddo
00394           enddo
00395 !
00396 ! Compute vertical flux [m^2/s^2] due to sloping terrain-following
00397 ! surfaces.
00398 !
00399           if (k.lt.N) then
00400             do j=Jstr,Jend
00401               do i=IstrU,Iend
00402                 cff1=pn_u(i,j)
00403                 cff2=pm_u(i,j)
00404                 cff=0.25*( dVdz(i,j,k2)+dVdz(i-1,j,k2)
00405      &                +dVdz(i,j+1,k2)+dVdz(i-1,j+1,k2))
00406                 dnUdz=cff1*dUdz(i,j,k2)
00407                 dmUdz=cff2*dUdz(i,j,k2)
00408                 dnVdz=cff1*cff
00409                 dmVdz=cff2*cff
00410  
00411                 cff1=min(dZdx_r(i-1,j,k1),0.)
00412                 cff2=min(dZdx_r(i  ,j,k2),0.)
00413                 cff3=max(dZdx_r(i-1,j,k2),0.)
00414                 cff4=max(dZdx_r(i  ,j,k1),0.)
00415                 cff5=min(dZde_r(i-1,j,k1),0.)
00416                 cff6=min(dZde_r(i  ,j,k2),0.)
00417                 cff7=max(dZde_r(i-1,j,k2),0.)
00418                 cff8=max(dZde_r(i  ,j,k1),0.)
00419  
00420                 cff=om_u(i,j)*( cff1*(cff1*dnUdz-dnUdx(i-1,j,k1))
00421      &                         +cff2*(cff2*dnUdz-dnUdx(i  ,j,k2))
00422      &                         +cff3*(cff3*dnUdz-dnUdx(i-1,j,k2))
00423      &                         +cff4*(cff4*dnUdz-dnUdx(i  ,j,k1))
00424      &                                                          )
00425      &             -om_u(i,j)*( cff1*(cff5*dmVdz-dmVde(i-1,j,k1))
00426      &                         +cff2*(cff6*dmVdz-dmVde(i  ,j,k2))
00427      &                         +cff3*(cff7*dmVdz-dmVde(i-1,j,k2))
00428      &                         +cff4*(cff8*dmVdz-dmVde(i  ,j,k1))
00429      &                                                          )
00430                 cff1=min(dZde_p(i,j  ,k1),0.)
00431                 cff2=min(dZde_p(i,j+1,k2),0.)
00432                 cff3=max(dZde_p(i,j  ,k2),0.)
00433                 cff4=max(dZde_p(i,j+1,k1),0.)
00434                 cff5=min(dZdx_p(i,j  ,k1),0.)
00435                 cff6=min(dZdx_p(i,j+1,k2),0.)
00436                 cff7=max(dZdx_p(i,j  ,k2),0.)
00437                 cff8=max(dZdx_p(i,j+1,k1),0.)
00438  
00439                 cff=cff + om_u(i,j)*(
00440      &                          cff1*(cff1*dmUdz-dmUde(i,j  ,k1))
00441      &                         +cff2*(cff2*dmUdz-dmUde(i,j+1,k2))
00442      &                         +cff3*(cff3*dmUdz-dmUde(i,j  ,k2))
00443      &                         +cff4*(cff4*dmUdz-dmUde(i,j+1,k1))
00444      &                                                          )
00445      &             +om_u(i,j)*( cff1*(cff5*dnVdz-dnVdx(i,j  ,k1))
00446      &                         +cff2*(cff6*dnVdz-dnVdx(i,j+1,k2))
00447      &                         +cff3*(cff7*dnVdz-dnVdx(i,j  ,k2))
00448      &                         +cff4*(cff8*dnVdz-dnVdx(i,j+1,k1))
00449      &                                                          )
00450                 UFs(i,j,k2)=0.25*(visc2_r(i-1,j)+visc2_r(i,j))*cff
00451               enddo
00452             enddo
00453  
00454             do j=JstrV,Jend
00455               do i=Istr,Iend
00456                 cff1=pn_v(i,j)
00457                 cff2=pm_v(i,j)
00458                 cff=0.25*( dUdz(i,j,k2)+dUdz(i+1,j,k2)
00459      &                +dUdz(i,j-1,k2)+dUdz(i+1,j-1,k2))
00460                 dnUdz=cff1*cff
00461                 dmUdz=cff2*cff
00462                 dnVdz=cff1*dVdz(i,j,k2)
00463                 dmVdz=cff2*dVdz(i,j,k2)
00464  
00465                 cff1=min(dZdx_p(i  ,j,k1),0.)
00466                 cff2=min(dZdx_p(i+1,j,k2),0.)
00467                 cff3=max(dZdx_p(i  ,j,k2),0.)
00468                 cff4=max(dZdx_p(i+1,j,k1),0.)
00469                 cff5=min(dZde_p(i  ,j,k1),0.)
00470                 cff6=min(dZde_p(i+1,j,k2),0.)
00471                 cff7=max(dZde_p(i  ,j,k2),0.)
00472                 cff8=max(dZde_p(i+1,j,k1),0.)
00473  
00474                 cff=om_v(i,j)*( cff1*(cff1*dnVdz-dnVdx(i  ,j,k1))
00475      &                         +cff2*(cff2*dnVdz-dnVdx(i+1,j,k2))
00476      &                         +cff3*(cff3*dnVdz-dnVdx(i  ,j,k2))
00477      &                         +cff4*(cff4*dnVdz-dnVdx(i+1,j,k1))
00478      &                                                          )
00479      &             +om_v(i,j)*( cff1*(cff5*dmUdz-dmUde(i  ,j,k1))
00480      &                         +cff2*(cff6*dmUdz-dmUde(i+1,j,k2))
00481      &                         +cff3*(cff7*dmUdz-dmUde(i  ,j,k2))
00482      &                         +cff4*(cff8*dmUdz-dmUde(i+1,j,k1))
00483      &                                                          )
00484                 cff1=min(dZde_r(i,j-1,k1),0.)
00485                 cff2=min(dZde_r(i,j  ,k2),0.)
00486                 cff3=max(dZde_r(i,j-1,k2),0.)
00487                 cff4=max(dZde_r(i,j  ,k1),0.)
00488                 cff5=min(dZdx_r(i,j-1,k1),0.)
00489                 cff6=min(dZdx_r(i,j  ,k2),0.)
00490                 cff7=max(dZdx_r(i,j-1,k2),0.)
00491                 cff8=max(dZdx_r(i,j  ,k1),0.)
00492  
00493                 cff=cff+om_v(i,j)*(
00494      &                          cff1*(cff1*dmVdz-dmVde(i,j-1,k1))
00495      &                         +cff2*(cff2*dmVdz-dmVde(i,j  ,k2))
00496      &                         +cff3*(cff3*dmVdz-dmVde(i,j-1,k2))
00497      &                         +cff4*(cff4*dmVdz-dmVde(i,j  ,k1))
00498      &                                                          )
00499      &             -om_v(i,j)*( cff1*(cff5*dnUdz-dnUdx(i,j-1,k1))
00500      &                         +cff2*(cff6*dnUdz-dnUdx(i,j  ,k2))
00501      &                         +cff3*(cff7*dnUdz-dnUdx(i,j-1,k2))
00502      &                         +cff4*(cff8*dnUdz-dnUdx(i,j  ,k1))
00503      &                                                          )
00504                 VFs(i,j,k2)=0.25*(visc2_r(i,j-1)+visc2_r(i,j))*cff
00505               enddo
00506             enddo
00507           endif
00508 !
00509 ! Apply viscous terms. Note that at this stage arrays u,v(...,3-nstp) 
00510 ! contain Hz*U and Hz*V with units of [m2/s].   Also compute vertical
00511 ! integral of viscous terms and add it into coupling terms for the 
00512 ! barotropic mode
00513 !
00514           do j=Jstr,Jend
00515             do i=IstrU,Iend
00516               cff=pn_u(i,j)*(UFx(i,j)-UFx(i-1,j))
00517      &           +pm_u(i,j)*(UFe(i,j+1)-UFe(i,j))
00518               cff1=pm_u(i,j)*pn_u(i,j)*cff
00519               rufrc(i,j)=rufrc(i,j) + cff
00520               u(i,j,k,indx)=u(i,j,k,indx) + dt*(cff1+UFs(i,j,k2)
00521      &                                              -UFs(i,j,k1))
00522 # ifdef DIAGNOSTICS_UV
00523               MHmix(i,j,k,1) = MHmix(i,j,k,1)+
00524      &                           2*(cff1+UFs(i,j,k2)-UFs(i,j,k1))
00525      &                                   /(Hz(i-1,j,k)+Hz(i,j,k))
00526 # endif
00527             enddo
00528           enddo
00529           do j=JstrV,Jend
00530             do i=Istr,Iend
00531               cff=pn_v(i,j)*(VFx(i+1,j)-VFx(i,j))
00532      &           -pm_v(i,j)*(VFe(i,j)-VFe(i,j-1))
00533               cff1=pm_v(i,j)*pn_v(i,j)*cff 
00534               rvfrc(i,j)=rvfrc(i,j) + cff
00535               v(i,j,k,indx)=v(i,j,k,indx) + dt*(cff1+VFs(i,j,k2)
00536      &                                              -VFs(i,j,k1))
00537 # ifdef DIAGNOSTICS_UV
00538               MHmix(i,j,k,2) = MHmix(i,j,k,2)+
00539      &                           2*(cff1+VFs(i,j,k2)-VFs(i,j,k1))
00540      &                                   /(Hz(i,j,k)+Hz(i,j-1,k))
00541 # endif
00542             enddo
00543           enddo
00544         endif
00545       enddo
00546       return
00547       end
00548 !====================================================================
00549 !
00550 # endif /* MIX_S_UV_SPG */
00551 
00552 #else
00553       subroutine uv3dmix_spg_empty
00554       end
00555 #endif
00556
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/uv3dmix_spg.F
