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

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

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