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

Last change on this file since 10518 was 10518, checked in by clem, 5 years ago

solve ticket #1997 as proposed
  • Property svn:keywords set to Id
File size: 13.8 KB
RevLine 
[3117]1MODULE bdydyn2d
2   !!======================================================================
3   !!                       ***  MODULE  bdydyn  ***
[3191]4   !! Unstructured Open Boundary Cond. :   Apply boundary conditions to barotropic solution
[3117]5   !!======================================================================
[3191]6   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
[3680]7   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
[4292]8   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
[3117]9   !!----------------------------------------------------------------------
[4292]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
[3117]15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE bdy_oce         ! ocean open boundary conditions
[4292]19   USE bdylib          ! BDY library routines
[3117]20   USE phycst          ! physical constants
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[9023]22   USE wet_dry         ! Use wet dry to get reference ssh level
[3117]23   USE in_out_manager  !
24
25   IMPLICIT NONE
26   PRIVATE
27
[4292]28   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn
29   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
[3117]30
31   !!----------------------------------------------------------------------
[9598]32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]33   !! $Id$
[10068]34   !! Software governed by the CeCILL license (see ./LICENSE)
[3117]35   !!----------------------------------------------------------------------
36CONTAINS
37
[4354]38   SUBROUTINE bdy_dyn2d( kt, pua2d, pva2d, pub2d, pvb2d, phur, phvr, pssh  )
[3117]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
[4354]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
[3117]50      !!
51      INTEGER                                  ::   ib_bdy ! Loop counter
52
53      DO ib_bdy=1, nb_bdy
54
[4292]55         SELECT CASE( cn_dyn2d(ib_bdy) )
56         CASE('none')
[3117]57            CYCLE
[4292]58         CASE('frs')
[4354]59            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d )
[4292]60         CASE('flather')
[4354]61            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr )
[4292]62         CASE('orlanski')
[4354]63            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
64                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.false.)
[4292]65         CASE('orlanski_npo')
[4354]66            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
67                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.true. )
[3117]68         CASE DEFAULT
[3191]69            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
[3117]70         END SELECT
71      ENDDO
72
73   END SUBROUTINE bdy_dyn2d
74
[4354]75   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
[3117]76      !!----------------------------------------------------------------------
77      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
78      !!
79      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
80      !!                at open boundaries.
81      !!
82      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
83      !!               a three-dimensional baroclinic ocean model with realistic
84      !!               topography. Tellus, 365-382.
85      !!----------------------------------------------------------------------
86      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
87      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
[3680]88      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
[4354]89      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
[3117]90      !!
91      INTEGER  ::   jb, jk         ! dummy loop indices
92      INTEGER  ::   ii, ij, igrd   ! local integers
93      REAL(wp) ::   zwgt           ! boundary weight
94      !!----------------------------------------------------------------------
95      !
96      igrd = 2                      ! Relaxation of zonal velocity
97      DO jb = 1, idx%nblen(igrd)
98         ii   = idx%nbi(jb,igrd)
99         ij   = idx%nbj(jb,igrd)
100         zwgt = idx%nbw(jb,igrd)
[4292]101         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
[3117]102      END DO
103      !
104      igrd = 3                      ! Relaxation of meridional velocity
105      DO jb = 1, idx%nblen(igrd)
106         ii   = idx%nbi(jb,igrd)
107         ij   = idx%nbj(jb,igrd)
108         zwgt = idx%nbw(jb,igrd)
[4292]109         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
[3117]110      END DO
[10425]111      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy ) 
112      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
[3117]113      !
114   END SUBROUTINE bdy_dyn2d_frs
115
116
[4354]117   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )
[3117]118      !!----------------------------------------------------------------------
119      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
120      !!             
121      !!              - Apply Flather boundary conditions on normal barotropic velocities
122      !!
123      !! ** WARNINGS about FLATHER implementation:
124      !!1. According to Palma and Matano, 1998 "after ssh" is used.
125      !!   In ROMS and POM implementations, it is "now ssh". In the current
126      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
127      !!   So I use "before ssh" in the following.
128      !!
129      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
130      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
131      !!   ssh in the code is not updated).
132      !!
133      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
134      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
135      !!----------------------------------------------------------------------
136      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
137      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
[3680]138      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
[4354]139      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
140      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh, phur, phvr 
[3117]141
142      INTEGER  ::   jb, igrd                         ! dummy loop indices
143      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
[4292]144      REAL(wp), POINTER :: flagu, flagv              ! short cuts
[3117]145      REAL(wp) ::   zcorr                            ! Flather correction
146      REAL(wp) ::   zforc                            ! temporary scalar
[4292]147      REAL(wp) ::   zflag, z1_2                      !    "        "
[3117]148      !!----------------------------------------------------------------------
149
[4292]150      z1_2 = 0.5_wp
151
[3117]152      ! ---------------------------------!
153      ! Flather boundary conditions     :!
154      ! ---------------------------------!
155     
156!!! REPLACE spgu with nemo_wrk work space
157
158      ! Fill temporary array with ssh data (here spgu):
159      igrd = 1
160      spgu(:,:) = 0.0
161      DO jb = 1, idx%nblenrim(igrd)
162         ii = idx%nbi(jb,igrd)
163         ij = idx%nbj(jb,igrd)
[9023]164         IF( ll_wd ) THEN
165            spgu(ii, ij) = dta%ssh(jb)  - ssh_ref 
166         ELSE
167            spgu(ii, ij) = dta%ssh(jb)
168         ENDIF
[3117]169      END DO
[4999]170
[10425]171      CALL lbc_bdy_lnk( 'bdydyn2d', spgu(:,:), 'T', 1., ib_bdy )
[3117]172      !
173      igrd = 2      ! Flather bc on u-velocity;
174      !             ! remember that flagu=-1 if normal velocity direction is outward
175      !             ! I think we should rather use after ssh ?
176      DO jb = 1, idx%nblenrim(igrd)
177         ii  = idx%nbi(jb,igrd)
178         ij  = idx%nbj(jb,igrd) 
[4292]179         flagu => idx%flagu(jb,igrd)
180         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
181         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
[3117]182         !
[4292]183         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
184
185         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
186         ! Use characteristics method instead
187         zflag = ABS(flagu)
188         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
189         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
[3117]190      END DO
191      !
192      igrd = 3      ! Flather bc on v-velocity
193      !             ! remember that flagv=-1 if normal velocity direction is outward
194      DO jb = 1, idx%nblenrim(igrd)
195         ii  = idx%nbi(jb,igrd)
196         ij  = idx%nbj(jb,igrd) 
[4292]197         flagv => idx%flagv(jb,igrd)
198         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
199         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
[3117]200         !
[4292]201         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
202
203         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
204         ! Use characteristics method instead
205         zflag = ABS(flagv)
206         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
207         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
[3117]208      END DO
[10425]209      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
210      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
[3117]211      !
212   END SUBROUTINE bdy_dyn2d_fla
[4292]213
214
[4354]215   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
[4292]216      !!----------------------------------------------------------------------
217      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
218      !!             
219      !!              - Apply Orlanski radiation condition adaptively:
220      !!                  - radiation plus weak nudging at outflow points
221      !!                  - no radiation and strong nudging at inflow points
222      !!
223      !!
224      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
225      !!----------------------------------------------------------------------
226      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
227      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
228      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
[4354]229      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
230      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
[4292]231      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
232
233      INTEGER  ::   ib, igrd                               ! dummy loop indices
234      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
235      !!----------------------------------------------------------------------
236      !
237      igrd = 2      ! Orlanski bc on u-velocity;
238      !           
239      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
240
241      igrd = 3      ! Orlanski bc on v-velocity
242      ! 
243      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
244      !
[10425]245      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
246      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
[4292]247      !
248   END SUBROUTINE bdy_dyn2d_orlanski
249
[9124]250
[4292]251   SUBROUTINE bdy_ssh( zssh )
252      !!----------------------------------------------------------------------
253      !!                  ***  SUBROUTINE bdy_ssh  ***
254      !!
255      !! ** Purpose : Duplicate sea level across open boundaries
256      !!
257      !!----------------------------------------------------------------------
258      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
259      !!
260      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
261      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
262
263      igrd = 1                       ! Everything is at T-points here
264
265      DO ib_bdy = 1, nb_bdy
266         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
267            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
268            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
269            ! Set gradient direction:
270            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
271            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
[10518]272            IF ( zcoef1+zcoef2 == 0 ) THEN   ! corner
273               zcoef = bdytmask(ii-1,ij-1) + bdytmask(ii+1,ij+1) + bdytmask(ii+1,ij-1) + bdytmask(ii-1,ij+1)
274               zssh(ii,ij) = zssh( ii-1, ij-1 ) * bdytmask( ii-1, ij-1) + &
275                 &           zssh( ii+1, ij+1 ) * bdytmask( ii+1, ij+1) + &
276                 &           zssh( ii+1, ij-1 ) * bdytmask( ii+1, ij-1) + &
277                 &           zssh( ii-1, ij+1 ) * bdytmask( ii-1, ij+1)
[4292]278               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
279            ELSE
280               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
281               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
282               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
283            ENDIF
284         END DO
285
286         ! Boundary points should be updated
[10425]287         CALL lbc_bdy_lnk( 'bdydyn2d', zssh(:,:), 'T', 1., ib_bdy )
[4292]288      END DO
289
290   END SUBROUTINE bdy_ssh
291
[3117]292   !!======================================================================
293END MODULE bdydyn2d
[4292]294
Note: See TracBrowser for help on using the repository browser.