source: branches/2013/dev_r3891_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 3902

Last change on this file since 3902 was 3902, checked in by davestorkey, 9 years ago

bug fixes

File size: 10.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   !!----------------------------------------------------------------------
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 )
60         CASE DEFAULT
61            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
62         END SELECT
63      ENDDO
64
65   END SUBROUTINE bdy_dyn2d
66
67   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy )
68      !!----------------------------------------------------------------------
69      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
70      !!
71      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
72      !!                at open boundaries.
73      !!
74      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
75      !!               a three-dimensional baroclinic ocean model with realistic
76      !!               topography. Tellus, 365-382.
77      !!----------------------------------------------------------------------
78      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
79      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
80      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
81      !!
82      INTEGER  ::   jb, jk         ! dummy loop indices
83      INTEGER  ::   ii, ij, igrd   ! local integers
84      REAL(wp) ::   zwgt           ! boundary weight
85      !!----------------------------------------------------------------------
86      !
87      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_frs')
88      !
89      igrd = 2                      ! Relaxation of zonal velocity
90      DO jb = 1, idx%nblen(igrd)
91         ii   = idx%nbi(jb,igrd)
92         ij   = idx%nbj(jb,igrd)
93         zwgt = idx%nbw(jb,igrd)
94         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
95      END DO
96      !
97      igrd = 3                      ! Relaxation of meridional velocity
98      DO jb = 1, idx%nblen(igrd)
99         ii   = idx%nbi(jb,igrd)
100         ij   = idx%nbj(jb,igrd)
101         zwgt = idx%nbw(jb,igrd)
102         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
103      END DO
104      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy ) 
105      CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
106      !
107      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
108      !
109
110   END SUBROUTINE bdy_dyn2d_frs
111
112
113   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy )
114      !!----------------------------------------------------------------------
115      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
116      !!             
117      !!              - Apply Flather boundary conditions on normal barotropic velocities
118      !!
119      !! ** WARNINGS about FLATHER implementation:
120      !!1. According to Palma and Matano, 1998 "after ssh" is used.
121      !!   In ROMS and POM implementations, it is "now ssh". In the current
122      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
123      !!   So I use "before ssh" in the following.
124      !!
125      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
126      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
127      !!   ssh in the code is not updated).
128      !!
129      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
130      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
131      !!----------------------------------------------------------------------
132      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
133      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
134      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
135
136      INTEGER  ::   jb, igrd                         ! dummy loop indices
137      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
138      REAL(wp) ::   zcorr                            ! Flather correction
139      REAL(wp) ::   zforc                            ! temporary scalar
140      !!----------------------------------------------------------------------
141
142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
143
144      ! ---------------------------------!
145      ! Flather boundary conditions     :!
146      ! ---------------------------------!
147     
148!!! REPLACE spgu with nemo_wrk work space
149
150      ! Fill temporary array with ssh data (here spgu):
151      igrd = 1
152      spgu(:,:) = 0.0
153      DO jb = 1, idx%nblenrim(igrd)
154         ii = idx%nbi(jb,igrd)
155         ij = idx%nbj(jb,igrd)
156         spgu(ii, ij) = dta%ssh(jb)
157      END DO
158      !
159      igrd = 2      ! Flather bc on u-velocity;
160      !             ! remember that flagu=-1 if normal velocity direction is outward
161      !             ! I think we should rather use after ssh ?
162      DO jb = 1, idx%nblenrim(igrd)
163         ii  = idx%nbi(jb,igrd)
164         ij  = idx%nbj(jb,igrd) 
165         flagu => idx%flagu(jb,igrd)
166         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
167         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
168         !
169         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
170         zforc = dta%u2d(jb)
171         pua2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
172      END DO
173      !
174      igrd = 3      ! Flather bc on v-velocity
175      !             ! remember that flagv=-1 if normal velocity direction is outward
176      DO jb = 1, idx%nblenrim(igrd)
177         ii  = idx%nbi(jb,igrd)
178         ij  = idx%nbj(jb,igrd) 
179         flagv => idx%flagv(jb,igrd)
180         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
181         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
182         !
183         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
184         zforc = dta%v2d(jb)
185         pva2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
186      END DO
187      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
188      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
189      !
190      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
191      !
192   END SUBROUTINE bdy_dyn2d_fla
193
194
195   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy )
196      !!----------------------------------------------------------------------
197      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
198      !!             
199      !!              - Apply Orlanski radiation condition adaptively:
200      !!                  - radiation plus weak nudging at outflow points
201      !!                  - no radiation and strong nudging at inflow points
202      !!
203      !!
204      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
205      !!----------------------------------------------------------------------
206      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
207      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
208      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
209
210      INTEGER  ::   jb, igrd                               ! dummy loop indices
211      !!----------------------------------------------------------------------
212
213      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
214      !
215      igrd = 2      ! Orlanski bc on u-velocity;
216      !           
217      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d )
218
219      igrd = 3      ! Orlanski bc on v-velocity
220     
221      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d )
222      !
223      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
224      !
225
226   END SUBROUTINE bdy_dyn2d_orlanski
227
228
229#else
230   !!----------------------------------------------------------------------
231   !!   Dummy module                   NO Unstruct Open Boundary Conditions
232   !!----------------------------------------------------------------------
233CONTAINS
234   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
235      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
236   END SUBROUTINE bdy_dyn2d
237#endif
238
239   !!======================================================================
240END MODULE bdydyn2d
Note: See TracBrowser for help on using the repository browser.