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.
bdydyn2d.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 4673

Last change on this file since 4673 was 4354, checked in by jchanut, 10 years ago

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