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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90 @ 11048

Last change on this file since 11048 was 11048, checked in by girrmann, 5 years ago

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