source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90 @ 3680

Last change on this file since 3680 was 3680, checked in by rblod, 8 years ago

First commit of the final branch for 2012 (future nemo_3_5), see ticket #1028

File size: 12.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 lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE in_out_manager  !
23   Use phycst
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn
29   PUBLIC   bdy_dyn3d_dmp ! routine called by step
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
35   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
36   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE bdy_dyn3d( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  SUBROUTINE bdy_dyn3d  ***
43      !!
44      !! ** Purpose : - Apply open boundary conditions for baroclinic velocities
45      !!
46      !!----------------------------------------------------------------------
47      INTEGER, INTENT( in ) :: kt     ! Main time step counter
48      !!
49      INTEGER               :: ib_bdy ! loop index
50      !!
51
52      DO ib_bdy=1, nb_bdy
53
54!!$         IF ( using Orlanski radiation conditions ) THEN
55!!$            CALL bdy_rad( kt,  bdyidx(ib_bdy) )
56!!$         ENDIF
57
58         SELECT CASE( nn_dyn3d(ib_bdy) )
59         CASE(jp_none)
60            CYCLE
61         CASE(jp_frs)
62            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
63         CASE(2)
64            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
65         CASE(3)
66            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
67         CASE DEFAULT
68            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
69         END SELECT
70      ENDDO
71
72   END SUBROUTINE bdy_dyn3d
73
74   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
75      !!----------------------------------------------------------------------
76      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
77      !!
78      !! ** Purpose : - Apply a specified value for baroclinic velocities
79      !!                at open boundaries.
80      !!
81      !!----------------------------------------------------------------------
82      INTEGER                     ::   kt
83      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
84      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
85      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
86      !!
87      INTEGER  ::   jb, jk         ! dummy loop indices
88      INTEGER  ::   ii, ij, igrd   ! local integers
89      REAL(wp) ::   zwgt           ! boundary weight
90      !!----------------------------------------------------------------------
91      !
92      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe')
93      !
94      igrd = 2                      ! Relaxation of zonal velocity
95      DO jb = 1, idx%nblenrim(igrd)
96         DO jk = 1, jpkm1
97            ii   = idx%nbi(jb,igrd)
98            ij   = idx%nbj(jb,igrd)
99            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
100         END DO
101      END DO
102      !
103      igrd = 3                      ! Relaxation of meridional velocity
104      DO jb = 1, idx%nblenrim(igrd)
105         DO jk = 1, jpkm1
106            ii   = idx%nbi(jb,igrd)
107            ij   = idx%nbj(jb,igrd)
108            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
109         END DO
110      END DO
111      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
112      !
113      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
114
115      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
116
117   END SUBROUTINE bdy_dyn3d_spe
118
119   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
120      !!----------------------------------------------------------------------
121      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
122      !!
123      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
124      !!
125      !!----------------------------------------------------------------------
126      INTEGER                     ::   kt
127      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
128      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
129      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
130      !!
131      INTEGER  ::   ib, ik         ! dummy loop indices
132      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
133      REAL(wp) ::   zwgt           ! boundary weight
134      !!----------------------------------------------------------------------
135      !
136      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro')
137      !
138      igrd = 2                       ! Everything is at T-points here
139      DO ib = 1, idx%nblenrim(igrd)
140         ii = idx%nbi(ib,igrd)
141         ij = idx%nbj(ib,igrd)
142         DO ik = 1, jpkm1
143            ua(ii,ij,ik) = 0._wp
144         END DO
145      END DO
146
147      igrd = 3                       ! Everything is at T-points here
148      DO ib = 1, idx%nblenrim(igrd)
149         ii = idx%nbi(ib,igrd)
150         ij = idx%nbj(ib,igrd)
151         DO ik = 1, jpkm1
152            va(ii,ij,ik) = 0._wp
153         END DO
154      END DO
155      !
156      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
157      !
158      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
159
160      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
161
162   END SUBROUTINE bdy_dyn3d_zro
163
164   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
165      !!----------------------------------------------------------------------
166      !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
167      !!
168      !! ** Purpose : - Apply the Flow Relaxation Scheme for baroclinic velocities
169      !!                at open boundaries.
170      !!
171      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
172      !!               a three-dimensional baroclinic ocean model with realistic
173      !!               topography. Tellus, 365-382.
174      !!----------------------------------------------------------------------
175      INTEGER                     ::   kt
176      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
177      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
178      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
179      !!
180      INTEGER  ::   jb, jk         ! dummy loop indices
181      INTEGER  ::   ii, ij, igrd   ! local integers
182      REAL(wp) ::   zwgt           ! boundary weight
183      !!----------------------------------------------------------------------
184      !
185      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_frs')
186      !
187      igrd = 2                      ! Relaxation of zonal velocity
188      DO jb = 1, idx%nblen(igrd)
189         DO jk = 1, jpkm1
190            ii   = idx%nbi(jb,igrd)
191            ij   = idx%nbj(jb,igrd)
192            zwgt = idx%nbw(jb,igrd)
193            ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk)
194         END DO
195      END DO
196      !
197      igrd = 3                      ! Relaxation of meridional velocity
198      DO jb = 1, idx%nblen(igrd)
199         DO jk = 1, jpkm1
200            ii   = idx%nbi(jb,igrd)
201            ij   = idx%nbj(jb,igrd)
202            zwgt = idx%nbw(jb,igrd)
203            va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk)
204         END DO
205      END DO
206      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
207      !
208      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
209
210      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_frs')
211
212   END SUBROUTINE bdy_dyn3d_frs
213
214   SUBROUTINE bdy_dyn3d_dmp( kt )
215      !!----------------------------------------------------------------------
216      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
217      !!
218      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
219      !!
220      !!----------------------------------------------------------------------
221      INTEGER                     ::   kt
222      !!
223      INTEGER  ::   jb, jk         ! dummy loop indices
224      INTEGER  ::   ii, ij, igrd   ! local integers
225      REAL(wp) ::   zwgt           ! boundary weight
226      INTEGER  ::  ib_bdy          ! loop index
227      !!----------------------------------------------------------------------
228      !
229      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
230      !
231      !-------------------------------------------------------
232      ! Remove barotropic part from before velocity
233      !-------------------------------------------------------
234      CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
235
236      pu2d(:,:) = 0.e0
237      pv2d(:,:) = 0.e0
238
239      DO jk = 1, jpkm1
240#if defined key_vvl
241         pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
242         pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
243#else
244         pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
245         pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
246#endif
247      END DO
248
249      IF( lk_vvl ) THEN
250         pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
251         pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
252      ELSE
253         pu2d(:,:) = pv2d(:,:) * hur(:,:)
254         pv2d(:,:) = pu2d(:,:) * hvr(:,:)
255      ENDIF
256
257      DO ib_bdy=1, nb_bdy
258         IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN
259            igrd = 2                      ! Relaxation of zonal velocity
260            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
261               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
262               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
263               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
264               DO jk = 1, jpkm1
265                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk)
266               END DO
267            END DO
268            !
269            igrd = 3                      ! Relaxation of meridional velocity
270            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
271               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
272               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
273               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
274               DO jk = 1, jpkm1
275                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk)
276               END DO
277            END DO
278         ENDIF
279      ENDDO
280      !
281      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
282      !
283      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
284      !
285      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
286
287   END SUBROUTINE bdy_dyn3d_dmp
288
289#else
290   !!----------------------------------------------------------------------
291   !!   Dummy module                   NO Unstruct Open Boundary Conditions
292   !!----------------------------------------------------------------------
293CONTAINS
294   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine
295      WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt
296   END SUBROUTINE bdy_dyn3d
297
298   SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine
299      WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt
300   END SUBROUTINE bdy_dyn3d_dmp
301
302#endif
303
304   !!======================================================================
305END MODULE bdydyn3d
Note: See TracBrowser for help on using the repository browser.