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 branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90 @ 4007

Last change on this file since 4007 was 4007, checked in by davestorkey, 11 years ago
  1. Bug fixes for flagu/flagv calculation in bdyini.F90.
  2. Introduce masking of derivatives in radiation velocity calculation in bdylib.F90
  3. Change relaxation term in Orlanski implementation to explicit timestepping in bdylib.F90.
  4. Remove bdyfmask (redundant).
File size: 14.4 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#if defined key_bdy 
10   !!----------------------------------------------------------------------
11   !!   'key_bdy' :                    Unstructured Open Boundary Condition
12   !!----------------------------------------------------------------------
13   !!   bdy_dyn3d        : apply open boundary conditions to baroclinic velocities
14   !!   bdy_dyn3d_frs    : apply Flow Relaxation Scheme
15   !!----------------------------------------------------------------------
16   USE timing          ! Timing
17   USE wrk_nemo        ! Memory Allocation
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE bdy_oce         ! ocean open boundary conditions
21   USE bdylib          ! for orlanski library routines
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE in_out_manager  !
24   Use phycst
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn
30   PUBLIC   bdy_dyn3d_dmp ! routine called by step
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE bdy_dyn3d( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  SUBROUTINE bdy_dyn3d  ***
44      !!
45      !! ** Purpose : - Apply open boundary conditions for baroclinic velocities
46      !!
47      !!----------------------------------------------------------------------
48      INTEGER, INTENT( in ) :: kt     ! Main time step counter
49      !!
50      INTEGER               :: ib_bdy ! loop index
51      !!
52
53      DO ib_bdy=1, nb_bdy
54
55         SELECT CASE( cn_dyn3d(ib_bdy) )
56         CASE('none')
57            CYCLE
58         CASE('frs')
59            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
60         CASE('specified')
61            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
62         CASE('zero')
63            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
64         CASE('orlanski')
65            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
66         CASE('orlanski_npo')
67            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
68         CASE DEFAULT
69            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
70         END SELECT
71      ENDDO
72
73   END SUBROUTINE bdy_dyn3d
74
75   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
76      !!----------------------------------------------------------------------
77      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
78      !!
79      !! ** Purpose : - Apply a specified value for baroclinic velocities
80      !!                at open boundaries.
81      !!
82      !!----------------------------------------------------------------------
83      INTEGER                     ::   kt
84      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
85      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
86      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
87      !!
88      INTEGER  ::   jb, jk         ! dummy loop indices
89      INTEGER  ::   ii, ij, igrd   ! local integers
90      REAL(wp) ::   zwgt           ! boundary weight
91      !!----------------------------------------------------------------------
92      !
93      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe')
94      !
95      igrd = 2                      ! Relaxation of zonal 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            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
101         END DO
102      END DO
103      !
104      igrd = 3                      ! Relaxation of meridional velocity
105      DO jb = 1, idx%nblenrim(igrd)
106         DO jk = 1, jpkm1
107            ii   = idx%nbi(jb,igrd)
108            ij   = idx%nbj(jb,igrd)
109            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
110         END DO
111      END DO
112      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated 
113      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
114      !
115      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
116
117      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
118
119   END SUBROUTINE bdy_dyn3d_spe
120
121   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
122      !!----------------------------------------------------------------------
123      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
124      !!
125      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
126      !!
127      !!----------------------------------------------------------------------
128      INTEGER                     ::   kt
129      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
130      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
131      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
132      !!
133      INTEGER  ::   ib, ik         ! dummy loop indices
134      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
135      REAL(wp) ::   zwgt           ! boundary weight
136      !!----------------------------------------------------------------------
137      !
138      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro')
139      !
140      igrd = 2                       ! Everything is at T-points here
141      DO ib = 1, idx%nblenrim(igrd)
142         ii = idx%nbi(ib,igrd)
143         ij = idx%nbj(ib,igrd)
144         DO ik = 1, jpkm1
145            ua(ii,ij,ik) = 0._wp
146         END DO
147      END DO
148
149      igrd = 3                       ! Everything is at T-points here
150      DO ib = 1, idx%nblenrim(igrd)
151         ii = idx%nbi(ib,igrd)
152         ij = idx%nbj(ib,igrd)
153         DO ik = 1, jpkm1
154            va(ii,ij,ik) = 0._wp
155         END DO
156      END DO
157      !
158      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
159      !
160      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
161
162      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
163
164   END SUBROUTINE bdy_dyn3d_zro
165
166   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
167      !!----------------------------------------------------------------------
168      !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
169      !!
170      !! ** Purpose : - Apply the Flow Relaxation Scheme for baroclinic velocities
171      !!                at open boundaries.
172      !!
173      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
174      !!               a three-dimensional baroclinic ocean model with realistic
175      !!               topography. Tellus, 365-382.
176      !!----------------------------------------------------------------------
177      INTEGER                     ::   kt
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  ::   jb, jk         ! dummy loop indices
183      INTEGER  ::   ii, ij, igrd   ! local integers
184      REAL(wp) ::   zwgt           ! boundary weight
185      !!----------------------------------------------------------------------
186      !
187      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_frs')
188      !
189      igrd = 2                      ! Relaxation of zonal velocity
190      DO jb = 1, idx%nblen(igrd)
191         DO jk = 1, jpkm1
192            ii   = idx%nbi(jb,igrd)
193            ij   = idx%nbj(jb,igrd)
194            zwgt = idx%nbw(jb,igrd)
195            ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk)
196         END DO
197      END DO
198      !
199      igrd = 3                      ! Relaxation of meridional velocity
200      DO jb = 1, idx%nblen(igrd)
201         DO jk = 1, jpkm1
202            ii   = idx%nbi(jb,igrd)
203            ij   = idx%nbj(jb,igrd)
204            zwgt = idx%nbw(jb,igrd)
205            va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk)
206         END DO
207      END DO
208      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
209      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
210      !
211      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
212
213      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_frs')
214
215   END SUBROUTINE bdy_dyn3d_frs
216
217   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo )
218      !!----------------------------------------------------------------------
219      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
220      !!             
221      !!              - Apply Orlanski radiation to baroclinic velocities.
222      !!              - Wrapper routine for bdy_orlanski_3d
223      !!
224      !!
225      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
226      !!----------------------------------------------------------------------
227      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
228      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
229      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
230      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
231
232      INTEGER  ::   jb, igrd                               ! dummy loop indices
233      !!----------------------------------------------------------------------
234
235      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_orlanski')
236      !
237      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.
238      !
239      igrd = 2      ! Orlanski bc on u-velocity;
240      !           
241      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo )
242
243      igrd = 3      ! Orlanski bc on v-velocity
244      ! 
245      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo )
246      !
247      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
248      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
249      !
250      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_orlanski')
251      !
252   END SUBROUTINE bdy_dyn3d_orlanski
253
254
255   SUBROUTINE bdy_dyn3d_dmp( kt )
256      !!----------------------------------------------------------------------
257      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
258      !!
259      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
260      !!
261      !!----------------------------------------------------------------------
262      INTEGER                     ::   kt
263      !!
264      INTEGER  ::   jb, jk         ! dummy loop indices
265      INTEGER  ::   ii, ij, igrd   ! local integers
266      REAL(wp) ::   zwgt           ! boundary weight
267      INTEGER  ::  ib_bdy          ! loop index
268      !!----------------------------------------------------------------------
269      !
270      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
271      !
272      !-------------------------------------------------------
273      ! Remove barotropic part from before velocity
274      !-------------------------------------------------------
275      CALL wrk_alloc(jpi,jpj,pub2d,pvb2d) 
276
277      pub2d(:,:) = 0.e0
278      pvb2d(:,:) = 0.e0
279
280      DO jk = 1, jpkm1
281#if defined key_vvl
282         pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
283         pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
284#else
285         pub2d(:,:) = pub2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
286         pvb2d(:,:) = pvb2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
287#endif
288      END DO
289
290      IF( lk_vvl ) THEN
291         pub2d(:,:) = pub2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
292         pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
293      ELSE
294         pub2d(:,:) = pvb2d(:,:) * hur(:,:)
295         pvb2d(:,:) = pub2d(:,:) * hvr(:,:)
296      ENDIF
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                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
307                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * 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                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
318                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk)
319               END DO
320            END DO
321         ENDIF
322      ENDDO
323      !
324      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 
325      !
326      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
327      !
328      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
329
330   END SUBROUTINE bdy_dyn3d_dmp
331
332#else
333   !!----------------------------------------------------------------------
334   !!   Dummy module                   NO Unstruct Open Boundary Conditions
335   !!----------------------------------------------------------------------
336CONTAINS
337   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine
338      WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt
339   END SUBROUTINE bdy_dyn3d
340
341   SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine
342      WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt
343   END SUBROUTINE bdy_dyn3d_dmp
344
345#endif
346
347   !!======================================================================
348END MODULE bdydyn3d
Note: See TracBrowser for help on using the repository browser.