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 @ 11067

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

dev_r10984_HPC-13 : new implementation of lbc_bdy_lnk in prevision of step 2, regroup communications, see #2285

  • Property svn:keywords set to Id
File size: 14.7 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      LOGICAL, DIMENSION(4) :: lsend2, lrecv2, lsend3, lrecv3  ! indicate how communications are to be carried out
54     
55      DO ib_bdy=1, nb_bdy
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      lsend2(:) = .false.
75      lrecv2(:) = .false.
76      lsend3(:) = .false.
77      lrecv3(:) = .false.
78      DO ib_bdy=1, nb_bdy
79         SELECT CASE( cn_dyn2d(ib_bdy) )
80         CASE('flather')
81            lsend2(:) = lsend2(:) .OR. lsend_bdy(ib_bdy,2,:)   ! to   every bdy neighbour, U points
82            lrecv2(:) = lrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:)   ! from every bdy neighbour, U points
83            lsend3(:) = lsend3(:) .OR. lsend_bdy(ib_bdy,3,:)   ! to   every bdy neighbour, V points
84            lrecv3(:) = lrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:)   ! from every bdy neighbour, V points
85         CASE('orlanski')
86            lsend2(:) = lsend2(:) .OR. lsend_bdy(ib_bdy,2,:)   ! to   every bdy neighbour, U points
87            lrecv2(:) = lrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:)   ! from every bdy neighbour, U points
88            lsend3(:) = lsend3(:) .OR. lsend_bdy(ib_bdy,3,:)   ! to   every bdy neighbour, V points
89            lrecv3(:) = lrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:)   ! from every bdy neighbour, V points
90         CASE('orlanski_npo')
91            lsend2(:) = lsend2(:) .OR. lsend_bdy(ib_bdy,2,:)   ! to   every bdy neighbour, U points
92            lrecv2(:) = lrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:)   ! from every bdy neighbour, U points
93            lsend3(:) = lsend3(:) .OR. lsend_bdy(ib_bdy,3,:)   ! to   every bdy neighbour, V points
94            lrecv3(:) = lrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:)   ! from every bdy neighbour, V points
95         END SELECT
96      END DO
97      IF( ANY(lsend2) .OR. ANY(lrecv2) ) THEN   ! if need to send/recv in at least one direction
98         CALL lbc_bdy_lnk( 'bdydyn2d', lsend2, lrecv2, pua2d, 'U', -1. )
99      END IF
100      IF( ANY(lsend3) .OR. ANY(lrecv3) ) THEN   ! if need to send/recv in at least one direction
101         CALL lbc_bdy_lnk( 'bdydyn2d', lsend3, lrecv3, pva2d, 'V', -1. )
102      END IF
103      !
104   END SUBROUTINE bdy_dyn2d
105
106   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
107      !!----------------------------------------------------------------------
108      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
109      !!
110      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
111      !!                at open boundaries.
112      !!
113      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
114      !!               a three-dimensional baroclinic ocean model with realistic
115      !!               topography. Tellus, 365-382.
116      !!----------------------------------------------------------------------
117      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
118      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
119      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
120      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
121      !!
122      INTEGER  ::   jb, jk         ! dummy loop indices
123      INTEGER  ::   ii, ij, igrd   ! local integers
124      REAL(wp) ::   zwgt           ! boundary weight
125      !!----------------------------------------------------------------------
126      !
127      igrd = 2                      ! Relaxation of zonal velocity
128      DO jb = 1, idx%nblen(igrd)
129         ii   = idx%nbi(jb,igrd)
130         ij   = idx%nbj(jb,igrd)
131         zwgt = idx%nbw(jb,igrd)
132         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
133      END DO
134      !
135      igrd = 3                      ! Relaxation of meridional velocity
136      DO jb = 1, idx%nblen(igrd)
137         ii   = idx%nbi(jb,igrd)
138         ij   = idx%nbj(jb,igrd)
139         zwgt = idx%nbw(jb,igrd)
140         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
141      END DO 
142      !
143   END SUBROUTINE bdy_dyn2d_frs
144
145
146   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )
147      !!----------------------------------------------------------------------
148      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
149      !!             
150      !!              - Apply Flather boundary conditions on normal barotropic velocities
151      !!
152      !! ** WARNINGS about FLATHER implementation:
153      !!1. According to Palma and Matano, 1998 "after ssh" is used.
154      !!   In ROMS and POM implementations, it is "now ssh". In the current
155      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
156      !!   So I use "before ssh" in the following.
157      !!
158      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
159      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
160      !!   ssh in the code is not updated).
161      !!
162      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
163      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
164      !!----------------------------------------------------------------------
165      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
166      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
167      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
168      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
169      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh, phur, phvr 
170
171      INTEGER  ::   jb, igrd                         ! dummy loop indices
172      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
173      REAL(wp), POINTER :: flagu, flagv              ! short cuts
174      REAL(wp) ::   zcorr                            ! Flather correction
175      REAL(wp) ::   zforc                            ! temporary scalar
176      REAL(wp) ::   zflag, z1_2                      !    "        "
177      !!----------------------------------------------------------------------
178
179      z1_2 = 0.5_wp
180
181      ! ---------------------------------!
182      ! Flather boundary conditions     :!
183      ! ---------------------------------!
184     
185!!! REPLACE spgu with nemo_wrk work space
186
187      ! Fill temporary array with ssh data (here spgu):
188      igrd = 1
189      spgu(:,:) = 0.0
190      DO jb = 1, idx%nblenrim(igrd)
191         ii = idx%nbi(jb,igrd)
192         ij = idx%nbj(jb,igrd)
193         IF( ll_wd ) THEN
194            spgu(ii, ij) = dta%ssh(jb)  - ssh_ref 
195         ELSE
196            spgu(ii, ij) = dta%ssh(jb)
197         ENDIF
198      END DO
199
200      !
201      igrd = 2      ! Flather bc on u-velocity;
202      !             ! remember that flagu=-1 if normal velocity direction is outward
203      !             ! I think we should rather use after ssh ?
204      DO jb = 1, idx%nblenrim(igrd)
205         ii  = idx%nbi(jb,igrd)
206         ij  = idx%nbj(jb,igrd)
207         flagu => idx%flagu(jb,igrd)
208         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
209         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
210         IF( iim1 > jpi .OR. iip1 > jpi )   CYCLE
211         !
212         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
213
214         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
215         ! Use characteristics method instead
216         zflag = ABS(flagu)
217         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
218         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
219      END DO
220      !
221      igrd = 3      ! Flather bc on v-velocity
222      !             ! remember that flagv=-1 if normal velocity direction is outward
223      DO jb = 1, idx%nblenrim(igrd)
224         ii  = idx%nbi(jb,igrd)
225         ij  = idx%nbj(jb,igrd)
226         flagv => idx%flagv(jb,igrd)
227         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
228         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
229         IF( ijm1 > jpj .OR. ijp1 > jpj )   CYCLE
230         !
231         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
232
233         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
234         ! Use characteristics method instead
235         zflag = ABS(flagv)
236         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
237         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
238      END DO
239      !
240   END SUBROUTINE bdy_dyn2d_fla
241
242
243   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
244      !!----------------------------------------------------------------------
245      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
246      !!             
247      !!              - Apply Orlanski radiation condition adaptively:
248      !!                  - radiation plus weak nudging at outflow points
249      !!                  - no radiation and strong nudging at inflow points
250      !!
251      !!
252      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
253      !!----------------------------------------------------------------------
254      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
255      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
256      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
257      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
258      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
259      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
260
261      INTEGER  ::   ib, igrd                               ! dummy loop indices
262      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
263      !!----------------------------------------------------------------------
264      !
265      igrd = 2      ! Orlanski bc on u-velocity;
266      !           
267      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
268
269      igrd = 3      ! Orlanski bc on v-velocity
270     
271      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
272      !
273   END SUBROUTINE bdy_dyn2d_orlanski
274
275
276   SUBROUTINE bdy_ssh( zssh )
277      !!----------------------------------------------------------------------
278      !!                  ***  SUBROUTINE bdy_ssh  ***
279      !!
280      !! ** Purpose : Duplicate sea level across open boundaries
281      !!
282      !!----------------------------------------------------------------------
283      REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) ::   zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
284      !!
285      INTEGER  ::   ib_bdy          ! bdy index
286      LOGICAL, DIMENSION(4) :: lsend1, lrecv1  ! indicate how communications are to be carried out
287      !!----------------------------------------------------------------------
288      lsend1(:) = .false.
289      lrecv1(:) = .false.
290      DO ib_bdy = 1, nb_bdy
291         CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh )   ! zssh is masked
292         lsend1(:) = lsend1(:) .OR. lsend_bdy(ib_bdy,1,:)   ! to   every bdy neighbour, T points
293         lrecv1(:) = lrecv1(:) .OR. lrecv_bdy(ib_bdy,1,:)   ! from every bdy neighbour, T points
294      END DO
295      IF( ANY(lsend1) .OR. ANY(lrecv1) ) THEN   ! if need to send/recv in at least one direction
296         CALL lbc_bdy_lnk( 'bdydyn2d', lsend1, lrecv1, zssh(:,:,1), 'T',  1. )
297      END IF
298      !
299   END SUBROUTINE bdy_ssh
300
301   !!======================================================================
302END MODULE bdydyn2d
303
Note: See TracBrowser for help on using the repository browser.