New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bdydyn3d.F90 in branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90 @ 6807

Last change on this file since 6807 was 6807, checked in by jamesharle, 9 years ago

adding addtional boundary condtions to BDY: 1) zero-grad tangential velocities and 2) Neumann boundary condtion for BC velocities
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1MODULE bdydyn3d
2   !!======================================================================
3   !!                       ***  MODULE  bdydyn3d  ***
4   !! Unstructured Open Boundary Cond. :   Flow relaxation scheme on baroclinic velocities
5   !!======================================================================
6   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
7   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
8   !!----------------------------------------------------------------------
9#if defined key_bdy 
10   !!----------------------------------------------------------------------
11   !!   'key_bdy' :                    Unstructured Open Boundary Condition
12   !!----------------------------------------------------------------------
13   !!   bdy_dyn3d        : apply open boundary conditions to baroclinic velocities
14   !!   bdy_dyn3d_frs    : apply Flow Relaxation Scheme
15   !!----------------------------------------------------------------------
16   USE timing          ! Timing
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE bdy_oce         ! ocean open boundary conditions
20   USE bdylib          ! for orlanski library routines
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE in_out_manager  !
23   Use phycst
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn
29   PUBLIC   bdy_dyn3d_dmp ! routine called by step
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
35   !! $Id$
36   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE bdy_dyn3d( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  SUBROUTINE bdy_dyn3d  ***
43      !!
44      !! ** Purpose : - Apply open boundary conditions for baroclinic velocities
45      !!
46      !!----------------------------------------------------------------------
47      INTEGER, INTENT( in ) :: kt     ! Main time step counter
48      !!
49      INTEGER               :: ib_bdy ! loop index
50      !!
51
52      DO ib_bdy=1, nb_bdy
53
54         SELECT CASE( cn_dyn3d(ib_bdy) )
55         CASE('none')
56            CYCLE
57         CASE('frs')
58            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
59         CASE('specified')
60            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
61         CASE('zerograd')
62            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
63         CASE('zero')
64            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
65         CASE('neumann')
66            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy )
67         CASE('orlanski')
68            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
69         CASE('orlanski_npo')
70            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
71         CASE DEFAULT
72            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
73         END SELECT
74      ENDDO
75
76   END SUBROUTINE bdy_dyn3d
77
78   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
79      !!----------------------------------------------------------------------
80      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
81      !!
82      !! ** Purpose : - Apply a specified value for baroclinic velocities
83      !!                at open boundaries.
84      !!
85      !!----------------------------------------------------------------------
86      INTEGER                     ::   kt
87      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
88      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
89      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
90      !!
91      INTEGER  ::   jb, jk         ! dummy loop indices
92      INTEGER  ::   ii, ij, igrd   ! local integers
93      REAL(wp) ::   zwgt           ! boundary weight
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe')
97      !
98      igrd = 2                      ! Relaxation of zonal velocity
99      DO jb = 1, idx%nblenrim(igrd)
100         DO jk = 1, jpkm1
101            ii   = idx%nbi(jb,igrd)
102            ij   = idx%nbj(jb,igrd)
103            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
104         END DO
105      END DO
106      !
107      igrd = 3                      ! Relaxation of meridional velocity
108      DO jb = 1, idx%nblenrim(igrd)
109         DO jk = 1, jpkm1
110            ii   = idx%nbi(jb,igrd)
111            ij   = idx%nbj(jb,igrd)
112            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
113         END DO
114      END DO
115      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated 
116      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
117      !
118      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
119
120      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
121
122   END SUBROUTINE bdy_dyn3d_spe
123
124   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )
125      !!----------------------------------------------------------------------
126      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
127      !!
128      !! ** Purpose : - Enforce a zero gradient of normal velocity
129      !!
130      !!----------------------------------------------------------------------
131      INTEGER                     ::   kt
132      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
133      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
134      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
135      !!
136      INTEGER  ::   jb, jk         ! dummy loop indices
137      INTEGER  ::   ii, ij, igrd   ! local integers
138      REAL(wp) ::   zwgt           ! boundary weight
139      INTEGER  ::   fu, fv
140      !!----------------------------------------------------------------------
141      !
142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad')
143      !
144      igrd = 2                      ! Copying tangential velocity into bdy points
145      DO jb = 1, idx%nblenrim(igrd)
146         DO jk = 1, jpkm1
147            ii   = idx%nbi(jb,igrd)
148            ij   = idx%nbj(jb,igrd)
149            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )
150            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) &
151                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )
152         END DO
153      END DO
154      !
155      igrd = 3                      ! Copying tangential velocity into bdy points
156      DO jb = 1, idx%nblenrim(igrd)
157         DO jk = 1, jpkm1
158            ii   = idx%nbi(jb,igrd)
159            ij   = idx%nbj(jb,igrd)
160            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )
161            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) &
162                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )
163         END DO
164      END DO
165      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated 
166      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
167      !
168      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
169
170      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad')
171
172   END SUBROUTINE bdy_dyn3d_zgrad
173
174   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
175      !!----------------------------------------------------------------------
176      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
177      !!
178      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
179      !!
180      !!----------------------------------------------------------------------
181      INTEGER                     ::   kt
182      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
183      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
184      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
185      !!
186      INTEGER  ::   ib, ik         ! dummy loop indices
187      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
188      REAL(wp) ::   zwgt           ! boundary weight
189      !!----------------------------------------------------------------------
190      !
191      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro')
192      !
193      igrd = 2                       ! Everything is at T-points here
194      DO ib = 1, idx%nblenrim(igrd)
195         ii = idx%nbi(ib,igrd)
196         ij = idx%nbj(ib,igrd)
197         DO ik = 1, jpkm1
198            ua(ii,ij,ik) = 0._wp
199         END DO
200      END DO
201
202      igrd = 3                       ! Everything is at T-points here
203      DO ib = 1, idx%nblenrim(igrd)
204         ii = idx%nbi(ib,igrd)
205         ij = idx%nbj(ib,igrd)
206         DO ik = 1, jpkm1
207            va(ii,ij,ik) = 0._wp
208         END DO
209      END DO
210      !
211      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
212      !
213      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
214
215      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
216
217   END SUBROUTINE bdy_dyn3d_zro
218
219   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
220      !!----------------------------------------------------------------------
221      !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
222      !!
223      !! ** Purpose : - Apply the Flow Relaxation Scheme for baroclinic velocities
224      !!                at open boundaries.
225      !!
226      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
227      !!               a three-dimensional baroclinic ocean model with realistic
228      !!               topography. Tellus, 365-382.
229      !!----------------------------------------------------------------------
230      INTEGER                     ::   kt
231      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
232      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
233      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
234      !!
235      INTEGER  ::   jb, jk         ! dummy loop indices
236      INTEGER  ::   ii, ij, igrd   ! local integers
237      REAL(wp) ::   zwgt           ! boundary weight
238      !!----------------------------------------------------------------------
239      !
240      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_frs')
241      !
242      igrd = 2                      ! Relaxation of zonal velocity
243      DO jb = 1, idx%nblen(igrd)
244         DO jk = 1, jpkm1
245            ii   = idx%nbi(jb,igrd)
246            ij   = idx%nbj(jb,igrd)
247            zwgt = idx%nbw(jb,igrd)
248            ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk)
249         END DO
250      END DO
251      !
252      igrd = 3                      ! Relaxation of meridional velocity
253      DO jb = 1, idx%nblen(igrd)
254         DO jk = 1, jpkm1
255            ii   = idx%nbi(jb,igrd)
256            ij   = idx%nbj(jb,igrd)
257            zwgt = idx%nbw(jb,igrd)
258            va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk)
259         END DO
260      END DO
261      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
262      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
263      !
264      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
265
266      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_frs')
267
268   END SUBROUTINE bdy_dyn3d_frs
269
270   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo )
271      !!----------------------------------------------------------------------
272      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
273      !!             
274      !!              - Apply Orlanski radiation to baroclinic velocities.
275      !!              - Wrapper routine for bdy_orlanski_3d
276      !!
277      !!
278      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
279      !!----------------------------------------------------------------------
280      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
281      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
282      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
283      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
284
285      INTEGER  ::   jb, igrd                               ! dummy loop indices
286      !!----------------------------------------------------------------------
287
288      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_orlanski')
289      !
290      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.
291      !
292      igrd = 2      ! Orlanski bc on u-velocity;
293      !           
294      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo )
295
296      igrd = 3      ! Orlanski bc on v-velocity
297      ! 
298      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo )
299      !
300      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
301      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
302      !
303      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_orlanski')
304      !
305   END SUBROUTINE bdy_dyn3d_orlanski
306
307
308   SUBROUTINE bdy_dyn3d_dmp( kt )
309      !!----------------------------------------------------------------------
310      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
311      !!
312      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
313      !!
314      !!----------------------------------------------------------------------
315      INTEGER                     ::   kt
316      !!
317      INTEGER  ::   jb, jk         ! dummy loop indices
318      INTEGER  ::   ii, ij, igrd   ! local integers
319      REAL(wp) ::   zwgt           ! boundary weight
320      INTEGER  ::  ib_bdy          ! loop index
321      !!----------------------------------------------------------------------
322      !
323      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
324      !
325      !-------------------------------------------------------
326
327      DO ib_bdy=1, nb_bdy
328         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN
329            igrd = 2                      ! Relaxation of zonal velocity
330            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
331               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
332               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
333               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
334               DO jk = 1, jpkm1
335                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
336                                   ub(ii,ij,jk) + ub_b(ii,ij)) ) * umask(ii,ij,jk)
337               END DO
338            END DO
339            !
340            igrd = 3                      ! Relaxation of meridional velocity
341            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
342               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
343               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
344               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
345               DO jk = 1, jpkm1
346                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
347                                   vb(ii,ij,jk) + vb_b(ii,ij)) ) * vmask(ii,ij,jk)
348               END DO
349            END DO
350         ENDIF
351      ENDDO
352      !
353      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
354      !
355      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
356
357   END SUBROUTINE bdy_dyn3d_dmp
358
359   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )
360      !!----------------------------------------------------------------------
361      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
362      !!             
363      !!              - Apply Neumann condition to baroclinic velocities.
364      !!              - Wrapper routine for bdy_nmn
365      !!
366      !!
367      !!----------------------------------------------------------------------
368      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
369      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
370
371      INTEGER  ::   jb, igrd                               ! dummy loop indices
372      !!----------------------------------------------------------------------
373
374      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn')
375      !
376      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.
377      !
378      igrd = 2      ! Neumann bc on u-velocity;
379      !           
380      CALL bdy_nmn( idx, igrd, ua )
381
382      igrd = 3      ! Neumann bc on v-velocity
383      ! 
384      CALL bdy_nmn( idx, igrd, va )
385      !
386      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
387      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )
388      !
389      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn')
390      !
391   END SUBROUTINE bdy_dyn3d_nmn
392
393#else
394   !!----------------------------------------------------------------------
395   !!   Dummy module                   NO Unstruct Open Boundary Conditions
396   !!----------------------------------------------------------------------
397CONTAINS
398   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine
399      WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt
400   END SUBROUTINE bdy_dyn3d
401
402   SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine
403      WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt
404   END SUBROUTINE bdy_dyn3d_dmp
405
406#endif
407
408   !!======================================================================
409END MODULE bdydyn3d
Note: See TracBrowser for help on using the repository browser.