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 NEMO/branches/UKMO/dev_r9950_GO6_mixing/src/OCE/BDY – NEMO

source: NEMO/branches/UKMO/dev_r9950_GO6_mixing/src/OCE/BDY/bdydyn2d.F90 @ 10323

Last change on this file since 10323 was 10323, checked in by davestorkey, 5 years ago

UKMO/dev_r9950_GO6_mixing: Update to be relative to rev 10321 of NEMO4_beta_mirror branch.

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