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_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdydyn3d.F90 @ 12065

Last change on this file since 12065 was 12065, checked in by smueller, 4 years ago

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