bdy_3d.F90 13.9 KB
Newer Older
bjb's avatar
bjb committed
1
!$Id: bdy_3d.F90,v 1.13 2009-09-04 07:58:19 bjb Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
6
! !MODULE:  bdy_3d - 3D boundary conditions \label{bdy-3d}
gotm's avatar
gotm committed
7 8 9 10 11
!
! !INTERFACE:
   module bdy_3d
!
! !DESCRIPTION:
12 13 14
!  
! Here, the three-dimensional boundary
! conditions for temperature and salinity are handled.
gotm's avatar
gotm committed
15 16
!
! !USES:
kbk's avatar
kbk committed
17
   use halo_zones, only : H_TAG,U_TAG,V_TAG
18
   use domain, only: imin,jmin,imax,jmax,kmax,H,az,au,av
kbk's avatar
kbk committed
19
   use domain, only: nsbv,NWB,NNB,NEB,NSB,bdy_index
gotm's avatar
gotm committed
20 21 22 23 24 25 26 27
   use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
   use variables_3d
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_bdy_3d, do_bdy_3d
kbk's avatar
kbk committed
28
   REALTYPE, public, allocatable       :: S_bdy(:,:),T_bdy(:,:)
bjb's avatar
bjb committed
29 30 31 32 33
   REALTYPE, public                    :: bdy3d_rlxucut=_ONE_/50
   REALTYPE, public                    :: bdy3d_rlxmax=_ONE_/4
   REALTYPE, public                    :: bdy3d_rlxmin=_ZERO_
   REALTYPE,         allocatable       :: bdyvertS(:), bdyvertT(:)
   REALTYPE,         allocatable       :: rlxcoef(:)
gotm's avatar
gotm committed
34 35 36
!
! !PRIVATE DATA MEMBERS:
!
kbk's avatar
kbk committed
37 38 39
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
40 41 42 43 44 45 46 47 48
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
49 50
! !IROUTINE: init_bdy_3d - initialising 3D boundary conditions
! \label{sec-init-bdy-3d}
gotm's avatar
gotm committed
51 52 53 54 55
!
! !INTERFACE:
   subroutine init_bdy_3d()
!
! !DESCRIPTION:
56 57 58
!
! Here, the necessary fields {\tt S\_bdy} and {\tt T\_bdy} for
! salinity and temperature, respectively, are allocated.
gotm's avatar
gotm committed
59 60 61 62 63 64 65 66 67 68 69
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
70
   integer                   :: rc,i,j,k,n
gotm's avatar
gotm committed
71 72 73 74 75 76 77 78 79 80
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_bdy_3d() # ',Ncall
#endif

   LEVEL2 'init_bdy_3d()'
81
   allocate(S_bdy(0:kmax,nsbv),stat=rc)
gotm's avatar
gotm committed
82 83
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (S_bdy)'

84
   allocate(T_bdy(0:kmax,nsbv),stat=rc)
gotm's avatar
gotm committed
85 86
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (T_bdy)'

bjb's avatar
bjb committed
87 88 89 90 91 92 93 94 95 96
   allocate(bdyvertS(0:kmax),stat=rc)
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (bdyvertS)'

   allocate(bdyvertT(0:kmax),stat=rc)
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (bdyvertT)'

   allocate(rlxcoef(0:kmax),stat=rc)
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (rlxcoef)'


gotm's avatar
gotm committed
97 98 99 100 101 102 103 104 105 106 107
#ifdef DEBUG
   write(debug,*) 'Leaving init_bdy_3d()'
   write(debug,*)
#endif
   return
   end subroutine init_bdy_3d
!EOC

!-----------------------------------------------------------------------
!BOP
!
108
! !IROUTINE:  do_bdy_3d  - updating 3D boundary conditions
hb's avatar
hb committed
109
! \label{sec-do-bdy-3d}
gotm's avatar
gotm committed
110 111 112 113 114
!
! !INTERFACE:
   subroutine do_bdy_3d(tag,field)
