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 branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 8981

Last change on this file since 8981 was 8981, checked in by deazer, 6 years ago

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