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 AGRIF
00006  
00007       subroutine uv3dmix_fine (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 uv3dmix_fine_tile (Istr,Iend,Jstr,Jend, 
00015      &                        A2d(1,1,trd), A2d(1,2,trd),
00016      &                        A2d(1,3,trd), A2d(1,4,trd))
00017       return
00018       end
00019  
00020       subroutine uv3dmix_fine_tile (Istr,Iend,Jstr,Jend, UFx,UFe,VFx,VFe)
00021 !
00022 ! Computes harmonic mixing of momentum, along constant S-surfaces
00023 ! as horizontal divergence of the stress tensor. Components of the
00024 ! stress tensor are:
00025 !                         du     dv
00026 !         s_xx = -s_yy = ---- - -----
00027 !                         dx     dy
00028 !
00029 !                         du     dv
00030 !         s_xy =  s_yx = ---- + ----
00031 !                         dy     dx
00032 !
00033 ! Reference:
00034 !
00035 ! Wajsowicz, R.C, 1993: A consistent formulation of the anisotropic
00036 !     stress tensor for use in models of the large-scale ocean
00037 !     circulation, JCP, 105, 333-338.
00038 !
00039 ! Sadourny, R. and K. Maynard, 1997: Formulations of lateral
00040 !     diffusion in geophysical fluid dynamics models, In "Numerical
00041 !     Methods of Atmospheric and Oceanic Modelling". Lin, Laprise,
00042 !     and Ritchie, Eds., NRC Research Press, 547-556.
00043 !
00044 ! Griffies, S.M. and R.W. Hallberg, 2000: Biharmonic friction with
00045 !     a Smagorinsky-like viscosity for use in large-scale eddy-
00046 !     permitting ocean models, Mon. Wea. Rev., 128, 8, 2935-2946.
00047 !
00048       implicit none
00049       integer Istr,Iend,Jstr,Jend, i,j,k, indx
00050       real UFe(PRIVATE_2D_SCRATCH_ARRAY),
00051      &     UFx(PRIVATE_2D_SCRATCH_ARRAY), cff,
00052      &     VFe(PRIVATE_2D_SCRATCH_ARRAY),
00053      &     VFx(PRIVATE_2D_SCRATCH_ARRAY)
00054      
00055 # include "param.h"
00056 # include "scalars.h"
00057 # include "grid.h"
00058 # include "ocean3d.h"
00059 # include "coupling.h"
00060 # include "mixing.h"
00061 # include "zoom.h"
00062 # ifdef DIAGNOSTICS_UV
00063 #  include "diagnostics.h"
00064 # endif
00065 !
00066 # ifdef MPI
00067 #  define LOCALLM Lmmpi
00068 #  define LOCALMM Mmmpi
00069 # else
00070 #  define LOCALLM Lm
00071 #  define LOCALMM Mm
00072 # endif
00073       real uspongeloc(PRIVATE_2D_SCRATCH_ARRAY,N)
00074       real vspongeloc(PRIVATE_2D_SCRATCH_ARRAY,N)
00075        integer decal
00076        real tinterp,onemtinterp, rrhot
00077        integer :: nold
00078        integer :: irhot
00079        
00080        external interpspongeu, interpspongev
00081 !
00082 # include "compute_auxiliary_bounds.h"
00083 
00084       irhot = Agrif_Irhot()
00085       rrhot = real(irhot)
00086       decal = 2*max(Agrif_Irhox(),Agrif_Irhoy())+1
00087                        
00088       If (nbcoarse == 1) THEN
00089       do k=1,N
00090       do j=JstrR,JendR
00091       do i=IstrR,IendR
00092       usponge(i,j,k) = 0.
00093       vsponge(i,j,k) = 0.
00094       enddo
00095       enddo
00096       enddo
00097       
00098 C$OMP BARRIER
00099 C$OMP MASTER
00100 
00101       Call Agrif_Set_bc(usponge,(/-decal,0/),
00102      &         InterpolationShouldbemade=.TRUE.)
00103       Call Agrif_Set_bc(vsponge,(/-decal,0/),
00104      &         InterpolationShouldbemade=.TRUE.)
00105 
00106        tinterp = 1.
00107                
00108       Call Agrif_Bc_Variable(usponge,usponge,calledweight=tinterp,
00109      &                         procname=interpspongeu)
00110       Call Agrif_Bc_Variable(vsponge,vsponge,calledweight=tinterp,
00111      &                         procname=interpspongev)
00112 
00113       UTimesponge = 3 - UTimesponge
00114 C$OMP END MASTER
00115 C$OMP BARRIER
00116 
00117        if (SOUTHERN_EDGE) then  
00118          do k=1,N
00119          do j=0,decal
00120          do i=Istr-1,Iend+1
00121           U_sponge_south(i,j,k,UTimesponge)=
00122      &              usponge(i,j,k)       
00123           V_sponge_south(i,j,k,UTimesponge)=
00124      &              vsponge(i,j,k)         
00125          enddo
00126          enddo
00127          enddo
00128        endif
00129                      
00130       if (NORTHERN_EDGE) then
00131          do k=1,N
00132          do j=Jend-decal,Jend+1
00133          do i=Istr-1,Iend+1
00134           U_sponge_north(i,j,k,UTimesponge)=
00135      &              usponge(i,j,k)          
00136           V_sponge_north(i,j,k,UTimesponge)=
00137      &              vsponge(i,j,k)     
00138          enddo
00139          enddo
00140          enddo
00141        endif               
00142 
00143 
00144       if (WESTERN_EDGE) then
00145          do k=1,N
00146          do j=Jstr-1,Jend+1
00147          do i=0,decal      
00148           U_sponge_west(i,j,k,UTimesponge)=
00149      &              usponge(i,j,k)       
00150           V_sponge_west(i,j,k,UTimesponge)=
00151      &              vsponge(i,j,k)        
00152          enddo
00153          enddo
00154          enddo
00155        endif
00156                     
00157       if (EASTERN_EDGE) then  
00158          do k=1,N
00159          do j=Jstr-1,Jend+1
00160          do i=Iend-decal,Iend+1    
00161           U_sponge_east(i,j,k,UTimesponge)=
00162      &              usponge(i,j,k)    
00163           V_sponge_east(i,j,k,UTimesponge)=
00164      &              vsponge(i,j,k)             
00165          enddo
00166          enddo
00167          enddo  
00168        endif
00169       
00170 C$OMP BARRIER  
00171       ENDIF    
00172 
00173       uspongeloc = 0.
00174       vspongeloc = 0.
00175 
00176        tinterp = real(nbcoarse-1)/rrhot
00177        IF (nbstep3d .LT. irhot) tinterp = 0.
00178        onemtinterp = -tinterp
00179        tinterp = 1.+tinterp
00180        
00181        nold = 3 - UTimesponge
00182 
00183        if (SOUTHERN_EDGE) then
00184        do k=1,N
00185        do j=0,decal
00186        do i=Istr-1,Iend+1
00187           uspongeloc(i,j,k) =
00188      &          (u(i,j,k,nstp)
00189      &       -onemtinterp*U_sponge_south(i,j,k,nold)
00190      &       -tinterp*U_sponge_south(i,j,k,UTimesponge))
00191      &           *umask(i,j)  
00192           vspongeloc(i,j,k) =
00193      &          (v(i,j,k,nstp)
00194      &       -onemtinterp*V_sponge_south(i,j,k,nold)
00195      &       -tinterp*V_sponge_south(i,j,k,UTimesponge))
00196      &           *vmask(i,j)            
00197        enddo
00198        enddo
00199        enddo
00200        endif
00201        
00202        if (NORTHERN_EDGE) then
00203        do k=1,N
00204        do j=Jend-decal,Jend+1
00205        do i=Istr-1,Iend+1
00206           uspongeloc(i,j,k) =
00207      &          (u(i,j,k,nstp)
00208      &       -onemtinterp*U_sponge_north(i,j,k,nold)
00209      &       -tinterp*U_sponge_north(i,j,k,UTimesponge))
00210      &           *umask(i,j)     
00211           vspongeloc(i,j,k) =
00212      &          (v(i,j,k,nstp)
00213      &       -onemtinterp*V_sponge_north(i,j,k,nold)
00214      &       -tinterp*V_sponge_north(i,j,k,UTimesponge))
00215      &           *vmask(i,j)        
00216        enddo
00217        enddo
00218        enddo           
00219        endif
00220        
00221       if (WESTERN_EDGE) then
00222          do k=1,N
00223          do j=Jstr-1,Jend+1
00224          do i=0,decal
00225           uspongeloc(i,j,k) =
00226      &          (u(i,j,k,nstp)
00227      &       -onemtinterp*U_sponge_west(i,j,k,nold)
00228      &       -tinterp*U_sponge_west(i,j,k,UTimesponge))
00229      &           *umask(i,j)
00230           vspongeloc(i,j,k) =
00231      &          (v(i,j,k,nstp)
00232      &       -onemtinterp*V_sponge_west(i,j,k,nold)
00233      &       -tinterp*V_sponge_west(i,j,k,UTimesponge))
00234      &           *vmask(i,j)     
00235          enddo
00236          enddo
00237          enddo
00238        endif      
00239        
00240       if (EASTERN_EDGE) then    
00241          do k=1,N
00242          do j=Jstr-1,Jend+1
00243          do i=Iend-decal,Iend+1    
00244           uspongeloc(i,j,k) =
00245      &          (u(i,j,k,nstp)
00246      &       -onemtinterp*U_sponge_east(i,j,k,nold)
00247      &       -tinterp*U_sponge_east(i,j,k,UTimesponge))
00248      &           *umask(i,j)    
00249           vspongeloc(i,j,k) =
00250      &          (v(i,j,k,nstp)
00251      &       -onemtinterp*V_sponge_east(i,j,k,nold)
00252      &       -tinterp*V_sponge_east(i,j,k,UTimesponge))
00253      &           *vmask(i,j)         
00254          enddo
00255          enddo
00256          enddo 
00257        endif
00258                           
00259       indx=3-nstp    !--> time index for target arrays;
00260 
00261             
00262 !
00263 !  Compute flux-components of the horizontal divergence of the stress
00264 !  tensor (m5/s2) in XI- and ETA-directions.
00265 !
00266       do k=1,N
00267         do j=JstrV-1,Jend
00268           do i=IstrU-1,Iend
00269             cff=0.5*Hz(i,j,k)*visc2_r(i,j)*(
00270      &    on_r(i,j)*pm(i,j)*( (pn(i  ,j)+pn(i+1,j))*uspongeloc(i+1,j,k)
00271      &                       -(pn(i-1,j)+pn(i  ,j))*uspongeloc(i  ,j,k)
00272      &                                                              )
00273      &   -om_r(i,j)*pn(i,j)*( (pm(i,j  )+pm(i,j+1))*vspongeloc(i,j+1,k)
00274      &                       -(pm(i,j-1)+pm(i,j  ))*vspongeloc(i,j  ,k)
00275      &                                                             ))
00276             UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
00277             VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
00278           enddo
00279         enddo
00280         do j=Jstr,Jend+1
00281           do i=Istr,Iend+1
00282             cff=0.125*visc2_p(i,j)*
00283      &      (Hz(i-1,j,k)+Hz(i,j,k)+Hz(i-1,j-1,k)+Hz(i,j-1,k))*(
00284      &      0.25*(pm(i-1,j)+pm(i,j)+pm(i-1,j-1)+pm(i,j-1))*on_p(i,j)
00285      &                  *( (pn(i  ,j-1)+pn(i  ,j))*vspongeloc(i  ,j,k)
00286      &                    -(pn(i-1,j-1)+pn(i-1,j))*vspongeloc(i-1,j,k)
00287      &                                                              )
00288      &     +0.25*(pn(i-1,j)+pn(i,j)+pn(i-1,j-1)+pn(i,j-1))*om_p(i,j)
00289      &                  *( (pm(i-1,j  )+pm(i,j  ))*uspongeloc(i,j  ,k)
00290      &                  -(pm(i-1,j-1)+pm(i,j-1))*uspongeloc(i,j-1,k)
00291      &                                                             ))
00292 # ifdef MASKING
00293      &                                                   *pmask(i,j)
00294 # endif
00295             UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
00296             VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
00297           enddo
00298         enddo
00299 !
00300 ! Apply viscous terms. Note that at this stage arrays u,v(...,3-nstp)
00301 ! contain Hz*U and Hz*V with units of [m2/s].   Also compute vertical
00302 ! integral of viscous terms and add it into coupling terms for the
00303 ! barotropic mode
00304 !
00305         do j=Jstr,Jend
00306           do i=IstrU,Iend
00307             cff=0.125*(pm(i-1,j)+pm(i,j))*(pn(i-1,j) +pn(i,j))
00308      &             *( (pn(i-1,j)+pn(i,j))*(UFx(i,j)-UFx(i-1,j))
00309      &               +(pm(i-1,j)+pm(i,j))*(UFe(i,j+1)-UFe(i,j))
00310      &                                                        )
00311             rufrc(i,j)=rufrc(i,j) + cff
00312             u(i,j,k,indx)=u(i,j,k,indx) + dt*cff
00313 # ifdef DIAGNOSTICS_UV
00314             MHmix(i,j,k,1) = cff*om_u(i,j)*on_u(i,j)
00315 # endif
00316           enddo
00317         enddo
00318         do j=JstrV,Jend
00319           do i=Istr,Iend
00320             cff=0.125*(pm(i,j)+pm(i,j-1))*(pn(i,j) +pn(i,j-1))
00321      &             *( (pn(i,j-1)+pn(i,j))*(VFx(i+1,j)-VFx(i,j))
00322      &               -(pm(i,j-1)+pm(i,j))*(VFe(i,j)-VFe(i,j-1))
00323      &                                                        )
00324             rvfrc(i,j)=rvfrc(i,j) + cff
00325             v(i,j,k,indx)=v(i,j,k,indx) + dt*cff
00326 # ifdef DIAGNOSTICS_UV
00327             MHmix(i,j,k,2) = cff*om_v(i,j)*on_v(i,j)
00328 # endif
00329           enddo
00330         enddo
00331       enddo
00332       return
00333       end
00334       
00335       subroutine interpspongeu(tabres,i1,i2,j1,j2,k1,k2)
00336       implicit none
00337 # include "param.h"
00338 # include "ocean3d.h"
00339 # include "scalars.h"
00340 
00341        integer i1,i2,j1,j2,k1,k2
00342        real tabres(i1:i2,j1:j2,k1:k2)
00343        
00344        tabres(i1:i2,j1:j2,k1:k2) = u(i1:i2,j1:j2,k1:k2,nstp)
00345              
00346       return
00347       end
00348       
00349       subroutine interpspongev(tabres,i1,i2,j1,j2,k1,k2)
00350       implicit none
00351 # include "param.h"
00352 # include "ocean3d.h"
00353 # include "scalars.h"
00354 
00355        integer i1,i2,j1,j2,k1,k2
00356        real tabres(i1:i2,j1:j2,k1:k2)
00357        
00358        tabres(i1:i2,j1:j2,k1:k2) = v(i1:i2,j1:j2,k1:k2,nstp)
00359              
00360       return
00361       end      
00362             
00363 #else
00364       subroutine nestingvisc3d_empty
00365       end
00366 #endif
00367
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/uv3dmix_fine.F
