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

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

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