source: branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 4007

Last change on this file since 4007 was 4007, checked in by davestorkey, 8 years ago
  1. Bug fixes for flagu/flagv calculation in bdyini.F90.
  2. Introduce masking of derivatives in radiation velocity calculation in bdylib.F90
  3. Change relaxation term in Orlanski implementation to explicit timestepping in bdylib.F90.
  4. Remove bdyfmask (redundant).
File size: 11.2 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   !!----------------------------------------------------------------------
9#if defined key_bdy 
10   !!----------------------------------------------------------------------
11   !!   'key_bdy' :                    Unstructured Open Boundary Condition
12   !!----------------------------------------------------------------------
13   !!   bdy_dyn2d      : Apply open boundary conditions to barotropic variables.
14   !!   bdy_dyn2d_fla    : Apply Flather condition
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 dynspg_oce      ! for barotropic variables
22   USE phycst          ! physical constants
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
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
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE bdy_dyn2d( kt )
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      !!
47      INTEGER                                  ::   ib_bdy ! Loop counter
48
49      DO ib_bdy=1, nb_bdy
50
51         SELECT CASE( cn_dyn2d(ib_bdy) )
52         CASE('none')
53            CYCLE
54         CASE('frs')
55            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
56         CASE('flather')
57            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
58         CASE('orlanski')
59            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
60         CASE('orlanski_npo')
61            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
62         CASE DEFAULT
63            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
64         END SELECT
65      ENDDO
66
67   END SUBROUTINE bdy_dyn2d
68
69   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy )
70      !!----------------------------------------------------------------------
71      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
72      !!
73      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
74      !!                at open boundaries.
75      !!
76      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
77      !!               a three-dimensional baroclinic ocean model with realistic
78      !!               topography. Tellus, 365-382.
79      !!----------------------------------------------------------------------
80      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
81      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
82      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
83      !!
84      INTEGER  ::   jb, jk         ! dummy loop indices
85      INTEGER  ::   ii, ij, igrd   ! local integers
86      REAL(wp) ::   zwgt           ! boundary weight
87      !!----------------------------------------------------------------------
88      !
89      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_frs')
90      !
91      igrd = 2                      ! Relaxation of zonal velocity
92      DO jb = 1, idx%nblen(igrd)
93         ii   = idx%nbi(jb,igrd)
94         ij   = idx%nbj(jb,igrd)
95         zwgt = idx%nbw(jb,igrd)
96         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
97      END DO
98      !
99      igrd = 3                      ! Relaxation of meridional 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         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
105      END DO
106      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) 
107      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
108      !
109      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
110      !
111
112   END SUBROUTINE bdy_dyn2d_frs
113
114
115   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy )
116      !!----------------------------------------------------------------------
117      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
118      !!             
119      !!              - Apply Flather boundary conditions on normal barotropic velocities
120      !!
121      !! ** WARNINGS about FLATHER implementation:
122      !!1. According to Palma and Matano, 1998 "after ssh" is used.
123      !!   In ROMS and POM implementations, it is "now ssh". In the current
124      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
125      !!   So I use "before ssh" in the following.
126      !!
127      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
128      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
129      !!   ssh in the code is not updated).
130      !!
131      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
132      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
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
138      INTEGER  ::   jb, igrd                         ! dummy loop indices
139      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
140      REAL(wp), POINTER :: flagu, flagv              ! short cuts
141      REAL(wp) ::   zcorr                            ! Flather correction
142      REAL(wp) ::   zforc                            ! temporary scalar
143      !!----------------------------------------------------------------------
144
145      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
146
147      ! ---------------------------------!
148      ! Flather boundary conditions     :!
149      ! ---------------------------------!
150     
151!!! REPLACE spgu with nemo_wrk work space
152
153      ! Fill temporary array with ssh data (here spgu):
154      igrd = 1
155      spgu(:,:) = 0.0
156      DO jb = 1, idx%nblenrim(igrd)
157         ii = idx%nbi(jb,igrd)
158         ij = idx%nbj(jb,igrd)
159         spgu(ii, ij) = dta%ssh(jb)
160      END DO
161      !
162      igrd = 2      ! Flather bc on u-velocity;
163      !             ! remember that flagu=-1 if normal velocity direction is outward
164      !             ! I think we should rather use after ssh ?
165      DO jb = 1, idx%nblenrim(igrd)
166         ii  = idx%nbi(jb,igrd)
167         ij  = idx%nbj(jb,igrd) 
168         flagu => idx%flagu(jb,igrd)
169         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
170         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
171         !
172         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
173         zforc = dta%u2d(jb)
174         pua2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
175      END DO
176      !
177      igrd = 3      ! Flather bc on v-velocity
178      !             ! remember that flagv=-1 if normal velocity direction is outward
179      DO jb = 1, idx%nblenrim(igrd)
180         ii  = idx%nbi(jb,igrd)
181         ij  = idx%nbj(jb,igrd) 
182         flagv => idx%flagv(jb,igrd)
183         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
184         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
185         !
186         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
187         zforc = dta%v2d(jb)
188         pva2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
189      END DO
190      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
191      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
192      !
193      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
194      !
195   END SUBROUTINE bdy_dyn2d_fla
196
197
198   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo )
199      !!----------------------------------------------------------------------
200      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
201      !!             
202      !!              - Apply Orlanski radiation condition adaptively:
203      !!                  - radiation plus weak nudging at outflow points
204      !!                  - no radiation and strong nudging at inflow points
205      !!
206      !!
207      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
208      !!----------------------------------------------------------------------
209      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
210      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
211      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
212      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
213
214      INTEGER  ::   ib, igrd                               ! dummy loop indices
215      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
216      !!----------------------------------------------------------------------
217
218      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
219      !
220      igrd = 2      ! Orlanski bc on u-velocity;
221      !           
222      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
223
224      igrd = 3      ! Orlanski bc on v-velocity
225     
226      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
227      !
228      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
229      !
230      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
231      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
232      !
233      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
234      !
235   END SUBROUTINE bdy_dyn2d_orlanski
236
237#else
238   !!----------------------------------------------------------------------
239   !!   Dummy module                   NO Unstruct Open Boundary Conditions
240   !!----------------------------------------------------------------------
241CONTAINS
242   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
243      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
244   END SUBROUTINE bdy_dyn2d
245#endif
246
247   !!======================================================================
248END MODULE bdydyn2d
Note: See TracBrowser for help on using the repository browser.