!
! !DESCRIPTION:
115 116 117 118
!
! Here, the boundary conditions for salinity and temperature are
! copied to the boundary points and relaxed to the near boundary points
! by means of the flow relaxation scheme by \cite{MARTINSENea87}.
gotm's avatar
gotm committed
119
!
bjb's avatar
bjb committed
120 121 122 123 124 125 126 127
! As an extention to the flow relaxation scheme, it is possible 
! to relax the boundary point values to the specified boundary 
! condition in time, thus giving more realistic situations 
! especially for outgoing flow conditions. This nudging is implemented
! to depend on the local (3D) current velocity perpendicular to 
! the boundary. For strong outflow, the boundary condition is turned 
! off, while for inflows it is given a high impact.
!
gotm's avatar
gotm committed
128 129 130 131
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
132
   integer, intent(in)                 :: tag
gotm's avatar
gotm committed
133 134
!
! !INPUT/OUTPUT PARAMETERS:
kbk's avatar
kbk committed
135
   REALTYPE, intent(inout)             :: field(I3DFIELD)
gotm's avatar
gotm committed
136 137 138 139
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
bjb's avatar
bjb committed
140
   integer                   :: i,j,k,l,n,ii,jj,kk
kbk's avatar
kbk committed
141
   REALTYPE                  :: sp(1:4),rat
bjb's avatar
bjb committed
142 143
   REALTYPE                  :: bdy3d_rlxumin
   REALTYPE                  :: wsum
gotm's avatar
gotm committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_bdy_3d() # ',Ncall
#endif

#if 0
   select case (tag)
      case (1)
      ! Lateral zero-gradient boundary condition (north & south)
         do k=1,kmax
158 159 160
            do i=imin,imax
               if (au(i,jmin) .eq. 3) field(i,jmin,k)=field(i,jmin+1,k)
               if (au(i,jmax) .eq. 3) field(i,jmax,k)=field(i,jmax-1,k)
gotm's avatar
gotm committed
161 162 163 164 165
            end do
         end do
      case (2)
      ! Lateral zero-gradient boundary conditions (west & east)
         do k=1,kmax
166 167 168
            do j=jmin,jmax
               if (av(imin,j) .eq. 3) field(imin,j,k)=field(imin+1,j,k)
               if (av(imax,j) .eq. 3) field(imax,j,k)=field(imax-1,j,k)
gotm's avatar
gotm committed
169 170 171 172
            end do
         end do
      case default
         FATAL 'Non valid tag'
kbk's avatar
kbk committed
173
         stop 'do_bdy_3d'
gotm's avatar
gotm committed
174 175 176
   end select
#endif

kbk's avatar
kbk committed
177
#ifndef NO_BAROCLINIC
gotm's avatar
gotm committed
178 179 180 181 182
! Sponge layer factors according to Martinsen and Engedahl, 1987.
   sp(1)=1.0
   sp(2)=0.5625
   sp(3)=0.25
   sp(4)=0.0625
bjb's avatar
bjb committed
183 184 185
! Hardcoding of lower limit of velocity cut-off for temporal relaxation.
! Linear variation between bdy3d_rlxumin and bdy3d_rlxucut.
   bdy3d_rlxumin = -0.25*bdy3d_rlxucut
gotm's avatar
gotm committed
186

kbk's avatar
kbk committed
187
   l = 0
gotm's avatar
gotm committed
188 189
   k = 0
   do n=1,NWB
kbk's avatar
kbk committed
190 191
      l = l+1
      k = bdy_index(l)
gotm's avatar
gotm committed
192
      i = wi(n)
gotm's avatar
gotm committed
193
      do j=wfj(n),wlj(n)
