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

00001 !
00002 ! $Id: pre_step3d.F,v 1.7 2005/10/11 08:37:32 pmarches Exp $
00003 !
00004 #include "cppdefs.h"
00005 #ifdef SOLVE3D
00006 
00007       subroutine pre_step3d (tile)
00008 !
00009 ! version pre_step3d4S
00010 !
00011       implicit none
00012       integer tile,  trd,omp_get_thread_num
00013 # include "param.h"
00014 # include "private_scratch.h"
00015 # include "ocean3d.h"
00016 # include "compute_tile_bounds.h"
00017       trd=omp_get_thread_num()
00018       call pre_step3d_tile (Istr,Iend,Jstr,Jend,
00019      &                              A3d(1,1,trd), A3d(1,2,trd),
00020      &                A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd),
00021      &                A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd),
00022      &                                            A3d(1,3,trd),
00023      &                u(START_2D_ARRAY,1,1),
00024      &                v(START_2D_ARRAY,1,1),
00025      &                t(START_2D_ARRAY,1,1,1),
00026      &                Hz(START_2D_ARRAY,1),
00027      &                Hz_bak(START_2D_ARRAY,1),
00028      &                Huon(START_2D_ARRAY,1),
00029      &                Hvom(START_2D_ARRAY,1),
00030      &                W(START_2D_ARRAY,0))
00031       return
00032       end
00033 
00034       subroutine pre_step3d_tile (Istr,Iend,Jstr,Jend, ru,rv,
00035      &                         FC,CF,DC, FX,FE,WORK, Hz_half,
00036      &                           u,v,t,Hz,Hz_bak,Huon,Hvom,W)
00037 !
00038 !--------------------------------------------------------------------
00039 ! Preliminary step: initialize computations of the new time step
00040 ! 3D primitive variables.
00041 !
00042 ! Since r.h.s. arrays ru,rv,rt(:,:,???[,:]), which at this moment
00043 ! contain r.h.s at time step n-2 will be overwritten by the
00044 ! subsequent routines within rhs3d driver, both [n-1 and n-2] old-
00045 ! time-step r.h.s. term in Adams-Bashforth stepping scheme are
00046 ! added at this time to the time step [n] fields and the result is
00047 ! stored as u,v,t(:,:,???[,:]).
00048 !
00049 ! The actual time step will be completed in step3d, after the
00050 ! time step [n] r.h.s. terms and new-time step Hz will be available
00051 ! after the completion rhs3d computations and the 2D (barotropic
00052 ! mode) computations.
00053 !--------------------------------------------------------------------
00054 !
00055       implicit none
00056 # include "param.h"
00057       integer Istr,Iend,Jstr,Jend, itrc, i,j,k, indx
00058 # ifdef PSOURCE
00059      &       ,is
00060 # endif
00061 # ifdef MPI
00062      &       ,imin,imax,jmin,jmax
00063 #endif
00064       real ru(PRIVATE_2D_SCRATCH_ARRAY,N),    cff,
00065      &     rv(PRIVATE_2D_SCRATCH_ARRAY,N),    cff1,
00066      &     FC(PRIVATE_1D_SCRATCH_ARRAY,0:N),  cff2,
00067      &     CF(PRIVATE_1D_SCRATCH_ARRAY,0:N),
00068      &     DC(PRIVATE_1D_SCRATCH_ARRAY,0:N),  gamma, 
00069      &         FX(PRIVATE_2D_SCRATCH_ARRAY),  epsil,
00070      &         FE(PRIVATE_2D_SCRATCH_ARRAY),
00071      &       WORK(PRIVATE_2D_SCRATCH_ARRAY),
00072      &      Hz_half(PRIVATE_2D_SCRATCH_ARRAY,N),
00073      &                     u(GLOBAL_2D_ARRAY,N,3),
00074      &                     v(GLOBAL_2D_ARRAY,N,3),
00075      &                     t(GLOBAL_2D_ARRAY,N,3,NT),
00076      &                     Hz(GLOBAL_2D_ARRAY,N),
00077      &                     Hz_bak(GLOBAL_2D_ARRAY,N),
00078      &                     Huon(GLOBAL_2D_ARRAY,N),
00079      &                     Hvom(GLOBAL_2D_ARRAY,N),
00080      &                     W(GLOBAL_2D_ARRAY,0:N)
00081       parameter (gamma=1./6., epsil=1.E-16)
00082 # include "grid.h"
00083 C--# include "ocean3d.h" --> u declared only once
00084 # include "coupling.h"
00085 # include "ocean2d.h"
00086 # include "forces.h"
00087 # include "mixing.h"
00088 # include "scalars.h"
00089 # ifdef PSOURCE
00090 #  include "sources.h"
00091 # endif
00092 !
00093 # include "compute_auxiliary_bounds.h"
00094 !
00095       indx=3-nstp
00096 
00097       if (FIRST_TIME_STEP) then
00098         cff=0.5*dt
00099         cff1=1.
00100         cff2=0.
00101       else
00102         cff=(1.-gamma)*dt
00103         cff1=0.5+gamma
00104         cff2=0.5-gamma
00105       endif
00106       do k=1,N
00107         do j=JstrV-1,Jend
00108           do i=IstrU-1,Iend
00109             Hz_half(i,j,k)=cff1*Hz(i,j,k)+cff2*Hz_bak(i,j,k)
00110      &        -cff*pm(i,j)*pn(i,j)*( Huon(i+1,j,k)-Huon(i,j,k)
00111      &                              +Hvom(i,j+1,k)-Hvom(i,j,k)
00112      &                                    +W(i,j,k)-W(i,j,k-1)
00113      &                                                       )
00114           enddo
00115         enddo
00116       enddo
00117 
00118 !
00119 !  Start computation of the auxiliary tracer field.
00120 !----------------------------------------------------
00121 !  Once it will be completed, t(:,:,???,:) is effectively halfway
00122 !  between time steps n and n+1. A high spatial order, centered,
00123 !  non-conservative [but constancy preserving!] scheme is used for
00124 !  this auxiliary step. This is done by introducing an artificial
00125 !  continuity equation [Easter, 1993].
00126 !
00127 !  Since this field will be used exclussively to compute the high-
00128 !  order fluxes during subsequent step3d_t operation, the final
00129 !  values of t(i,j,k,??,itrc) alfer step3d_t will be computed in
00130 !  a flux-conservative manner. The overall time step will be both 
00131 !  conservative and constancy preserving.
00132 !
00133 !  This preliminary step shall be done before field t(:,:,:,???,:)
00134 !  loses its meaningful values during the pre-step operation. 
00135 !
00136 !
00137 !---------------------------------------------------------------
00138 ! Compute horizontal advection
00139 !---------------------------------------------------------------
00140 !
00141 # ifdef HADV_UPSTREAM_TS
00142 #  define curv WORK
00143 # else
00144 #  define grad WORK
00145 # endif
00146       do k=1,N
00147         do itrc=1,NT
00148 # ifdef EW_PERIODIC
00149 #  define I_EXT_RANGE Istr-1,Iend+2
00150 # else
00151 #  ifdef MPI
00152           if (WEST_INTER) then
00153             imin=Istr-1
00154           else
00155             imin=max(Istr-1,1)
00156           endif
00157           if (EAST_INTER) then
00158             imax=Iend+2
00159           else
00160             imax=min(Iend+2,Lmmpi+1)
00161         endif
00162 #   define I_EXT_RANGE imin,imax
00163 #  else
00164 #   define I_EXT_RANGE max(Istr-1,1),min(Iend+2,Lm+1)
00165 #  endif
00166 # endif
00167           do j=Jstr,Jend
00168             do i=I_EXT_RANGE
00169               FX(i,j)=(t(i,j,k,nstp,itrc)-t(i-1,j,k,nstp,itrc))
00170 # ifdef MASKING
00171      &                                               *umask(i,j)
00172 # endif
00173             enddo
00174           enddo
00175 #  undef I_EXT_RANGE
00176 # ifndef EW_PERIODIC 
00177           if (WESTERN_EDGE) then
00178             do j=Jstr,Jend
00179               FX(0,j)=FX(1,j)
00180             enddo
00181           endif
00182           if (EASTERN_EDGE) then
00183 #  ifdef MPI
00184             do j=Jstr,Jend
00185               FX(Lmmpi+2,j)=FX(Lmmpi+1,j)
00186             enddo
00187 #  else
00188              do j=Jstr,Jend
00189               FX(Lm+2,j)=FX(Lm+1,j)
00190             enddo
00191 #  endif
00192           endif
00193 # endif
00194           do j=Jstr,Jend
00195             do i=Istr-1,Iend+1
00196 # if defined HADV_UPSTREAM_TS
00197               curv(i,j)=FX(i+1,j)-FX(i,j)
00198 # elif defined HADV_AKIMA_TS
00199               cff=2.*FX(i+1,j)*FX(i,j)
00200               if (cff.gt.epsil) then
00201                 grad(i,j)=cff/(FX(i+1,j)+FX(i,j))
00202               else
00203                 grad(i,j)=0.
00204               endif
00205 # else
00206               grad(i,j)=0.5*(FX(i+1,j)+FX(i,j))
00207 # endif
00208             enddo
00209           enddo             !--> discard FX
00210 
00211           do j=Jstr,Jend
00212             do i=Istr,Iend+1
00213 # ifdef HADV_UPSTREAM_TS
00214               if (Huon(i,j,k) .gt. 0.) then
00215                 cff=curv(i-1,j)
00216               else
00217                 cff=curv(i,j)
00218               endif
00219               FX(i,j)=0.5*( t(i,j,k,nstp,itrc)+t(i-1,j,k,nstp,itrc)
00220      &                           -0.333333333333*cff )*Huon(i,j,k)
00221 # else
00222               FX(i,j)=0.5*( t(i,j,k,nstp,itrc)+t(i-1,j,k,nstp,itrc)
00223      &                     -0.333333333333*(grad(i,j)-grad(i-1,j))
00224      &                                               )*Huon(i,j,k)
00225 # endif
00226             enddo
00227           enddo            !--> discard grad
00228 
00229 # ifdef NS_PERIODIC
00230 #  define J_EXT_RANGE Jstr-1,Jend+2
00231 # else
00232 #  ifdef MPI
00233           if (SOUTH_INTER) then
00234             jmin=Jstr-1
00235           else
00236             jmin=max(Jstr-1,1)
00237           endif
00238           if (NORTH_INTER) then
00239             jmax=Jend+2
00240           else
00241             jmax=min(Jend+2,Mmmpi+1)
00242           endif
00243 #   define J_EXT_RANGE jmin,jmax
00244 #  else
00245 #   define J_EXT_RANGE max(Jstr-1,1),min(Jend+2,Mm+1)
00246 #  endif
00247 # endif
00248 
00249           do j=J_EXT_RANGE
00250             do i=Istr,Iend
00251               FE(i,j)=(t(i,j,k,nstp,itrc)-t(i,j-1,k,nstp,itrc))
00252 # ifdef MASKING
00253      &                                               *vmask(i,j)
00254 # endif
00255             enddo
00256           enddo
00257 # undef J_EXT_RANGE
00258 # ifndef NS_PERIODIC
00259           if (SOUTHERN_EDGE) then
00260             do i=Istr,Iend
00261               FE(i,0)=FE(i,1)
00262             enddo
00263           endif
00264           if (NORTHERN_EDGE) then
00265 #  ifdef MPI
00266             do i=Istr,Iend
00267               FE(i,Mmmpi+2)=FE(i,Mmmpi+1)
00268             enddo
00269 #  else
00270             do i=Istr,Iend
00271               FE(i,Mm+2)=FE(i,Mm+1)
00272             enddo
00273 #  endif
00274           endif
00275 # endif
00276 
00277           do j=Jstr-1,Jend+1
00278             do i=Istr,Iend
00279 # if defined HADV_UPSTREAM_TS
00280               curv(i,j)=FE(i,j+1)-FE(i,j)
00281 # elif defined HADV_AKIMA_TS
00282               cff=2.*FE(i,j+1)*FE(i,j)
00283               if (cff.gt.epsil) then
00284                 grad(i,j)=cff/(FE(i,j+1)+FE(i,j))
00285               else
00286                 grad(i,j)=0.
00287               endif
00288 # else
00289               grad(i,j)=0.5*(FE(i,j+1)+FE(i,j))
00290 # endif
00291             enddo
00292           enddo            !--> discard FE
00293 
00294           do j=Jstr,Jend+1
00295             do i=Istr,Iend
00296 # ifdef HADV_UPSTREAM_TS
00297               if (Hvom(i,j,k) .gt. 0.) then
00298                 cff=curv(i,j-1)
00299               else
00300                 cff=curv(i,j)
00301               endif
00302               FE(i,j)=0.5*( t(i,j,k,nstp,itrc)+t(i,j-1,k,nstp,itrc)
00303      &                          -0.333333333333*cff )*Hvom(i,j,k)
00304 #  undef curv
00305 # else
00306               FE(i,j)=0.5*( t(i,j,k,nstp,itrc)+t(i,j-1,k,nstp,itrc)
00307      &                     -0.333333333333*(grad(i,j)-grad(i,j-1))
00308      &                                               )*Hvom(i,j,k)
00309 #  undef grad
00310 # endif
00311             enddo
00312           enddo            !--> discard grad
00313 
00314 # ifdef PSOURCE
00315 !
00316 !-------------------------------------------------------------------
00317 ! Apply point sources for river runoff simulations
00318 !-------------------------------------------------------------------
00319 !
00320           do is=1,Nsrc
00321             i=Isrc(is)
00322             j=Jsrc(is)
00323             if (Istr.le.i .and. i.le.Iend+1
00324      &                   .and. Jstr.le.j .and. j.le.Jend+1) then
00325               if (Dsrc(is).eq.0) then
00326                 if (Lsrc(is,itrc)) then
00327                   FX(i,j)=Huon(i,j,k)*Tsrc(is,k,itrc)
00328 #   ifdef MASKING
00329                 else
00330                   if (rmask(i,j).eq.0 .and. rmask(i-1,j).eq.1) then
00331                     FX(i,j)=Huon(i,j,k)*t(i-1,j,k,nstp,itrc)
00332                   elseif (rmask(i,j).eq.1. .and. rmask(i-1,j).eq.0) then
00333                     FX(i,j)=Huon(i,j,k)*t(i  ,j,k,nstp,itrc)
00334                   endif
00335 #   endif
00336                 endif
00337               else
00338                 if (Lsrc(is,itrc)) then
00339                   FE(i,j)=Hvom(i,j,k)*Tsrc(is,k,itrc)
00340 #   ifdef MASKING
00341                 else
00342                   if (rmask(i,j).eq.0 .and. rmask(i,j-1).eq.1) then
00343                     FE(i,j)=Hvom(i,j,k)*t(i,j-1,k,nstp,itrc)
00344                   elseif (rmask(i,j).eq.1 .and. rmask(i,j-1).eq.0) then
00345                     FE(i,j)=Hvom(i,j,k)*t(i,j  ,k,nstp,itrc)
00346                   endif
00347 #   endif
00348                 endif
00349               endif
00350             endif
00351           enddo
00352 # endif
00353 
00354           if (FIRST_TIME_STEP) then
00355             cff=0.5*dt
00356             do j=Jstr,Jend
00357               do i=Istr,Iend
00358                 t(i,j,k,nnew,itrc)=Hz(i,j,k)*t(i,j,k,nstp,itrc)
00359      &                      -cff*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
00360      &                                            +FE(i,j+1)-FE(i,j))
00361               enddo
00362             enddo
00363           else
00364             cff=(1.-gamma)*dt
00365             cff1=0.5+gamma
00366             cff2=0.5-gamma
00367             do j=Jstr,Jend
00368               do i=Istr,Iend
00369                 t(i,j,k,nnew,itrc)=cff1*Hz(i,j,k)*t(i,j,k,nstp,itrc)
00370      &                         +cff2*Hz_bak(i,j,k)*t(i,j,k,indx,itrc)
00371      &                      -cff*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
00372      &                                            +FE(i,j+1)-FE(i,j))
00373               enddo
00374             enddo
00375           endif
00376         enddo   ! <-- itrc
00377       enddo       ! <-- k
00378 
00379 !
00380 ! Continue computation of the auxiliary tracer field: auxiliary
00381 ! continuity equation (the same for all tracers): DC=1/Hz_half_new,
00382 ! where Hz_half_new is Hz at time step n+1/2 as it would be computed
00383 ! from three-dimensional divergence of volume fluxes Huon,Hvom and W.
00384 !
00385 !---------------------------------------------------------------
00386 ! Compute vertical advection
00387 !---------------------------------------------------------------
00388 ! Finalize computation of the auxiliary tracer field: compute its
00389 ! change due to vertical advection. Computation of vertical advective
00390 ! fluxes requires interpolation of tracer values to the verical grid-
00391 ! box interfaces (W-points). This can be done by either using
00392 ! parabolic spline interpolation or, more simple local cubic
00393 ! polynomial [linear interpolation is considered obsolete].
00394 !
00395       do j=Jstr,Jend
00396         do k=1,N
00397           do i=Istr,Iend
00398             DC(i,k)=1./Hz_half(i,j,k)
00399           enddo
00400         enddo
00401         do itrc=1,NT
00402 # ifdef VADV_SPLINES_TS
00403                               ! Construct parabolic splines: here
00404           do i=Istr,Iend      ! CF is the set of vertical derivatives
00405             FC(i,0)=0.        ! of the tracer field t(:,:,:,nstp,:),
00406             CF(i,0)=0.        ! FC is an auxiliary scratch variable. 
00407           enddo
00408           do k=1,N-1,+1
00409             do i=Istr,Iend
00410               cff=1./(2.*Hz(i,j,k+1)+Hz(i,j,k)*(2.-FC(i,k-1)))
00411               FC(i,k)=cff*Hz(i,j,k+1)
00412               CF(i,k)=cff*( 6.*( t(i,j,k+1,nstp,itrc)
00413      &               -t(i,j,k,nstp,itrc))-Hz(i,j,k)*CF(i,k-1))
00414             enddo
00415           enddo
00416           do i=Istr,Iend
00417             CF(i,N)=0.
00418           enddo
00419           do k=N-1,1,-1       !<-- irreversible
00420             do i=Istr,Iend
00421               CF(i,k)=CF(i,k)-FC(i,k)*CF(i,k+1)
00422             enddo
00423           enddo               !--> discard FC, keep CF
00424 
00425           cff=1./3.           ! Compute vertical advective fluxes 
00426           do k=1,N-1          ! FC=W*[spline-interpolated tracer]
00427             do i=Istr,Iend
00428               FC(i,k)=W(i,j,k)*( t(i,j,k,nstp,itrc)+cff*Hz(i,j,k)
00429      &                                  *(CF(i,k)+0.5*CF(i,k-1)) )
00430             enddo
00431           enddo               !--> discard CF
00432           do i=Istr,Iend
00433             FC(i,N)=0.
00434             FC(i,0)=0.
00435           enddo
00436 # elif defined VADV_AKIMA_TS
00437 !
00438 ! Compute vertical advective fluxes using 4th-order Akima scheme
00439 !
00440           do k=1,N-1
00441             do i=istr,iend
00442               FC(i,k)=t(i,j,k+1,nstp,itrc)-t(i,j,k,nstp,itrc)
00443             enddo
00444           enddo
00445           do i=istr,iend
00446             FC(i,0)=FC(i,1)
00447             FC(i,N)=FC(i,N-1)
00448           enddo
00449           do k=1,N
00450             do i=istr,iend
00451               cff=2.*FC(i,k)*FC(i,k-1)
00452               if (cff.gt.epsil) then
00453                 CF(i,k)=cff/(FC(i,k)+FC(i,k-1))
00454               else
00455                 CF(i,k)=0.
00456               endif
00457             enddo
00458           enddo            !--> discard FC
00459           do k=1,N-1
00460             do i=istr,iend
00461               FC(i,k)=0.5*( t(i,j,k,nstp,itrc)+t(i,j,k+1,nstp,itrc)
00462      &                -0.333333333333*(CF(i,k+1)-CF(i,k)) )*W(i,j,k)
00463             enddo
00464           enddo            !--> discard CF
00465           do i=istr,iend
00466             FC(i,0)=0.
00467             FC(i,N)=0.
00468           enddo
00469 # elif VADV_C2_TS
00470 !
00471 ! Compute vertical advective fluxes using 2th-order centered scheme
00472 !
00473           do k=1,N-1
00474            do i=Istr,Iend
00475              FC(i,k)=0.5*W(i,j,k)*(t( i,j,k,nstp,itrc)
00476      &                           +t(i,j,k+1,nstp,itrc))
00477            enddo
00478          enddo
00479          do i=Istr,Iend
00480            FC(i, 0)=0.
00481            FC(i,N )=0.
00482          enddo
00483 # else
00484 !
00485 ! Compute vertical advective fluxes using 4th-order centered scheme
00486 !
00487           do k=2,N-2
00488             do i=Istr,Iend
00489               FC(i,k)=W(i,j,k)*(
00490      &                     0.58333333333333*( t(i,j,k  ,nstp,itrc)
00491      &                                       +t(i,j,k+1,nstp,itrc))
00492      &                    -0.08333333333333*( t(i,j,k-1,nstp,itrc)
00493      &                                       +t(i,j,k+2,nstp,itrc))
00494      &                                                            )
00495             enddo
00496           enddo
00497           do i=Istr,Iend
00498             FC(i, 0)=0.0
00499             FC(i,  1)=W(i,j,  1)*(     0.5*t(i,j,  1,nstp,itrc)
00500      &                       +0.58333333333333*t(i,j,  2,nstp,itrc)
00501      &                       -0.08333333333333*t(i,j,  3,nstp,itrc)
00502      &                                                            )
00503             FC(i,N-1)=W(i,j,N-1)*(     0.5*t(i,j,N  ,nstp,itrc)
00504      &                       +0.58333333333333*t(i,j,N-1,nstp,itrc)
00505      &                       -0.08333333333333*t(i,j,N-2,nstp,itrc)
00506      &                                                            )
00507             FC(i,N )=0.0
00508           enddo
00509 # endif
00510 
00511           if (FIRST_TIME_STEP) then
00512             cff=0.5*dt
00513           else
00514             cff=(1.-gamma)*dt
00515           endif
00516           do k=1,N
00517             do i=Istr,Iend
00518               t(i,j,k,nnew,itrc)=DC(i,k)*( t(i,j,k,nnew,itrc)
00519      &               -cff*pm(i,j)*pn(i,j)*(FC(i,k)-FC(i,k-1)))
00520 c**
00521 c**   if (itrc.eq.2) t(i,j,k,nnew,2)=1.! CONSTANCY PRESERVATION CHECK
00522 c**
00523 #ifdef CONST_TRACERS 
00524               t(i,j,k,nnew,itrc)=t(i,j,k,nstp,itrc)
00525 #endif
00526             enddo
00527           enddo            !--> discard FC
00528         enddo   !<-- itrc  !--> discard DC
00529 
00530 !
00531 !---------------------------------------------------------------
00532 ! Momentum equations: time stepping to time n+1/2
00533 !---------------------------------------------------------------
00534 !
00535         do i=IstrU,Iend
00536           DC(i,0)=0.25*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
00537         enddo
00538         if (FIRST_TIME_STEP) then
00539           do k=1,N
00540             do i=IstrU,Iend
00541               u(i,j,k,nnew)=( u(i,j,k,nstp)*(Hz(i,j,k)+Hz(i-1,j,k))
00542      &                                         +dt*DC(i,0)*ru(i,j,k)
00543      &                           )/(Hz_half(i,j,k)+Hz_half(i-1,j,k))
00544 
00545               u(i,j,k,indx)=u(i,j,k,nstp)*0.5*(Hz(i,j,k)+
00546      &                                         Hz(i-1,j,k))
00547             enddo
00548           enddo
00549         else
00550           cff=2.*(1.-gamma)*dt
00551           cff1=0.5+gamma
00552           cff2=0.5-gamma
00553           do k=1,N
00554             do i=IstrU,Iend
00555               u(i,j,k,nnew)=( cff1*u(i,j,k,nstp)*(Hz(i,j,k)+
00556      &                                            Hz(i-1,j,k))
00557      &                       +cff2*u(i,j,k,indx)*(Hz_bak(i,j,k)+
00558      &                                            Hz_bak(i-1,j,k))
00559      &                       +cff*DC(i,0)*ru(i,j,k)
00560      &                         )/(Hz_half(i,j,k)+Hz_half(i-1,j,k))
00561 
00562               u(i,j,k,indx)=u(i,j,k,nstp)*0.5*(Hz(i,j,k)+
00563      &                                         Hz(i-1,j,k))
00564             enddo
00565           enddo
00566         endif
00567 
00568         if (j.ge.JstrV) then
00569           do i=Istr,Iend
00570             DC(i,0)=0.25*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
00571           enddo
00572           if (FIRST_TIME_STEP) then
00573             do k=1,N
00574               do i=Istr,Iend
00575                 v(i,j,k,nnew)=(v(i,j,k,nstp)*(Hz(i,j,k)+Hz(i,j-1,k))
00576      &                                          +dt*DC(i,0)*rv(i,j,k)
00577      &                            )/(Hz_half(i,j,k)+Hz_half(i,j-1,k))
00578 
00579                 v(i,j,k,indx)=v(i,j,k,nstp)*0.5*(Hz(i,j,k)+
00580      &                                           Hz(i,j-1,k))
00581               enddo
00582             enddo
00583           else
00584             cff=2.*(1.-gamma)*dt
00585             cff1=0.5+gamma
00586             cff2=0.5-gamma
00587             do k=1,N
00588               do i=Istr,Iend
00589                 v(i,j,k,nnew)=( cff1*v(i,j,k,nstp)*(Hz(i,j,k)+
00590      &                                              Hz(i,j-1,k))
00591      &                         +cff2*v(i,j,k,indx)*(Hz_bak(i,j,k)+
00592      &                                              Hz_bak(i,j-1,k))
00593      &                         +cff*DC(i,0)*rv(i,j,k)
00594      &                           )/(Hz_half(i,j,k)+Hz_half(i,j-1,k))
00595 
00596                 v(i,j,k,indx)=v(i,j,k,nstp)*0.5*(Hz(i,j,k)+
00597      &                                           Hz(i,j-1,k))
00598               enddo
00599             enddo               !--> discard DC(:,0) 
00600           endif
00601         endif
00602       enddo     !<-- j
00603 
00604 !
00605 !---------------------------------------------------------------
00606 ! Set PHYSICAL lateral boundary conditions.
00607 !---------------------------------------------------------------
00608 !
00609       do itrc=1,NT
00610         call t3dbc_tile (Istr,Iend,Jstr,Jend, nnew,itrc, WORK)
00611 # if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
00612         call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
00613      &                          t(START_2D_ARRAY,1,nnew,itrc))
00614 # endif
00615       enddo
00616 
00617       call u3dbc_tile (Istr,Iend,Jstr,Jend, WORK)
00618       call v3dbc_tile (Istr,Iend,Jstr,Jend, WORK)
00619 
00620 !
00621 !---------------------------------------------------------------
00622 ! Coupling, include ghost points associated with PHYSICAL
00623 ! boundaries ONLY. Do not touch periodic ghost points or
00624 ! internal computational margins (MPI code). 
00625 !---------------------------------------------------------------
00626 !
00627 # ifdef EW_PERIODIC
00628 #  define IU_RANGE Istr,Iend
00629 #  define IV_RANGE Istr,Iend
00630 # else
00631 #  define IU_RANGE Istr,IendR
00632 #  define IV_RANGE IstrR,IendR
00633 # endif
00634 
00635 # ifdef NS_PERIODIC
00636 #  define JU_RANGE Jstr,Jend
00637 #  define JV_RANGE Jstr,Jend
00638 # else
00639 #  define JU_RANGE JstrR,JendR
00640 #  define JV_RANGE Jstr,JendR
00641 # endif
00642 !
00643       do j=JU_RANGE                      ! XI-component
00644         do i=IU_RANGE
00645           DC(i,0)=0.
00646           CF(i,0)=0.
00647         enddo
00648         do k=1,N,+1
00649           do i=IU_RANGE
00650             DC(i,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*on_u(i,j)
00651             DC(i,0)=DC(i,0)+DC(i,k)
00652             CF(i,0)=CF(i,0)+DC(i,k)*u(i,j,k,nnew)
00653           enddo
00654         enddo
00655         if (FIRST_TIME_STEP) then
00656           cff1=1.
00657           cff2=0.
00658         else
00659           cff1=1.5
00660           cff2=0.5
00661         endif
00662         do i=IU_RANGE
00663           CF(i,0)=( CF(i,0)-cff1*DU_avg1(i,j,nstp)
00664      &             +cff2*DU_avg1(i,j,indx) )/DC(i,0)
00665         enddo
00666         do k=N,1,-1
00667           do i=IU_RANGE
00668             u(i,j,k,nnew)=(u(i,j,k,nnew)-CF(i,0))
00669 # ifdef MASKING
00670      &                                 *umask(i,j)
00671 # endif
00672           enddo
00673         enddo
00674       enddo
00675 
00676       do j=JV_RANGE
00677           do i=IV_RANGE                   ! ETA-component
00678             DC(i,0)=0.
00679             CF(i,0)=0.
00680           enddo
00681           do k=1,N,+1
00682             do i=IV_RANGE
00683               DC(i,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*om_v(i,j)
00684               DC(i,0)=DC(i,0)+DC(i,k)
00685               CF(i,0)=CF(i,0)+DC(i,k)*v(i,j,k,nnew)
00686             enddo
00687           enddo
00688           if (FIRST_TIME_STEP) then
00689             cff1=1.
00690             cff2=0.
00691           else
00692             cff1=1.5
00693             cff2=0.5
00694           endif
00695           do i=IV_RANGE
00696             CF(i,0)=( CF(i,0)-cff1*DV_avg1(i,j,nstp)
00697      &               +cff2*DV_avg1(i,j,indx) )/DC(i,0)
00698           enddo
00699           do k=N,1,-1
00700             do i=IV_RANGE
00701               v(i,j,k,nnew)=(v(i,j,k,nnew)-CF(i,0))
00702 # ifdef MASKING
00703      &                                  *vmask(i,j)
00704 # endif
00705             enddo
00706           enddo
00707       enddo
00708 
00709 # if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
00710       call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
00711      &                                   u(START_2D_ARRAY,1,nnew))
00712       call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
00713      &                                   v(START_2D_ARRAY,1,nnew))
00714 # endif
00715 
00716 !
00717 !---------------------------------------------------------------
00718 ! Prepare to start two-dimensional time stepping:
00719 ! set the initial values of the fast-time-step free-surface
00720 ! field to its fast-time-averaged values corresponding
00721 ! o the time step n (nstp).
00722 !---------------------------------------------------------------
00723 !
00724 # ifdef MPI
00725       do j=JU_RANGE
00726         do i=IV_RANGE
00727           zeta(i,j,knew)=Zt_avg1(i,j)
00728         enddo
00729       enddo
00730       call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
00731      &                          zeta(START_2D_ARRAY,knew))
00732 # else
00733       do j=JstrR,JendR
00734         do i=IstrR,IendR
00735           zeta(i,j,knew)=Zt_avg1(i,j)
00736         enddo
00737       enddo
00738 # endif
00739 !
00740 # undef IU_RANGE
00741 # undef JU_RANGE
00742 # undef IV_RANGE
00743 # undef JV_RANGE
00744 #else
00745       subroutine pre_step3d_empty
00746 #endif
00747       return
00748       end
/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/pre_step3d.F
