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

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

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