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

Last change on this file since 15354 was 15354, checked in by smasson, 12 months ago

trunk: BDY compliant with corner communications, see #2731

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