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

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90 @ 7753

Last change on this file since 7753 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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