source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 7 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

File size: 14.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#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      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points
269      !!----------------------------------------------------------------------
270      !
271      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
272      !
273      !-------------------------------------------------------
274      ! Remove barotropic part from before velocity
275      !-------------------------------------------------------
276      CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 
277
278      pub2d(:,:) = 0.e0
279      pvb2d(:,:) = 0.e0
280
281      phur1(:,:) = 0.
282      phvr1(:,:) = 0.
283      DO jk = 1, jpkm1
284#if defined key_vvl
285         phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
286         phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
287         pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
288         pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
289#else
290         pub2d(:,:) = pub2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
291         pvb2d(:,:) = pvb2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
292#endif
293      END DO
294
295      IF( lk_vvl ) THEN
296         phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) )
297         phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) )
298         pub2d(:,:) = pub2d(:,:) * umask(:,:,1) * phur1(:,:)
299         pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) * phvr1(:,:)
300      ELSE
301         pub2d(:,:) = pvb2d(:,:) * hur(:,:)
302         pvb2d(:,:) = pub2d(:,:) * hvr(:,:)
303      ENDIF
304
305      DO ib_bdy=1, nb_bdy
306         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN
307            igrd = 2                      ! Relaxation of zonal velocity
308            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
309               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
310               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
311               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
312               DO jk = 1, jpkm1
313                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
314                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk)
315               END DO
316            END DO
317            !
318            igrd = 3                      ! Relaxation of meridional velocity
319            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
320               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
321               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
322               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
323               DO jk = 1, jpkm1
324                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
325                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk)
326               END DO
327            END DO
328         ENDIF
329      ENDDO
330      !
331      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 
332      !
333      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
334      !
335      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
336
337   END SUBROUTINE bdy_dyn3d_dmp
338
339#else
340   !!----------------------------------------------------------------------
341   !!   Dummy module                   NO Unstruct Open Boundary Conditions
342   !!----------------------------------------------------------------------
343CONTAINS
344   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine
345      WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt
346   END SUBROUTINE bdy_dyn3d
347
348   SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine
349      WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt
350   END SUBROUTINE bdy_dyn3d_dmp
351
352#endif
353
354   !!======================================================================
355END MODULE bdydyn3d
Note: See TracBrowser for help on using the repository browser.