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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90 @ 11049

Last change on this file since 11049 was 11049, checked in by girrmann, 5 years ago

dev_r10984_HPC-13 : CYCLE instruction is not systematic anymore, computation is done on the halo whenever possible and overwritten by lbc_bdy instruction, see #2285
  • Property svn:keywords set to Id
File size: 12.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         IF( iim1 > jpi .OR. iip1 > jpi )   CYCLE
184         !
185         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
186
187         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
188         ! Use characteristics method instead
189         zflag = ABS(flagu)
190         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
191         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
192      END DO
193      !
194      igrd = 3      ! Flather bc on v-velocity
195      !             ! remember that flagv=-1 if normal velocity direction is outward
196      DO jb = 1, idx%nblenrim(igrd)
197         ii  = idx%nbi(jb,igrd)
198         ij  = idx%nbj(jb,igrd)
199         flagv => idx%flagv(jb,igrd)
200         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
201         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
202         IF( ijm1 > jpj .OR. ijp1 > jpj )   CYCLE
203         !
204         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
205
206         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
207         ! Use characteristics method instead
208         zflag = ABS(flagv)
209         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
210         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
211      END DO
212      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
213      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
214      !
215   END SUBROUTINE bdy_dyn2d_fla
216
217
218   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
219      !!----------------------------------------------------------------------
220      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
221      !!             
222      !!              - Apply Orlanski radiation condition adaptively:
223      !!                  - radiation plus weak nudging at outflow points
224      !!                  - no radiation and strong nudging at inflow points
225      !!
226      !!
227      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
228      !!----------------------------------------------------------------------
229      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
230      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
231      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
232      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
233      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
234      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
235
236      INTEGER  ::   ib, igrd                               ! dummy loop indices
237      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
238      !!----------------------------------------------------------------------
239      !
240      igrd = 2      ! Orlanski bc on u-velocity;
241      !           
242      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
243
244      igrd = 3      ! Orlanski bc on v-velocity
245      ! 
246      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
247      !
248      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
249      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
250      !
251   END SUBROUTINE bdy_dyn2d_orlanski
252
253
254   SUBROUTINE bdy_ssh( zssh )
255      !!----------------------------------------------------------------------
256      !!                  ***  SUBROUTINE bdy_ssh  ***
257      !!
258      !! ** Purpose : Duplicate sea level across open boundaries
259      !!
260      !!----------------------------------------------------------------------
261      REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) ::   zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
262      !!
263      INTEGER  ::   ib_bdy          ! bdy index
264      !!----------------------------------------------------------------------
265      DO ib_bdy = 1, nb_bdy
266         CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh )   ! zssh is masked
267         CALL lbc_bdy_lnk( 'bdydyn2d', zssh(:,:,1), 'T', 1., ib_bdy )
268      END DO
269      !
270   END SUBROUTINE bdy_ssh
271
272   !!======================================================================
273END MODULE bdydyn2d
274
Note: See TracBrowser for help on using the repository browser.