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 branches/UKMO/dev_r5107_kara_mld/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/UKMO/dev_r5107_kara_mld/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 5244

Last change on this file since 5244 was 5244, checked in by davestorkey, 9 years ago

UKMO Kara MLD branch: remove svn keyword property and clear keywords.

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$
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      CALL lbc_bdy_lnk( spgu(:,:), 'T', 1., ib_bdy )
180      !
181      igrd = 2      ! Flather bc on u-velocity;
182      !             ! remember that flagu=-1 if normal velocity direction is outward
183      !             ! I think we should rather use after ssh ?
184      DO jb = 1, idx%nblenrim(igrd)
185         ii  = idx%nbi(jb,igrd)
186         ij  = idx%nbj(jb,igrd) 
187         flagu => idx%flagu(jb,igrd)
188         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
189         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
190         !
191         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
192
193         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
194         ! Use characteristics method instead
195         zflag = ABS(flagu)
196         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
197         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
198      END DO
199      !
200      igrd = 3      ! Flather bc on v-velocity
201      !             ! remember that flagv=-1 if normal velocity direction is outward
202      DO jb = 1, idx%nblenrim(igrd)
203         ii  = idx%nbi(jb,igrd)
204         ij  = idx%nbj(jb,igrd) 
205         flagv => idx%flagv(jb,igrd)
206         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
207         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
208         !
209         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
210
211         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
212         ! Use characteristics method instead
213         zflag = ABS(flagv)
214         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
215         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
216      END DO
217      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
218      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
219      !
220      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
221      !
222   END SUBROUTINE bdy_dyn2d_fla
223
224
225   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
226      !!----------------------------------------------------------------------
227      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
228      !!             
229      !!              - Apply Orlanski radiation condition adaptively:
230      !!                  - radiation plus weak nudging at outflow points
231      !!                  - no radiation and strong nudging at inflow points
232      !!
233      !!
234      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
235      !!----------------------------------------------------------------------
236      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
237      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
238      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
239      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
240      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
241      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
242
243      INTEGER  ::   ib, igrd                               ! dummy loop indices
244      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
245      !!----------------------------------------------------------------------
246
247      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
248      !
249      igrd = 2      ! Orlanski bc on u-velocity;
250      !           
251      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
252
253      igrd = 3      ! Orlanski bc on v-velocity
254     
255      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
256      !
257      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
258      !
259      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
260      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
261      !
262      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
263      !
264   END SUBROUTINE bdy_dyn2d_orlanski
265
266   SUBROUTINE bdy_ssh( zssh )
267      !!----------------------------------------------------------------------
268      !!                  ***  SUBROUTINE bdy_ssh  ***
269      !!
270      !! ** Purpose : Duplicate sea level across open boundaries
271      !!
272      !!----------------------------------------------------------------------
273      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
274      !!
275      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
276      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
277
278      igrd = 1                       ! Everything is at T-points here
279
280      DO ib_bdy = 1, nb_bdy
281         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
282            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
283            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
284            ! Set gradient direction:
285            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
286            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
287            IF ( zcoef1+zcoef2 == 0 ) THEN
288               ! corner
289!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
290!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
291!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
292!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
293!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
294               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
295               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
296                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
297                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
298                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
299               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
300            ELSE
301               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
302               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
303               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
304            ENDIF
305         END DO
306
307         ! Boundary points should be updated
308         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
309      END DO
310
311   END SUBROUTINE bdy_ssh
312
313#else
314   !!----------------------------------------------------------------------
315   !!   Dummy module                   NO Unstruct Open Boundary Conditions
316   !!----------------------------------------------------------------------
317CONTAINS
318   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
319      INTEGER, intent(in) :: kt
320      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt
321   END SUBROUTINE bdy_dyn2d
322
323#endif
324
325   !!======================================================================
326END MODULE bdydyn2d
327
Note: See TracBrowser for help on using the repository browser.