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 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

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