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/trunk/src/OCE/BDY – NEMO

source: NEMO/trunk/src/OCE/BDY/bdydyn2d.F90 @ 15360

Last change on this file since 15360 was 15360, checked in by smasson, 3 years ago

trunk: bugfix following [15354], see #2731
  • Property svn:keywords set to Id
File size: 17.5 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 dom_oce         ! ocean space and time domain
17   USE bdy_oce         ! ocean open boundary conditions
18   USE bdylib          ! BDY library routines
19   USE phycst          ! physical constants
20   USE lib_mpp, ONLY: jpfillnothing
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
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, ir     ! BDY set index, rim index
53      LOGICAL               ::   llrim0         ! indicate if rim 0 is treated
54      LOGICAL, DIMENSION(8) ::   llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out
55      !!----------------------------------------------------------------------
56     
57      llsend2(:) = .false.   ;   llrecv2(:) = .false.
58      llsend3(:) = .false.   ;   llrecv3(:) = .false.
59      DO ir = 1, 0, -1   ! treat rim 1 before rim 0
60         IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
61         ELSE                 ;   llrim0 = .FALSE.
62         END IF
63         DO ib_bdy=1, nb_bdy
64            SELECT CASE( cn_dyn2d(ib_bdy) )
65            CASE('none')
66               CYCLE
67            CASE('frs')   ! treat the whole boundary at once
68               IF( llrim0 )   CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d )
69            CASE('flather')
70               CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 )
71            CASE('orlanski')
72               CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
73                    & pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.false. )
74            CASE('orlanski_npo')
75               CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
76                    & pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.true.  )
77            CASE DEFAULT
78               CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
79            END SELECT
80         ENDDO
81         !
82         IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
83         IF( nn_hls == 1 ) THEN
84            llsend2(:) = .false.   ;   llrecv2(:) = .false.
85            llsend3(:) = .false.   ;   llrecv3(:) = .false.
86         END IF
87         DO ib_bdy=1, nb_bdy
88            SELECT CASE( cn_dyn2d(ib_bdy) )
89            CASE('flather')
90               llsend2(:) = llsend2(:) .OR. lsend_bdyext(ib_bdy,2,:,ir) .OR. lsend_bdyint(ib_bdy,2,:,ir)
91               llrecv2(:) = llrecv2(:) .OR. lrecv_bdyext(ib_bdy,2,:,ir) .OR. lrecv_bdyint(ib_bdy,2,:,ir)
92               llsend3(:) = llsend3(:) .OR. lsend_bdyext(ib_bdy,3,:,ir) .OR. lsend_bdyint(ib_bdy,3,:,ir)
93               llrecv3(:) = llrecv3(:) .OR. lrecv_bdyext(ib_bdy,3,:,ir) .OR. lrecv_bdyint(ib_bdy,3,:,ir)
94            CASE('orlanski', 'orlanski_npo')
95               llsend2(:) = llsend2(:) .OR. lsend_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points
96               llrecv2(:) = llrecv2(:) .OR. lrecv_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points
97               llsend3(:) = llsend3(:) .OR. lsend_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points
98               llrecv3(:) = llrecv3(:) .OR. lrecv_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points
99            END SELECT
100         END DO
101         IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
102            CALL lbc_lnk( 'bdydyn2d', pua2d, 'U', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
103         END IF
104         IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
105            CALL lbc_lnk( 'bdydyn2d', pva2d, 'V', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
106         END IF
107         !
108      END DO   ! ir
109      !
110   END SUBROUTINE bdy_dyn2d
111
112   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
113      !!----------------------------------------------------------------------
114      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
115      !!
116      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
117      !!                at open boundaries.
118      !!
119      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
120      !!               a three-dimensional baroclinic ocean model with realistic
121      !!               topography. Tellus, 365-382.
122      !!----------------------------------------------------------------------
123      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
124      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
125      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
126      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
127      !!
128      INTEGER  ::   jb             ! dummy loop indices
129      INTEGER  ::   ii, ij, igrd   ! local integers
130      REAL(wp) ::   zwgt           ! boundary weight
131      !!----------------------------------------------------------------------
132      !
133      igrd = 2                      ! Relaxation of zonal velocity
134      DO jb = 1, idx%nblen(igrd)
135         ii   = idx%nbi(jb,igrd)
136         ij   = idx%nbj(jb,igrd)
137         zwgt = idx%nbw(jb,igrd)
138         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
139      END DO
140      !
141      igrd = 3                      ! Relaxation of meridional velocity
142      DO jb = 1, idx%nblen(igrd)
143         ii   = idx%nbi(jb,igrd)
144         ij   = idx%nbj(jb,igrd)
145         zwgt = idx%nbw(jb,igrd)
146         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
147      END DO 
148      !
149   END SUBROUTINE bdy_dyn2d_frs
150
151
152   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 )
153      !!----------------------------------------------------------------------
154      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
155      !!             
156      !!              - Apply Flather boundary conditions on normal barotropic velocities
157      !!
158      !! ** WARNINGS about FLATHER implementation:
159      !!1. According to Palma and Matano, 1998 "after ssh" is used.
160      !!   In ROMS and POM implementations, it is "now ssh". In the current
161      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
162      !!   So I use "before ssh" in the following.
163      !!
164      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
165      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
166      !!   ssh in the code is not updated).
167      !!
168      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
169      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
170      !!----------------------------------------------------------------------
171      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
172      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
173      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
174      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
175      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh, phur, phvr
176      LOGICAL                     , INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
177      INTEGER  ::   ibeg, iend                       ! length of rim to be treated (rim 0 or rim 1)
178      INTEGER  ::   jb, igrd                         ! dummy loop indices
179      INTEGER  ::   ii, ij                           ! 2D addresses
180      INTEGER  ::   iiTrim, ijTrim                   ! T pts i/j-indice on the rim
181      INTEGER  ::   iiToce, ijToce, iiUoce, ijVoce   ! T, U and V pts i/j-indice of the ocean next to the rim
182      REAL(wp) ::   flagu, flagv                     ! short cuts
183      REAL(wp) ::   zfla                             ! Flather correction
184      REAL(wp) ::   z1_2                             !
185      REAL(wp), DIMENSION(jpi,jpj) ::   sshdta       ! 2D version of dta%ssh
186      !!----------------------------------------------------------------------
187
188      z1_2 = 0.5_wp
189
190      ! ---------------------------------!
191      ! Flather boundary conditions     :!
192      ! ---------------------------------!
193
194      ! Fill temporary array with ssh data (here we use spgu with the alias sshdta):
195      igrd = 1
196      IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)
197      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
198      END IF
199      !
200      DO jb = ibeg, iend
201         ii = idx%nbi(jb,igrd)
202         ij = idx%nbj(jb,igrd)
203         IF( ll_wd ) THEN   ;   sshdta(ii, ij) = dta%ssh(jb) - ssh_ref 
204         ELSE               ;   sshdta(ii, ij) = dta%ssh(jb)
205         ENDIF
206      END DO
207      !
208      igrd = 2      ! Flather bc on u-velocity
209      !             ! remember that flagu=-1 if normal velocity direction is outward
210      !             ! I think we should rather use after ssh ?
211      IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)
212      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
213      END IF
214      DO jb = ibeg, iend
215         ii    = idx%nbi(jb,igrd)
216         ij    = idx%nbj(jb,igrd)
217         flagu = idx%flagu(jb,igrd)
218         IF( flagu == 0. ) THEN
219            pua2d(ii,ij) = dta%u2d(jb)
220         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and U points
221            IF( flagu ==  1. ) THEN   ;   iiTrim = ii     ;   iiToce = ii+1   ;   iiUoce = ii+1   ;   ENDIF
222            IF( flagu == -1. ) THEN   ;   iiTrim = ii+1   ;   iiToce = ii     ;   iiUoce = ii-1   ;   ENDIF
223            !
224            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
225            IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 )   CYCLE   
226            !
227            zfla = dta%u2d(jb) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) )
228            !
229            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
230            ! mix Flather scheme with velocity of the ocean next to the rim
231            pua2d(ii,ij) =  z1_2 * ( pua2d(iiUoce,ij) + zfla )
232         END IF
233      END DO
234      !
235      igrd = 3      ! Flather bc on v-velocity
236      !             ! remember that flagv=-1 if normal velocity direction is outward
237      IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)
238      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
239      END IF
240      DO jb = ibeg, iend
241         ii    = idx%nbi(jb,igrd)
242         ij    = idx%nbj(jb,igrd)
243         flagv = idx%flagv(jb,igrd)
244         IF( flagv == 0. ) THEN
245            pva2d(ii,ij) = dta%v2d(jb)
246         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and V points
247            IF( flagv ==  1. ) THEN   ;   ijTrim = ij     ;   ijToce = ij+1   ;   ijVoce = ij+1   ;   ENDIF
248            IF( flagv == -1. ) THEN   ;   ijTrim = ij+1   ;   ijToce = ij     ;   ijVoce = ij-1   ;   ENDIF
249            !
250            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
251            IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 )   CYCLE
252            !
253            zfla = dta%v2d(jb) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) )
254            !
255            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
256            ! mix Flather scheme with velocity of the ocean next to the rim
257            pva2d(ii,ij) =  z1_2 * ( pva2d(ii,ijVoce) + zfla )
258         END IF
259      END DO
260      !
261   END SUBROUTINE bdy_dyn2d_fla
262
263
264   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo )
265      !!----------------------------------------------------------------------
266      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
267      !!             
268      !!              - Apply Orlanski radiation condition adaptively:
269      !!                  - radiation plus weak nudging at outflow points
270      !!                  - no radiation and strong nudging at inflow points
271      !!
272      !!
273      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
274      !!----------------------------------------------------------------------
275      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
276      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
277      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
278      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
279      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
280      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
281      LOGICAL,                      INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
282      INTEGER  ::   ib, igrd                               ! dummy loop indices
283      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
284      !!----------------------------------------------------------------------
285      !
286      igrd = 2      ! Orlanski bc on u-velocity;
287      !           
288      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, llrim0, ll_npo )
289
290      igrd = 3      ! Orlanski bc on v-velocity
291      ! 
292      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, llrim0, ll_npo )
293      !
294   END SUBROUTINE bdy_dyn2d_orlanski
295
296
297   SUBROUTINE bdy_ssh( zssh )
298      !!----------------------------------------------------------------------
299      !!                  ***  SUBROUTINE bdy_ssh  ***
300      !!
301      !! ** Purpose : Duplicate sea level across open boundaries
302      !!
303      !!----------------------------------------------------------------------
304      REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) ::   zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
305      !!
306      INTEGER ::   ib_bdy, ir      ! bdy index, rim index
307      INTEGER ::   ibeg, iend      ! length of rim to be treated (rim 0 or rim 1)
308      LOGICAL ::   llrim0          ! indicate if rim 0 is treated
309      LOGICAL, DIMENSION(8) :: llsend1, llrecv1  ! indicate how communications are to be carried out
310      !!----------------------------------------------------------------------
311      llsend1(:) = .false.   ;   llrecv1(:) = .false.
312      DO ir = 1, 0, -1   ! treat rim 1 before rim 0
313         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   END IF
314         IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
315         ELSE                 ;   llrim0 = .FALSE.
316         END IF
317         DO ib_bdy = 1, nb_bdy
318            CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh, llrim0 )   ! zssh is masked
319            llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points
320            llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points
321         END DO
322         IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
323         IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
324            CALL lbc_lnk( 'bdydyn2d', zssh(:,:,1), 'T',  1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
325         END IF
326      END DO
327      !
328   END SUBROUTINE bdy_ssh
329
330   !!======================================================================
331END MODULE bdydyn2d
332
Note: See TracBrowser for help on using the repository browser.