source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn3d.F90 @ 10957

Last change on this file since 10957 was 10957, checked in by davestorkey, 2 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : BDY module (including some removal of redundant code). Tested in AMM12 and ORCA1 (not full SETTE test at this stage).

  • 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, Kbb, puu, pvv, Kaa )
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      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
44      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
45      !
46      INTEGER ::   ib_bdy   ! loop index
47      !!----------------------------------------------------------------------
48      !
49      DO ib_bdy=1, nb_bdy
50         !
51         SELECT CASE( cn_dyn3d(ib_bdy) )
52         CASE('none')        ;   CYCLE
53         CASE('frs' )        ;   CALL bdy_dyn3d_frs(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
54         CASE('specified')   ;   CALL bdy_dyn3d_spe(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
55         CASE('zero')        ;   CALL bdy_dyn3d_zro(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
56         CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
57         CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
58         CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad(         puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
59         CASE('neumann')     ;   CALL bdy_dyn3d_nmn(           puu, pvv, Kaa, idx_bdy(ib_bdy), ib_bdy )
60         CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
61         END SELECT
62      END DO
63      !
64   END SUBROUTINE bdy_dyn3d
65
66
67   SUBROUTINE bdy_dyn3d_spe( puu, pvv, Kaa, idx, dta, ib_bdy )
68      !!----------------------------------------------------------------------
69      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
70      !!
71      !! ** Purpose : - Apply a specified value for baroclinic velocities
72      !!                at open boundaries.
73      !!
74      !!----------------------------------------------------------------------
75      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
76      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
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      igrd = 2                      ! Relaxation of zonal velocity
87      DO jb = 1, idx%nblenrim(igrd)
88         DO jk = 1, jpkm1
89            ii   = idx%nbi(jb,igrd)
90            ij   = idx%nbj(jb,igrd)
91            puu(ii,ij,jk,Kaa) = dta%u3d(jb,jk) * umask(ii,ij,jk)
92         END DO
93      END DO
94      !
95      igrd = 3                      ! Relaxation of meridional velocity
96      DO jb = 1, idx%nblenrim(igrd)
97         DO jk = 1, jpkm1
98            ii   = idx%nbi(jb,igrd)
99            ij   = idx%nbj(jb,igrd)
100            pvv(ii,ij,jk,Kaa) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
101         END DO
102      END DO
103      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ! Boundary points should be updated 
104      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
105      !
106   END SUBROUTINE bdy_dyn3d_spe
107
108
109   SUBROUTINE bdy_dyn3d_zgrad( puu, pvv, Kaa, idx, dta, ib_bdy )
110      !!----------------------------------------------------------------------
111      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
112      !!
113      !! ** Purpose : - Enforce a zero gradient of normal velocity
114      !!
115      !!----------------------------------------------------------------------
116      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
117      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
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      igrd = 2                      ! Copying tangential velocity into bdy points
129      DO jb = 1, idx%nblenrim(igrd)
130         DO jk = 1, jpkm1
131            ii   = idx%nbi(jb,igrd)
132            ij   = idx%nbj(jb,igrd)
133            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )
134            puu(ii,ij,jk,Kaa) = puu(ii,ij,jk,Kaa) * REAL( 1 - fu) + ( puu(ii,ij+fu,jk,Kaa) * umask(ii,ij+fu,jk) &
135                        &+ puu(ii,ij-fu,jk,Kaa) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )
136         END DO
137      END DO
138      !
139      igrd = 3                      ! Copying tangential velocity into bdy points
140      DO jb = 1, idx%nblenrim(igrd)
141         DO jk = 1, jpkm1
142            ii   = idx%nbi(jb,igrd)
143            ij   = idx%nbj(jb,igrd)
144            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )
145            pvv(ii,ij,jk,Kaa) = pvv(ii,ij,jk,Kaa) * REAL( 1 - fv ) + ( pvv(ii+fv,ij,jk,Kaa) * vmask(ii+fv,ij,jk) &
146                        &+ pvv(ii-fv,ij,jk,Kaa) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )
147         END DO
148      END DO
149      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ! Boundary points should be updated 
150      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
151      !
152   END SUBROUTINE bdy_dyn3d_zgrad
153
154
155   SUBROUTINE bdy_dyn3d_zro( puu, pvv, Kaa, idx, dta, ib_bdy )
156      !!----------------------------------------------------------------------
157      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
158      !!
159      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
160      !!
161      !!----------------------------------------------------------------------
162      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
163      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
164      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
165      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
166      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
167      !
168      INTEGER  ::   ib, ik         ! dummy loop indices
169      INTEGER  ::   ii, ij, igrd   ! local integers
170      REAL(wp) ::   zwgt           ! boundary weight
171      !!----------------------------------------------------------------------
172      !
173      igrd = 2                       ! Everything is at T-points here
174      DO ib = 1, idx%nblenrim(igrd)
175         ii = idx%nbi(ib,igrd)
176         ij = idx%nbj(ib,igrd)
177         DO ik = 1, jpkm1
178            puu(ii,ij,ik,Kaa) = 0._wp
179         END DO
180      END DO
181
182      igrd = 3                       ! Everything is at T-points here
183      DO ib = 1, idx%nblenrim(igrd)
184         ii = idx%nbi(ib,igrd)
185         ij = idx%nbj(ib,igrd)
186         DO ik = 1, jpkm1
187            pvv(ii,ij,ik,Kaa) = 0._wp
188         END DO
189      END DO
190      !
191      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1.,ib_bdy )   ! Boundary points should be updated
192      !
193   END SUBROUTINE bdy_dyn3d_zro
194
195
196   SUBROUTINE bdy_dyn3d_frs( puu, pvv, Kaa, idx, dta, ib_bdy )
197      !!----------------------------------------------------------------------
198      !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
199      !!
200      !! ** Purpose : - Apply the Flow Relaxation Scheme for baroclinic velocities
201      !!                at open boundaries.
202      !!
203      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
204      !!               a three-dimensional baroclinic ocean model with realistic
205      !!               topography. Tellus, 365-382.
206      !!----------------------------------------------------------------------
207      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
208      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
209      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
210      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
211      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
212      !
213      INTEGER  ::   jb, jk         ! dummy loop indices
214      INTEGER  ::   ii, ij, igrd   ! local integers
215      REAL(wp) ::   zwgt           ! boundary weight
216      !!----------------------------------------------------------------------
217      !
218      igrd = 2                      ! Relaxation of zonal velocity
219      DO jb = 1, idx%nblen(igrd)
220         DO jk = 1, jpkm1
221            ii   = idx%nbi(jb,igrd)
222            ij   = idx%nbj(jb,igrd)
223            zwgt = idx%nbw(jb,igrd)
224            puu(ii,ij,jk,Kaa) = ( puu(ii,ij,jk,Kaa) + zwgt * ( dta%u3d(jb,jk) - puu(ii,ij,jk,Kaa) ) ) * umask(ii,ij,jk)
225         END DO
226      END DO
227      !
228      igrd = 3                      ! Relaxation of meridional velocity
229      DO jb = 1, idx%nblen(igrd)
230         DO jk = 1, jpkm1
231            ii   = idx%nbi(jb,igrd)
232            ij   = idx%nbj(jb,igrd)
233            zwgt = idx%nbw(jb,igrd)
234            pvv(ii,ij,jk,Kaa) = ( pvv(ii,ij,jk,Kaa) + zwgt * ( dta%v3d(jb,jk) - pvv(ii,ij,jk,Kaa) ) ) * vmask(ii,ij,jk)
235         END DO
236      END DO
237      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
238      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
239      !
240   END SUBROUTINE bdy_dyn3d_frs
241
242
243   SUBROUTINE bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx, dta, ib_bdy, ll_npo )
244      !!----------------------------------------------------------------------
245      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
246      !!             
247      !!              - Apply Orlanski radiation to baroclinic velocities.
248      !!              - Wrapper routine for bdy_orlanski_3d
249      !!
250      !!
251      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
252      !!----------------------------------------------------------------------
253      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
254      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
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      LOGICAL                             , INTENT( in    ) ::   ll_npo  ! switch for NPO version
259
260      INTEGER  ::   jb, igrd                               ! dummy loop indices
261      !!----------------------------------------------------------------------
262      !
263      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities.
264      !
265      igrd = 2      ! Orlanski bc on u-velocity;
266      !           
267      CALL bdy_orlanski_3d( idx, igrd, puu(:,:,:,Kbb), puu(:,:,:,Kaa), dta%u3d, ll_npo )
268
269      igrd = 3      ! Orlanski bc on v-velocity
270     
271      CALL bdy_orlanski_3d( idx, igrd, pvv(:,:,:,Kbb), pvv(:,:,:,Kaa), dta%v3d, ll_npo )
272      !
273      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
274      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
275      !
276   END SUBROUTINE bdy_dyn3d_orlanski
277
278
279   SUBROUTINE bdy_dyn3d_dmp( kt, Kbb, puu, pvv, Krhs )
280      !!----------------------------------------------------------------------
281      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
282      !!
283      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
284      !!
285      !!----------------------------------------------------------------------
286      INTEGER                             , INTENT( in    ) ::   kt        ! time step
287      INTEGER                             , INTENT( in    ) ::   Kbb, Krhs ! Time level indices
288      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities and trends (to be updated at open boundaries)
289      !
290      INTEGER  ::   jb, jk         ! dummy loop indices
291      INTEGER  ::   ib_bdy         ! loop index
292      INTEGER  ::   ii, ij, igrd   ! local integers
293      REAL(wp) ::   zwgt           ! boundary weight
294      !!----------------------------------------------------------------------
295      !
296      IF( ln_timing )   CALL timing_start('bdy_dyn3d_dmp')
297      !
298      DO ib_bdy=1, nb_bdy
299         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN
300            igrd = 2                      ! Relaxation of zonal velocity
301            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
302               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
303               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
304               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
305               DO jk = 1, jpkm1
306                  puu(ii,ij,jk,Krhs) = ( puu(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
307                                   puu(ii,ij,jk,Kbb) + uu_b(ii,ij,Kbb)) ) * umask(ii,ij,jk)
308               END DO
309            END DO
310            !
311            igrd = 3                      ! Relaxation of meridional velocity
312            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
313               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
314               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
315               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
316               DO jk = 1, jpkm1
317                  pvv(ii,ij,jk,Krhs) = ( pvv(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
318                                   pvv(ii,ij,jk,Kbb) + vv_b(ii,ij,Kbb)) ) * vmask(ii,ij,jk)
319               END DO
320            END DO
321         ENDIF
322      END DO
323      !
324      CALL lbc_lnk_multi( 'bdydyn3d', puu(:,:,:,Krhs), 'U', -1.,  pvv(:,:,:,Krhs), 'V', -1. )   ! Boundary points should be updated
325      !
326      IF( ln_timing )   CALL timing_stop('bdy_dyn3d_dmp')
327      !
328   END SUBROUTINE bdy_dyn3d_dmp
329
330
331   SUBROUTINE bdy_dyn3d_nmn( puu, pvv, Kaa, idx, ib_bdy )
332      !!----------------------------------------------------------------------
333      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
334      !!             
335      !!              - Apply Neumann condition to baroclinic velocities.
336      !!              - Wrapper routine for bdy_nmn
337      !!
338      !!
339      !!----------------------------------------------------------------------
340      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
341      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
342      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
343      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
344
345      INTEGER  ::   jb, igrd                               ! dummy loop indices
346      !!----------------------------------------------------------------------
347      !
348      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities.
349      !
350      igrd = 2      ! Neumann bc on u-velocity;
351      !           
352      CALL bdy_nmn( idx, igrd, puu(:,:,:,Kaa) )
353
354      igrd = 3      ! Neumann bc on v-velocity
355     
356      CALL bdy_nmn( idx, igrd, pvv(:,:,:,Kaa) )
357      !
358      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
359      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )
360      !
361   END SUBROUTINE bdy_dyn3d_nmn
362
363   !!======================================================================
364END MODULE bdydyn3d
Note: See TracBrowser for help on using the repository browser.