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

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

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