source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90 @ 11048

Last change on this file since 11048 was 11048, checked in by girrmann, 19 months ago

dev_r10984_HPC-13 : Step 1, boundary is now detected all over the local domain, this does not change the result. Improve bdy treatment for bdy_rnf in bdytra.F90, this changes the result when keyword runoff is specified in namelist

  • Property svn:keywords set to Id
File size: 13.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   !!            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 oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE bdy_oce         ! ocean open boundary conditions
19   USE bdylib          ! BDY library routines
20   USE phycst          ! physical constants
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   USE lib_mpp, ONLY: ctl_stop
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/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
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      igrd = 2                      ! Relaxation of zonal velocity
98      DO jb = 1, idx%nblen(igrd)
99         ii   = idx%nbi(jb,igrd)
100         ij   = idx%nbj(jb,igrd)
101         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
102         zwgt = idx%nbw(jb,igrd)
103         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
104      END DO
105      !
106      igrd = 3                      ! Relaxation of meridional velocity
107      DO jb = 1, idx%nblen(igrd)
108         ii   = idx%nbi(jb,igrd)
109         ij   = idx%nbj(jb,igrd)
110         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
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( 'bdydyn2d', pua2d, 'U', -1., ib_bdy ) 
115      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
116      !
117   END SUBROUTINE bdy_dyn2d_frs
118
119
120   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )
121      !!----------------------------------------------------------------------
122      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
123      !!             
124      !!              - Apply Flather boundary conditions on normal barotropic velocities
125      !!
126      !! ** WARNINGS about FLATHER implementation:
127      !!1. According to Palma and Matano, 1998 "after ssh" is used.
128      !!   In ROMS and POM implementations, it is "now ssh". In the current
129      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
130      !!   So I use "before ssh" in the following.
131      !!
132      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
133      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
134      !!   ssh in the code is not updated).
135      !!
136      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
137      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
138      !!----------------------------------------------------------------------
139      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
140      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
141      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
142      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
143      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh, phur, phvr 
144
145      INTEGER  ::   jb, igrd                         ! dummy loop indices
146      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
147      REAL(wp), POINTER :: flagu, flagv              ! short cuts
148      REAL(wp) ::   zcorr                            ! Flather correction
149      REAL(wp) ::   zforc                            ! temporary scalar
150      REAL(wp) ::   zflag, z1_2                      !    "        "
151      !!----------------------------------------------------------------------
152
153      z1_2 = 0.5_wp
154
155      ! ---------------------------------!
156      ! Flather boundary conditions     :!
157      ! ---------------------------------!
158     
159!!! REPLACE spgu with nemo_wrk work space
160
161      ! Fill temporary array with ssh data (here spgu):
162      igrd = 1
163      spgu(:,:) = 0.0
164      DO jb = 1, idx%nblenrim(igrd)
165         ii = idx%nbi(jb,igrd)
166         ij = idx%nbj(jb,igrd)
167         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
168         IF( ll_wd ) THEN
169            spgu(ii, ij) = dta%ssh(jb)  - ssh_ref 
170         ELSE
171            spgu(ii, ij) = dta%ssh(jb)
172         ENDIF
173      END DO
174
175      CALL lbc_bdy_lnk( 'bdydyn2d', spgu(:,:), 'T', 1., ib_bdy )
176      !
177      igrd = 2      ! Flather bc on u-velocity;
178      !             ! remember that flagu=-1 if normal velocity direction is outward
179      !             ! I think we should rather use after ssh ?
180      DO jb = 1, idx%nblenrim(igrd)
181         ii  = idx%nbi(jb,igrd)
182         ij  = idx%nbj(jb,igrd) 
183         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
184         flagu => idx%flagu(jb,igrd)
185         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
186         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary
187         !
188         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
189
190         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
191         ! Use characteristics method instead
192         zflag = ABS(flagu)
193         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
194         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
195      END DO
196      !
197      igrd = 3      ! Flather bc on v-velocity
198      !             ! remember that flagv=-1 if normal velocity direction is outward
199      DO jb = 1, idx%nblenrim(igrd)
200         ii  = idx%nbi(jb,igrd)
201         ij  = idx%nbj(jb,igrd) 
202         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
203         flagv => idx%flagv(jb,igrd)
204         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
205         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary
206         !
207         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
208
209         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
210         ! Use characteristics method instead
211         zflag = ABS(flagv)
212         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
213         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
214      END DO
215      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
216      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
217      !
218   END SUBROUTINE bdy_dyn2d_fla
219
220
221   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
222      !!----------------------------------------------------------------------
223      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
224      !!             
225      !!              - Apply Orlanski radiation condition adaptively:
226      !!                  - radiation plus weak nudging at outflow points
227      !!                  - no radiation and strong nudging at inflow points
228      !!
229      !!
230      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
231      !!----------------------------------------------------------------------
232      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
233      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
234      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
235      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
236      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
237      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
238
239      INTEGER  ::   ib, igrd                               ! dummy loop indices
240      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
241      !!----------------------------------------------------------------------
242      !
243      igrd = 2      ! Orlanski bc on u-velocity;
244      !           
245      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
246
247      igrd = 3      ! Orlanski bc on v-velocity
248     
249      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
250      !
251      CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
252      CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy )   !
253      !
254   END SUBROUTINE bdy_dyn2d_orlanski
255
256
257   SUBROUTINE bdy_ssh( zssh )
258      !!----------------------------------------------------------------------
259      !!                  ***  SUBROUTINE bdy_ssh  ***
260      !!
261      !! ** Purpose : Duplicate sea level across open boundaries
262      !!
263      !!----------------------------------------------------------------------
264      REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) ::   zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
265      !!
266      INTEGER  ::   ib_bdy          ! bdy index
267      !!----------------------------------------------------------------------
268      DO ib_bdy = 1, nb_bdy
269         CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh )   ! zssh is masked
270         CALL lbc_bdy_lnk( 'bdydyn2d', zssh(:,:,1), 'T', 1., ib_bdy )
271      END DO
272      !
273   END SUBROUTINE bdy_ssh
274
275   !!======================================================================
276END MODULE bdydyn2d
277
Note: See TracBrowser for help on using the repository browser.