bjb's avatar
bjb committed
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
         if (au(i,j).gt.0) then
            ! Local temporal relaxation coeficient depends on 
            ! local current just *inside* domain:
            do kk=1,kmax
               if (uu(i,j,kk).ge.bdy3d_rlxucut) then
                  rlxcoef(kk) = bdy3d_rlxmax
               else if (uu(i,j,kk).le.bdy3d_rlxumin) then
                  rlxcoef(kk) = bdy3d_rlxmin
               else
                  rlxcoef(kk) = (bdy3d_rlxmax-bdy3d_rlxmin)*(uu(i,j,kk)-bdy3d_rlxumin) &
                       /(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
               end if
            end do
         else
            rlxcoef(:)=bdy3d_rlxmax
         end if
         ! Temporal relaxation: Weight inner (actual) solution near boundary 
         ! with boundary condition (outer solution.)
         wsum= MIN(az(i-1+2,j),1)*sp(2)+MIN(az(i-1+3,j),1)*sp(3)+MIN(az(i-1+4,j),1)*sp(4)
         if (wsum>_ZERO_) then
            ! Get (weighted avr of) inner near-bdy solution (sponge) cells:
            bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i-1+2,j),1) * sp(2) * S(i-1+2,j,:) &
                                        +MIN(az(i-1+3,j),1) * sp(3) * S(i-1+3,j,:) &
                                        +MIN(az(i-1+4,j),1) * sp(4) * S(i-1+4,j,:) )
            bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i-1+2,j),1) * sp(2) * T(i-1+2,j,:) &
                                        +MIN(az(i-1+3,j),1) * sp(3) * T(i-1+3,j,:) &
                                        +MIN(az(i-1+4,j),1) * sp(4) * T(i-1+4,j,:) )
            ! Weight inner and outer (bc) solutions for use 
            ! in spatial relaxation/sponge
            bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
            bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
         else
            ! No near-bdy points. Just clamp bdy temporally:
            bdyvertS(:) = S_bdy(:,k)
            bdyvertT(:) = T_bdy(:,k)
         end if
gotm's avatar
gotm committed
230
         do ii=1,4
bjb's avatar
bjb committed
231
            ! Spatial relaxation / sponge:
kbk's avatar
kbk committed
232
            if (az(i-1+ii,j).gt.0) then
bjb's avatar
bjb committed
233 234
               S(i-1+ii,j,:) = sp(ii)*bdyvertS(:)+(_ONE_-sp(ii))*S(i-1+ii,j,:)
               T(i-1+ii,j,:) = sp(ii)*bdyvertT(:)+(_ONE_-sp(ii))*T(i-1+ii,j,:)
kbk's avatar
kbk committed
235
            end if
236 237 238 239 240
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i+1,j,:)
            end if
#endif
gotm's avatar
gotm committed
241
         end do
kbk's avatar
kbk committed
242
         k = k+1
gotm's avatar
gotm committed
243 244
      end do
   end do
kbk's avatar
kbk committed
245

gotm's avatar
gotm committed
246
   do n = 1,NNB
kbk's avatar
kbk committed
247 248
      l = l+1
      k = bdy_index(l)
gotm's avatar
gotm committed
249 250
      j = nj(n)
      do i = nfi(n),nli(n)
bjb's avatar
bjb committed
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
         if (av(i,j-1).gt.0) then
            do kk=1,kmax
               if (vv(i,j-1,kk).le.-bdy3d_rlxucut) then
                  rlxcoef(kk) = bdy3d_rlxmax
               else if (vv(i,j-1,kk).ge.-bdy3d_rlxumin) then
                  rlxcoef(kk) = bdy3d_rlxmin
               else
                  rlxcoef(kk) = -(bdy3d_rlxmax-bdy3d_rlxmin)*(vv(i,j-1,kk)+bdy3d_rlxumin) &
                       /(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
               end if
            end do
         else
            rlxcoef(:)=bdy3d_rlxmax
         end if
         wsum= MIN(az(i,j+1-2),1)*sp(2)+MIN(az(i,j+1-3),1)*sp(3)+MIN(az(i,j+1-4),1)*sp(4)
         if (wsum>_ZERO_) then
            bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i,j+1-2),1) * sp(2) * S(i,j+1-2,:) &
                                        +MIN(az(i,j+1-3),1) * sp(3) * S(i,j+1-3,:) &
                                        +MIN(az(i,j+1-4),1) * sp(4) * S(i,j+1-4,:) )
            bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i,j+1-2),1) * sp(2) * T(i,j+1-2,:) &
                                        +MIN(az(i,j+1-3),1) * sp(3) * T(i,j+1-3,:) &
                                        +MIN(az(i,j+1-4),1) * sp(4) * T(i,j+1-4,:) )
            bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
            bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
         else
            bdyvertS(:) = S_bdy(:,k)
            bdyvertT(:) = T_bdy(:,k)
         end if
