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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 6043

Last change on this file since 6043 was 6043, checked in by timgraham, 8 years ago

Merged head of trunk into branch

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