source: NEMO/releases/release-4.0/src/OCE/BDY/bdydyn2d.F90 @ 11073

Last change on this file since 11073 was 11073, checked in by girrmann, 18 months ago

v4 : fix bug in flather for bdy, see #2287

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