source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 4317

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

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

File size: 14.3 KB
Line 
1MODULE bdydyn2d
2   !!======================================================================
3   !!                       ***  MODULE  bdydyn  ***
4   !! Unstructured Open Boundary Cond. :   Apply boundary conditions to barotropic solution
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   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
9   !!----------------------------------------------------------------------
10#if defined key_bdy 
11   !!----------------------------------------------------------------------
12   !!   'key_bdy' :                    Unstructured Open Boundary Condition
13   !!----------------------------------------------------------------------
14   !!   bdy_dyn2d          : Apply open boundary conditions to barotropic variables.
15   !!   bdy_dyn2d_frs      : Apply Flow Relaxation Scheme
16   !!   bdy_dyn2d_fla      : Apply Flather condition
17   !!   bdy_dyn2d_orlanski : Orlanski Radiation
18   !!   bdy_ssh            : Duplicate sea level across open boundaries
19   !!----------------------------------------------------------------------
20   USE timing          ! Timing
21   USE oce             ! ocean dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE bdy_oce         ! ocean open boundary conditions
24   USE bdylib          ! BDY library routines
25   USE dynspg_oce      ! for barotropic variables
26   USE phycst          ! physical constants
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE in_out_manager  !
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn
34   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
39   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE bdy_dyn2d( kt )
44      !!----------------------------------------------------------------------
45      !!                  ***  SUBROUTINE bdy_dyn2d  ***
46      !!
47      !! ** Purpose : - Apply open boundary conditions for barotropic variables
48      !!
49      !!----------------------------------------------------------------------
50      INTEGER,                      INTENT(in) ::   kt   ! Main time step counter
51      !!
52      INTEGER                                  ::   ib_bdy ! Loop counter
53
54      DO ib_bdy=1, nb_bdy
55
56         SELECT CASE( cn_dyn2d(ib_bdy) )
57         CASE('none')
58            CYCLE
59         CASE('frs')
60            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
61         CASE('flather')
62            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
63         CASE('orlanski')
64            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
65         CASE('orlanski_npo')
66            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
67         CASE DEFAULT
68            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
69         END SELECT
70      ENDDO
71
72   END SUBROUTINE bdy_dyn2d
73
74   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy )
75      !!----------------------------------------------------------------------
76      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
77      !!
78      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
79      !!                at open boundaries.
80      !!
81      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
82      !!               a three-dimensional baroclinic ocean model with realistic
83      !!               topography. Tellus, 365-382.
84      !!----------------------------------------------------------------------
85      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
86      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
87      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
88      !!
89      INTEGER  ::   jb, jk         ! dummy loop indices
90      INTEGER  ::   ii, ij, igrd   ! local integers
91      REAL(wp) ::   zwgt           ! boundary weight
92      !!----------------------------------------------------------------------
93      !
94      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_frs')
95      !
96      igrd = 2                      ! Relaxation of zonal velocity
97      DO jb = 1, idx%nblen(igrd)
98         ii   = idx%nbi(jb,igrd)
99         ij   = idx%nbj(jb,igrd)
100         zwgt = idx%nbw(jb,igrd)
101         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
102      END DO
103      !
104      igrd = 3                      ! Relaxation of meridional velocity
105      DO jb = 1, idx%nblen(igrd)
106         ii   = idx%nbi(jb,igrd)
107         ij   = idx%nbj(jb,igrd)
108         zwgt = idx%nbw(jb,igrd)
109         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
110      END DO
111      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) 
112      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
113      !
114      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
115      !
116
117   END SUBROUTINE bdy_dyn2d_frs
118
119
120   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy )
121      !!----------------------------------------------------------------------
122      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
123      !!             
124      !!              - Apply Flather boundary conditions on normal barotropic velocities
125      !!
126      !! ** WARNINGS about FLATHER implementation:
127      !!1. According to Palma and Matano, 1998 "after ssh" is used.
128      !!   In ROMS and POM implementations, it is "now ssh". In the current
129      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
130      !!   So I use "before ssh" in the following.
131      !!
132      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
133      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
134      !!   ssh in the code is not updated).
135      !!
136      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
137      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
138      !!----------------------------------------------------------------------
139      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
140      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
141      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
142
143      INTEGER  ::   jb, igrd                         ! dummy loop indices
144      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
145      REAL(wp), POINTER :: flagu, flagv              ! short cuts
146      REAL(wp) ::   zcorr                            ! Flather correction
147      REAL(wp) ::   zforc                            ! temporary scalar
148      REAL(wp) ::   zflag, z1_2                      !    "        "
149      !!----------------------------------------------------------------------
150
151      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
152
153      z1_2 = 0.5_wp
154
155      ! ---------------------------------!
156      ! Flather boundary conditions     :!
157      ! ---------------------------------!
158     
159!!! REPLACE spgu with nemo_wrk work space
160
161      ! Fill temporary array with ssh data (here spgu):
162      igrd = 1
163      spgu(:,:) = 0.0
164      DO jb = 1, idx%nblenrim(igrd)
165         ii = idx%nbi(jb,igrd)
166         ij = idx%nbj(jb,igrd)
167         spgu(ii, ij) = dta%ssh(jb)
168      END DO
169      !
170      igrd = 2      ! Flather bc on u-velocity;
171      !             ! remember that flagu=-1 if normal velocity direction is outward
172      !             ! I think we should rather use after ssh ?
173      DO jb = 1, idx%nblenrim(igrd)
174         ii  = idx%nbi(jb,igrd)
175         ij  = idx%nbj(jb,igrd) 
176         flagu => idx%flagu(jb,igrd)
177         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
178         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
179         !
180         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
181
182         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
183         ! Use characteristics method instead
184         zflag = ABS(flagu)
185         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
186         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
187      END DO
188      !
189      igrd = 3      ! Flather bc on v-velocity
190      !             ! remember that flagv=-1 if normal velocity direction is outward
191      DO jb = 1, idx%nblenrim(igrd)
192         ii  = idx%nbi(jb,igrd)
193         ij  = idx%nbj(jb,igrd) 
194         flagv => idx%flagv(jb,igrd)
195         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
196         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
197         !
198         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
199
200         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
201         ! Use characteristics method instead
202         zflag = ABS(flagv)
203         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
204         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
205      END DO
206      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
207      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
208      !
209      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
210      !
211   END SUBROUTINE bdy_dyn2d_fla
212
213
214   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo )
215      !!----------------------------------------------------------------------
216      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
217      !!             
218      !!              - Apply Orlanski radiation condition adaptively:
219      !!                  - radiation plus weak nudging at outflow points
220      !!                  - no radiation and strong nudging at inflow points
221      !!
222      !!
223      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
224      !!----------------------------------------------------------------------
225      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
226      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
227      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
228      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
229
230      INTEGER  ::   ib, igrd                               ! dummy loop indices
231      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
232      !!----------------------------------------------------------------------
233
234      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
235      !
236      igrd = 2      ! Orlanski bc on u-velocity;
237      !           
238      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
239
240      igrd = 3      ! Orlanski bc on v-velocity
241     
242      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
243      !
244      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
245      !
246      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
247      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
248      !
249      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
250      !
251   END SUBROUTINE bdy_dyn2d_orlanski
252
253   SUBROUTINE bdy_ssh( zssh )
254      !!----------------------------------------------------------------------
255      !!                  ***  SUBROUTINE bdy_ssh  ***
256      !!
257      !! ** Purpose : Duplicate sea level across open boundaries
258      !!
259      !!----------------------------------------------------------------------
260      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
261      !!
262      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
263      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
264
265      igrd = 1                       ! Everything is at T-points here
266
267      DO ib_bdy = 1, nb_bdy
268         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
269            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
270            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
271            ! Set gradient direction:
272            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
273            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
274            IF ( zcoef1+zcoef2 == 0 ) THEN
275               ! corner
276!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
277!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
278!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
279!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
280!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
281               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
282               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
283                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
284                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
285                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
286               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
287            ELSE
288               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
289               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
290               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
291            ENDIF
292         END DO
293
294         ! Boundary points should be updated
295         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
296      END DO
297
298   END SUBROUTINE bdy_ssh
299
300#else
301   !!----------------------------------------------------------------------
302   !!   Dummy module                   NO Unstruct Open Boundary Conditions
303   !!----------------------------------------------------------------------
304CONTAINS
305   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
306      INTEGER, intent(in) :: kt
307      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt
308   END SUBROUTINE bdy_dyn2d
309
310#endif
311
312   !!======================================================================
313END MODULE bdydyn2d
314
Note: See TracBrowser for help on using the repository browser.