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/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 3970

Last change on this file since 3970 was 3970, checked in by cbricaud, 11 years ago

Time splitting update, see ticket #1079
File size: 12.0 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#if defined key_bdy 
11   !!----------------------------------------------------------------------
12   !!   'key_bdy' :                    Unstructured Open Boundary Condition
13   !!----------------------------------------------------------------------
14   !!   bdy_dyn2d      : Apply open boundary conditions to barotropic variables.
15   !!   bdy_dyn2d_fla    : Apply Flather condition
16   !!----------------------------------------------------------------------
17   USE timing          ! Timing
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE bdy_oce         ! ocean open boundary conditions
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   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
34   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
35   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE bdy_dyn2d( kt )
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      !!
48      INTEGER                                  ::   ib_bdy ! Loop counter
49
50      DO ib_bdy=1, nb_bdy
51
52         SELECT CASE( nn_dyn2d(ib_bdy) )
53         CASE(jp_none)
54            CYCLE
55         CASE(jp_frs)
56            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
57         CASE(jp_flather)
58            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
59         CASE DEFAULT
60            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
61         END SELECT
62      ENDDO
63
64   END SUBROUTINE bdy_dyn2d
65
66   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy )
67      !!----------------------------------------------------------------------
68      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
69      !!
70      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
71      !!                at open boundaries.
72      !!
73      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
74      !!               a three-dimensional baroclinic ocean model with realistic
75      !!               topography. Tellus, 365-382.
76      !!----------------------------------------------------------------------
77      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
78      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
79      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
80      !!
81      INTEGER  ::   jb, jk         ! dummy loop indices
82      INTEGER  ::   ii, ij, igrd   ! local integers
83      REAL(wp) ::   zwgt           ! boundary weight
84      !!----------------------------------------------------------------------
85      !
86      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_frs')
87      !
88      igrd = 2                      ! Relaxation of zonal velocity
89      DO jb = 1, idx%nblen(igrd)
90         ii   = idx%nbi(jb,igrd)
91         ij   = idx%nbj(jb,igrd)
92         zwgt = idx%nbw(jb,igrd)
93         pu2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1)
94      END DO
95      !
96      igrd = 3                      ! Relaxation of meridional 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)
101         pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1)
102      END DO
103      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy ) 
104      CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
105      !
106      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
107      !
108
109   END SUBROUTINE bdy_dyn2d_frs
110
111
112   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy )
113      !!----------------------------------------------------------------------
114      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
115      !!             
116      !!              - Apply Flather boundary conditions on normal barotropic velocities
117      !!
118      !! ** WARNINGS about FLATHER implementation:
119      !!1. According to Palma and Matano, 1998 "after ssh" is used.
120      !!   In ROMS and POM implementations, it is "now ssh". In the current
121      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.
122      !!   So I use "before ssh" in the following.
123      !!
124      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
125      !!   fact, the model ssh just inside the dynamical boundary is used (the outside 
126      !!   ssh in the code is not updated).
127      !!
128      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
129      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
130      !!----------------------------------------------------------------------
131      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
132      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
133      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
134
135      INTEGER  ::   jb, igrd                         ! dummy loop indices
136      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
137      REAL(wp) ::   zcorr                            ! Flather correction
138      REAL(wp) ::   zforc                            ! temporary scalar
139      REAL(wp) ::   zflag, z1_2                      !    "        "
140      !!----------------------------------------------------------------------
141
142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
143
144      z1_2 = 0.5_wp
145
146      ! ---------------------------------!
147      ! Flather boundary conditions     :!
148      ! ---------------------------------!
149     
150!!! REPLACE spgu with nemo_wrk work space
151
152      ! Fill temporary array with ssh data (here spgu):
153      igrd = 1
154      spgu(:,:) = 0.0
155      DO jb = 1, idx%nblenrim(igrd)
156         ii = idx%nbi(jb,igrd)
157         ij = idx%nbj(jb,igrd)
158         spgu(ii, ij) = dta%ssh(jb)
159      END DO
160      !
161      igrd = 2      ! Flather bc on u-velocity;
162      !             ! remember that flagu=-1 if normal velocity direction is outward
163      !             ! I think we should rather use after ssh ?
164      DO jb = 1, idx%nblenrim(igrd)
165         ii  = idx%nbi(jb,igrd)
166         ij  = idx%nbj(jb,igrd) 
167         iim1 = ii + MAX( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice inside the boundary
168         iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice outside the boundary
169         !
170         zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
171         ! bg jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
172!!         zforc = dta%u2d(jb)
173         zflag = ABS(idx%flagu(jb))
174         iim1 = ii + idx%flagu(jb)
175         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pu2d(iim1,ij)
176         pu2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
177         ! end jchanut tschanges
178      END DO
179      !
180      igrd = 3      ! Flather bc on v-velocity
181      !             ! remember that flagv=-1 if normal velocity direction is outward
182      DO jb = 1, idx%nblenrim(igrd)
183         ii  = idx%nbi(jb,igrd)
184         ij  = idx%nbj(jb,igrd) 
185         ijm1 = ij + MAX( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice inside the boundary
186         ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice outside the boundary
187         !
188         zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
189         ! bg jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
190!!         zforc = dta%v2d(jb)
191         zflag = ABS(idx%flagv(jb))
192         ijm1 = ij + idx%flagv(jb)
193         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pv2d(ii,ijm1)
194         pv2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
195         ! end jchanut tschanges
196      END DO
197      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
198      CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy )   !
199      !
200      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
201      !
202   END SUBROUTINE bdy_dyn2d_fla
203
204   SUBROUTINE bdy_ssh( zssh )
205      !!----------------------------------------------------------------------
206      !!                  ***  SUBROUTINE bdy_ssh  ***
207      !!
208      !! ** Purpose : Duplicate sea level across open boundaries
209      !!
210      !!----------------------------------------------------------------------
211      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
212      !!
213      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
214      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
215
216      igrd = 1                       ! Everything is at T-points here
217
218      DO ib_bdy = 1, nb_bdy
219         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
220            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
221            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
222            ! Set gradient direction:
223            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
224            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
225            IF ( zcoef1+zcoef2 == 0 ) THEN
226               ! corner
227!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
228!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
229!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
230!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
231!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
232               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
233               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
234                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
235                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
236                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
237               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
238            ELSE
239               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
240               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
241               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
242            ENDIF
243         END DO
244
245         ! Boundary points should be updated
246         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
247      END DO
248
249   END SUBROUTINE bdy_ssh
250#else
251   !!----------------------------------------------------------------------
252   !!   Dummy module                   NO Unstruct Open Boundary Conditions
253   !!----------------------------------------------------------------------
254CONTAINS
255   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
256      INTEGER, intent(in) :: kt
257      WRITE(*,*) 'bdy_dyn2: You should not have seen this print! error?', kt
258   END SUBROUTINE bdy_dyn2d
259
260#endif
261
262   !!======================================================================
263END MODULE bdydyn2d
264
Note: See TracBrowser for help on using the repository browser.