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

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

dev_r10984_HPC-13 : bug fix and cleaning of flather, see #2287 and #2285
  • Property svn:keywords set to Id
File size: 15.6 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, only : sshdta => spgu             ! 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                           ! 2D addresses
172      INTEGER  ::   iiTrim, ijTrim                   ! T pts i/j-indice on the rim
173      INTEGER  ::   iiToce, ijToce, iiUoce, ijVoce   ! T, U and V pts i/j-indice of the ocean next to the rim
174      REAL(wp), POINTER :: flagu, flagv              ! short cuts
175      REAL(wp) ::   zfla                             ! Flather correction
176      REAL(wp) ::   z1_2                             !
177      !!----------------------------------------------------------------------
178
179      z1_2 = 0.5_wp
180
181      ! ---------------------------------!
182      ! Flather boundary conditions     :!
183      ! ---------------------------------!
184
185      ! Fill temporary array with ssh data (here we use spgu with the alias sshdta):
186      igrd = 1
187      sshdta(:,:) = 0.0
188      DO jb = 1, idx%nblenrim(igrd)
189         ii = idx%nbi(jb,igrd)
190         ij = idx%nbj(jb,igrd)
191         IF( ll_wd ) THEN
192            sshdta(ii, ij) = dta%ssh(jb)  - ssh_ref 
193         ELSE
194            sshdta(ii, ij) = dta%ssh(jb)
195         ENDIF
196      END DO
197      !
198      igrd = 2      ! Flather bc on u-velocity;
199      !             ! remember that flagu=-1 if normal velocity direction is outward
200      !             ! I think we should rather use after ssh ?
201      DO jb = 1, idx%nblenrim(igrd)
202         ii = idx%nbi(jb,igrd)
203         ij = idx%nbj(jb,igrd)
204         flagu => idx%flagu(jb,igrd)
205         IF( flagu == 0. ) THEN
206            pua2d(ii,ij) = dta%u2d(jb)
207         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and U points
208            IF( flagu ==  1. ) THEN   ;   iiTrim = ii     ;   iiToce = ii+1   ;   iiUoce = ii+1   ;   ENDIF
209            IF( flagu == -1. ) THEN   ;   iiTrim = ii+1   ;   iiToce = ii     ;   iiUoce = ii-1   ;   ENDIF
210            !
211            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
212            IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 )   CYCLE   
213            !
214            zfla = pua2d(iiUoce,ij) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) )
215            !
216            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
217            ! mix Flather scheme with velocity of the ocean next to the rim
218            pua2d(ii,ij) =  z1_2 * ( dta%u2d(jb) + zfla )
219         END IF
220      END DO
221      !
222      igrd = 3      ! Flather bc on v-velocity
223      !             ! remember that flagv=-1 if normal velocity direction is outward
224      DO jb = 1, idx%nblenrim(igrd)
225         ii = idx%nbi(jb,igrd)
226         ij = idx%nbj(jb,igrd)
227         flagv => idx%flagv(jb,igrd)
228         IF( flagv == 0. ) THEN
229            pva2d(ii,ij) = dta%v2d(jb)
230         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and V points
231            IF( flagv ==  1. ) THEN   ;   ijTrim = ij     ;   ijToce = ij+1   ;   ijVoce = ij+1   ;   ENDIF
232            IF( flagv == -1. ) THEN   ;   ijTrim = ij+1   ;   ijToce = ij     ;   ijVoce = ij-1   ;   ENDIF
233            !
234            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
235            IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 )   CYCLE
236            !
237            zfla = pva2d(ii,ijVoce) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) )
238            !
239            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
240            ! mix Flather scheme with velocity of the ocean next to the rim
241            pva2d(ii,ij) =  z1_2 * ( dta%v2d(jb) + zfla )
242         END IF
243      END DO
244      !
245   END SUBROUTINE bdy_dyn2d_fla
246
247
248   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
249      !!----------------------------------------------------------------------
250      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
251      !!             
252      !!              - Apply Orlanski radiation condition adaptively:
253      !!                  - radiation plus weak nudging at outflow points
254      !!                  - no radiation and strong nudging at inflow points
255      !!
256      !!
257      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
258      !!----------------------------------------------------------------------
259      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
260      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
261      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
262      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
263      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
264      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
265
266      INTEGER  ::   ib, igrd                               ! dummy loop indices
267      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
268      !!----------------------------------------------------------------------
269      !
270      igrd = 2      ! Orlanski bc on u-velocity;
271      !           
272      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
273
274      igrd = 3      ! Orlanski bc on v-velocity
275      ! 
276      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
277      !
278   END SUBROUTINE bdy_dyn2d_orlanski
279
280
281   SUBROUTINE bdy_ssh( zssh )
282      !!----------------------------------------------------------------------
283      !!                  ***  SUBROUTINE bdy_ssh  ***
284      !!
285      !! ** Purpose : Duplicate sea level across open boundaries
286      !!
287      !!----------------------------------------------------------------------
288      REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) ::   zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
289      !!
290      INTEGER  ::   ib_bdy          ! bdy index
291      LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how communications are to be carried out
292      !!----------------------------------------------------------------------
293      llsend1(:) = .false.
294      llrecv1(:) = .false.
295      DO ib_bdy = 1, nb_bdy
296         CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh )   ! zssh is masked
297         llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:)   ! possibly every direction, T points
298         llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:)   ! possibly every direction, T points
299      END DO
300      IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
301         CALL lbc_bdy_lnk( 'bdydyn2d', llsend1, llrecv1, zssh(:,:,1), 'T',  1. )
302      END IF
303      !
304   END SUBROUTINE bdy_ssh
305
306   !!======================================================================
307END MODULE bdydyn2d
308
Note: See TracBrowser for help on using the repository browser.