gotm's avatar
gotm committed
279
         do jj=1,4
kbk's avatar
kbk committed
280
            if (az(i,j+1-jj).gt.0) then
bjb's avatar
bjb committed
281 282
               S(i,j+1-jj,:) = sp(jj)*bdyvertS(:)+(_ONE_-sp(jj))*S(i,j+1-jj,:)
               T(i,j+1-jj,:) = sp(jj)*bdyvertT(:)+(_ONE_-sp(jj))*T(i,j+1-jj,:)
kbk's avatar
kbk committed
283
            end if
284 285 286 287 288
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i,j-1,:)
            end if
#endif
gotm's avatar
gotm committed
289
         end do
kbk's avatar
kbk committed
290
         k = k+1
gotm's avatar
gotm committed
291 292
      end do
   end do
kbk's avatar
kbk committed
293

gotm's avatar
gotm committed
294
   do n=1,NEB
kbk's avatar
kbk committed
295 296
      l = l+1
      k = bdy_index(l)
gotm's avatar
gotm committed
297
      i = ei(n)
gotm's avatar
gotm committed
298
      do j=efj(n),elj(n)
bjb's avatar
bjb committed
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
         if (au(i-1,j).gt.0) then
            do kk=1,kmax
               if (uu(i-1,j,kk).le.-bdy3d_rlxucut) then
                  rlxcoef(kk) = bdy3d_rlxmax
               else if (uu(i-1,j,kk).ge.-bdy3d_rlxumin) then
                  rlxcoef(kk) = bdy3d_rlxmin
               else
                  rlxcoef(kk) = -(bdy3d_rlxmax-bdy3d_rlxmin)*(uu(i-1,j,kk)+bdy3d_rlxumin) &
                       /(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
               end if
            end do
         else
            rlxcoef(:)=bdy3d_rlxmax
         end if
         wsum= MIN(az(i+1-2,j),1)*sp(2)+MIN(az(i+1-3,j),1)*sp(3)+MIN(az(i+1-4,j),1)*sp(4)
         if (wsum>_ZERO_) then
            bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i+1-2,j),1) * sp(2) * S(i+1-2,j,:) &
                                        +MIN(az(i+1-3,j),1) * sp(3) * S(i+1-3,j,:) &
                                        +MIN(az(i+1-4,j),1) * sp(4) * S(i+1-4,j,:) )
            bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i+1-2,j),1) * sp(2) * T(i+1-2,j,:) &
                                        +MIN(az(i+1-3,j),1) * sp(3) * T(i+1-3,j,:) &
                                        +MIN(az(i+1-4,j),1) * sp(4) * T(i+1-4,j,:) )
            bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
            bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
         else
            bdyvertS(:) = S_bdy(:,k)
            bdyvertT(:) = T_bdy(:,k)
         end if
gotm's avatar
gotm committed
327
         do ii=1,4
kbk's avatar
kbk committed
328
            if (az(i+1-ii,j).gt.0) then
bjb's avatar
bjb committed
329 330
               S(i+1-ii,j,:) = sp(ii)*bdyvertS(:)+(_ONE_-sp(ii))*S(i+1-ii,j,:)
               T(i+1-ii,j,:) = sp(ii)*bdyvertT(:)+(_ONE_-sp(ii))*T(i+1-ii,j,:)
