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 NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90 @ 11350

Last change on this file since 11350 was 11350, checked in by jcastill, 5 years ago

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