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

Last change on this file since 11536 was 11536, checked in by smasson, 13 months ago

trunk: merge dev_r10984_HPC-13 into the trunk

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