kbk's avatar
kbk committed
331
            end if
332 333 334 335 336
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i-1,j,:)
            end if
#endif
gotm's avatar
gotm committed
337
         end do
kbk's avatar
kbk committed
338
         k = k+1
gotm's avatar
gotm committed
339 340
      end do
   end do
kbk's avatar
kbk committed
341

gotm's avatar
gotm committed
342
   do n = 1,NSB
kbk's avatar
kbk committed
343 344
      l = l+1
      k = bdy_index(l)
gotm's avatar
gotm committed
345 346
      j = sj(n)
      do i = sfi(n),sli(n)
bjb's avatar
bjb committed
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
         if (av(i,j).gt.0) then
            do kk=1,kmax
               if (vv(i,j,kk).ge.bdy3d_rlxucut) then
                  rlxcoef(kk) = bdy3d_rlxmax
               else if (vv(i,j,kk).le.bdy3d_rlxumin) then
                  rlxcoef(kk) = bdy3d_rlxmin
               else
                  rlxcoef(kk) = (bdy3d_rlxmax-bdy3d_rlxmin)*(vv(i,j,kk)-bdy3d_rlxumin) &
                       /(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
               end if
            end do
         else
            rlxcoef(:)=bdy3d_rlxmax
         end if
         wsum= MIN(az(i,j-1+2),1)*sp(2)+MIN(az(i,j-1+3),1)*sp(3)+MIN(az(i,j-1+4),1)*sp(4)
         if (wsum>_ZERO_) then
            bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i,j-1+2),1) * sp(2) * S(i,j-1+2,:) &
                                        +MIN(az(i,j-1+3),1) * sp(3) * S(i,j-1+3,:) &
                                        +MIN(az(i,j-1+4),1) * sp(4) * S(i,j-1+4,:) )
            bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i,j-1+2),1) * sp(2) * T(i,j-1+2,:) &
                                        +MIN(az(i,j-1+3),1) * sp(3) * T(i,j-1+3,:) &
                                        +MIN(az(i,j-1+4),1) * sp(4) * T(i,j-1+4,:) )
            bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
            bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
         else
            bdyvertS(:) = S_bdy(:,k)
            bdyvertT(:) = T_bdy(:,k)
         end if
gotm's avatar
gotm committed
375
         do jj=1,4
kbk's avatar
kbk committed
376
            if (az(i,j-1+jj).gt.0) then
bjb's avatar
bjb committed
377 378
               S(i,j-1+jj,:) = sp(jj)*bdyvertS(:)+(_ONE_-sp(jj))*S(i,j-1+jj,:)
               T(i,j-1+jj,:) = sp(jj)*bdyvertT(:)+(_ONE_-sp(jj))*T(i,j-1+jj,:)
kbk's avatar
kbk committed
379
            end if
380 381 382 383 384
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i,j+1,:)
            end if
#endif
gotm's avatar
gotm committed
385
         end do
kbk's avatar
kbk committed
386
         k = k+1
gotm's avatar
gotm committed
387 388
      end do
   end do
kbk's avatar
kbk committed
389

390 391 392 393 394 395 396 397 398
#ifdef GETM_BIO
   if ( allocated(cc3d) ) then
      do n=1,size(cc3d,1)
         call mirror_bdy_3d(cc3d(n,:,:,:),H_TAG)
      end do
   end if
#endif


kbk's avatar
kbk committed
399 400
#endif

gotm's avatar
gotm committed
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415
#ifdef DEBUG
   write(debug,*) 'leaving do_bdy_3d()'
   write(debug,*)
#endif
   return
   end subroutine do_bdy_3d
!EOC

!-----------------------------------------------------------------------

   end module bdy_3d

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Karsten Bolding and Hans Burchard               !
!-----------------------------------------------------------------------