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

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

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