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.
bdyini.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3424

Last change on this file since 3424 was 3424, checked in by davestorkey, 10 years ago

Bug fix for bdyini.F90. See ticket #978.

  • Property svn:keywords set to Id
File size: 34.8 KB
Line 
1MODULE bdyini
2   !!======================================================================
3   !!                       ***  MODULE  bdyini  ***
4   !! Unstructured open boundaries : initialisation
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-01  (D. Storkey) Tidal forcing
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
13   !!----------------------------------------------------------------------
14#if defined key_bdy
15   !!----------------------------------------------------------------------
16   !!   'key_bdy'                     Unstructured Open Boundary Conditions
17   !!----------------------------------------------------------------------
18   !!   bdy_init       : Initialization of unstructured open boundaries
19   !!----------------------------------------------------------------------
20   USE timing          ! Timing
21   USE oce             ! ocean dynamics and tracers variables
22   USE dom_oce         ! ocean space and time domain
23   USE bdy_oce         ! unstructured open boundary conditions
24   USE in_out_manager  ! I/O units
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE lib_mpp         ! for mpp_sum 
27   USE iom             ! I/O
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   bdy_init   ! routine called in nemo_init
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40   
41   SUBROUTINE bdy_init
42      !!----------------------------------------------------------------------
43      !!                 ***  ROUTINE bdy_init  ***
44      !!         
45      !! ** Purpose :   Initialization of the dynamics and tracer fields with
46      !!              unstructured open boundaries.
47      !!
48      !! ** Method  :   Read initialization arrays (mask, indices) to identify
49      !!              an unstructured open boundary
50      !!
51      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
52      !!----------------------------------------------------------------------     
53      ! namelist variables
54      !-------------------
55      INTEGER, PARAMETER          :: jp_nseg = 100
56      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
57      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
58      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
59      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
60      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft
61
62      ! local variables
63      !-------------------
64      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
65      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
66      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
67      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
68      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
69      REAL   , POINTER  ::  flagu, flagv                   !    -   -
70      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
71      INTEGER, DIMENSION (2)                ::   kdimsz
72      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays
73      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
74      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
75      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
76      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
77      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
78      !!
79      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
80         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
81         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         & 
82#if defined key_lim2
83         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
84#endif
85         &             ln_vol, nn_volctl, nn_rimwidth
86      !!
87      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             &
88                             nbdysegw, jpiwob, jpjwdt, jpjwft,             &
89                             nbdysegn, jpjnob, jpindt, jpinft,             &
90                             nbdysegs, jpjsob, jpisdt, jpisft
91
92      !!----------------------------------------------------------------------
93
94      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
95
96      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
97
98      IF(lwp) WRITE(numout,*)
99      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
100      IF(lwp) WRITE(numout,*) '~~~~~~~~'
101      !
102
103      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
104         &                               ' and general open boundary condition are not compatible' )
105
106      cgrid= (/'t','u','v'/)
107
108      ! -----------------------------------------
109      ! Initialise and read namelist parameters
110      ! -----------------------------------------
111
112      nb_bdy            = 0
113      ln_coords_file(:) = .false.
114      cn_coords_file(:) = ''
115      ln_mask_file      = .false.
116      cn_mask_file(:)   = ''
117      nn_dyn2d(:)       = 0
118      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
119      nn_dyn3d(:)       = 0
120      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
121      nn_tra(:)         = 0
122      nn_tra_dta(:)     = -1  ! uninitialised flag
123#if defined key_lim2
124      nn_ice_lim2(:)    = 0
125      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
126#endif
127      ln_vol            = .false.
128      nn_volctl         = -1  ! uninitialised flag
129      nn_rimwidth(:)    = -1  ! uninitialised flag
130
131      REWIND( numnam )                   
132      READ  ( numnam, nambdy )
133
134      ! -----------------------------------------
135      ! Check and write out namelist parameters
136      ! -----------------------------------------
137
138      !                                   ! control prints
139      IF(lwp) WRITE(numout,*) '         nambdy'
140
141      IF( nb_bdy .eq. 0 ) THEN
142        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
143      ELSE
144        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
145      ENDIF
146
147      DO ib_bdy = 1,nb_bdy
148        IF(lwp) WRITE(numout,*) ' ' 
149        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
150
151        IF( ln_coords_file(ib_bdy) ) THEN
152           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
153        ELSE
154           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
155        ENDIF
156        IF(lwp) WRITE(numout,*)
157
158        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
159        SELECT CASE( nn_dyn2d(ib_bdy) )                 
160          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
161          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
162          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
163          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
164        END SELECT
165        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
166           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
167              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
168              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
169              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
170              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
171              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
172           END SELECT
173        ENDIF
174        IF(lwp) WRITE(numout,*)
175
176        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
177        SELECT CASE( nn_dyn3d(ib_bdy) )                 
178          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
179          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
181        END SELECT
182        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
183           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
184              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
185              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
186              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
187           END SELECT
188        ENDIF
189        IF(lwp) WRITE(numout,*)
190
191        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
192        SELECT CASE( nn_tra(ib_bdy) )                 
193          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
194          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
195          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
196        END SELECT
197        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
198           SELECT CASE( nn_tra_dta(ib_bdy) )                   !
199              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
200              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
201              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
202           END SELECT
203        ENDIF
204        IF(lwp) WRITE(numout,*)
205
206#if defined key_lim2
207        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
208        SELECT CASE( nn_ice_lim2(ib_bdy) )                 
209          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
210          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
211          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
212        END SELECT
213        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN
214           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !
215              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
216              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
217              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
218           END SELECT
219        ENDIF
220        IF(lwp) WRITE(numout,*)
221#endif
222
223        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
224        IF(lwp) WRITE(numout,*)
225
226      ENDDO
227
228     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
229       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
230       IF(lwp) WRITE(numout,*)
231       SELECT CASE ( nn_volctl )
232         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
233         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
234         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
235       END SELECT
236       IF(lwp) WRITE(numout,*)
237     ELSE
238       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
239       IF(lwp) WRITE(numout,*)
240     ENDIF
241
242      ! -------------------------------------------------
243      ! Initialise indices arrays for open boundaries
244      ! -------------------------------------------------
245
246      ! Work out global dimensions of boundary data
247      ! ---------------------------------------------
248      REWIND( numnam )                   
249      DO ib_bdy = 1, nb_bdy
250
251         jpbdta = 1
252         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
253 
254            ! No REWIND here because may need to read more than one nambdy_index namelist.
255            READ  ( numnam, nambdy_index )
256
257            ! Automatic boundary definition: if nbdysegX = -1
258            ! set boundary to whole side of model domain.
259            IF( nbdysege == -1 ) THEN
260               nbdysege = 1
261               jpieob(1) = jpiglo - 1
262               jpjedt(1) = 2
263               jpjeft(1) = jpjglo - 1
264            ENDIF
265            IF( nbdysegw == -1 ) THEN
266               nbdysegw = 1
267               jpiwob(1) = 2
268               jpjwdt(1) = 2
269               jpjwft(1) = jpjglo - 1
270            ENDIF
271            IF( nbdysegn == -1 ) THEN
272               nbdysegn = 1
273               jpjnob(1) = jpjglo - 1
274               jpindt(1) = 2
275               jpinft(1) = jpiglo - 1
276            ENDIF
277            IF( nbdysegs == -1 ) THEN
278               nbdysegs = 1
279               jpjsob(1) = 2
280               jpisdt(1) = 2
281               jpisft(1) = jpiglo - 1
282            ENDIF
283
284            nblendta(:,ib_bdy) = 0
285            DO iseg = 1, nbdysege
286               igrd = 1
287               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
288               igrd = 2
289               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
290               igrd = 3
291               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)               
292            ENDDO
293            DO iseg = 1, nbdysegw
294               igrd = 1
295               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
296               igrd = 2
297               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
298               igrd = 3
299               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)               
300            ENDDO
301            DO iseg = 1, nbdysegn
302               igrd = 1
303               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1               
304               igrd = 2
305               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)               
306               igrd = 3
307               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
308            ENDDO
309            DO iseg = 1, nbdysegs
310               igrd = 1
311               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
312               igrd = 2
313               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
314               igrd = 3
315               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
316            ENDDO
317
318            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
319            jpbdta = MAXVAL(nblendta(:,ib_bdy))               
320
321
322         ELSE            ! Read size of arrays in boundary coordinates file.
323
324
325            CALL iom_open( cn_coords_file(ib_bdy), inum )
326            jpbdta = 1
327            DO igrd = 1, jpbgrd
328               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
329               nblendta(igrd,ib_bdy) = kdimsz(1)
330               jpbdta = MAX(jpbdta, kdimsz(1))
331            ENDDO
332
333         ENDIF
334
335      ENDDO ! ib_bdy
336
337      ! Allocate arrays
338      !---------------
339      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
340         &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
341
342      ALLOCATE( dta_global(jpbdta, 1, jpk) )
343
344      ! Calculate global boundary index arrays or read in from file
345      !------------------------------------------------------------
346      REWIND( numnam )                   
347      DO ib_bdy = 1, nb_bdy
348
349         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
350
351            ! No REWIND here because may need to read more than one nambdy_index namelist.
352            READ  ( numnam, nambdy_index )
353
354            ! Automatic boundary definition: if nbdysegX = -1
355            ! set boundary to whole side of model domain.
356            IF( nbdysege == -1 ) THEN
357               nbdysege = 1
358               jpieob(1) = jpiglo - 1
359               jpjedt(1) = 2
360               jpjeft(1) = jpjglo - 1
361            ENDIF
362            IF( nbdysegw == -1 ) THEN
363               nbdysegw = 1
364               jpiwob(1) = 2
365               jpjwdt(1) = 2
366               jpjwft(1) = jpjglo - 1
367            ENDIF
368            IF( nbdysegn == -1 ) THEN
369               nbdysegn = 1
370               jpjnob(1) = jpjglo - 1
371               jpindt(1) = 2
372               jpinft(1) = jpiglo - 1
373            ENDIF
374            IF( nbdysegs == -1 ) THEN
375               nbdysegs = 1
376               jpjsob(1) = 2
377               jpisdt(1) = 2
378               jpisft(1) = jpiglo - 1
379            ENDIF
380
381            ! ------------ T points -------------
382            igrd = 1 
383            icount = 0
384            DO ir = 1, nn_rimwidth(ib_bdy)
385               ! east
386               DO iseg = 1, nbdysege
387                  DO ij = jpjedt(iseg), jpjeft(iseg)
388                     icount = icount + 1
389                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
390                     nbjdta(icount, igrd, ib_bdy) = ij
391                     nbrdta(icount, igrd, ib_bdy) = ir
392                  ENDDO
393               ENDDO
394               ! west
395               DO iseg = 1, nbdysegw
396                  DO ij = jpjwdt(iseg), jpjwft(iseg)
397                     icount = icount + 1
398                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
399                     nbjdta(icount, igrd, ib_bdy) = ij
400                     nbrdta(icount, igrd, ib_bdy) = ir
401                  ENDDO
402               ENDDO
403               ! north
404               DO iseg = 1, nbdysegn
405                  DO ii = jpindt(iseg), jpinft(iseg)
406                     icount = icount + 1
407                     nbidta(icount, igrd, ib_bdy) = ii
408                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
409                     nbrdta(icount, igrd, ib_bdy) = ir
410                  ENDDO
411               ENDDO
412               ! south
413               DO iseg = 1, nbdysegs
414                  DO ii = jpisdt(iseg), jpisft(iseg)
415                     icount = icount + 1
416                     nbidta(icount, igrd, ib_bdy) = ii
417                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
418                     nbrdta(icount, igrd, ib_bdy) = ir
419                  ENDDO
420               ENDDO
421            ENDDO
422
423            ! ------------ U points -------------
424            igrd = 2 
425            icount = 0
426            DO ir = 1, nn_rimwidth(ib_bdy)
427               ! east
428               DO iseg = 1, nbdysege
429                  DO ij = jpjedt(iseg), jpjeft(iseg)
430                     icount = icount + 1
431                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
432                     nbjdta(icount, igrd, ib_bdy) = ij
433                     nbrdta(icount, igrd, ib_bdy) = ir
434                  ENDDO
435               ENDDO
436               ! west
437               DO iseg = 1, nbdysegw
438                  DO ij = jpjwdt(iseg), jpjwft(iseg)
439                     icount = icount + 1
440                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
441                     nbjdta(icount, igrd, ib_bdy) = ij
442                     nbrdta(icount, igrd, ib_bdy) = ir
443                  ENDDO
444               ENDDO
445               ! north
446               DO iseg = 1, nbdysegn
447                  DO ii = jpindt(iseg), jpinft(iseg) - 1
448                     icount = icount + 1
449                     nbidta(icount, igrd, ib_bdy) = ii
450                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
451                     nbrdta(icount, igrd, ib_bdy) = ir
452                  ENDDO
453               ENDDO
454               ! south
455               DO iseg = 1, nbdysegs
456                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
457                     icount = icount + 1
458                     nbidta(icount, igrd, ib_bdy) = ii
459                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
460                     nbrdta(icount, igrd, ib_bdy) = ir
461                  ENDDO
462               ENDDO
463            ENDDO
464
465            ! ------------ V points -------------
466            igrd = 3 
467            icount = 0
468            DO ir = 1, nn_rimwidth(ib_bdy)
469               ! east
470               DO iseg = 1, nbdysege
471                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
472                     icount = icount + 1
473                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
474                     nbjdta(icount, igrd, ib_bdy) = ij
475                     nbrdta(icount, igrd, ib_bdy) = ir
476                  ENDDO
477               ENDDO
478               ! west
479               DO iseg = 1, nbdysegw
480                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
481                     icount = icount + 1
482                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
483                     nbjdta(icount, igrd, ib_bdy) = ij
484                     nbrdta(icount, igrd, ib_bdy) = ir
485                  ENDDO
486               ENDDO
487               ! north
488               DO iseg = 1, nbdysegn
489                  DO ii = jpindt(iseg), jpinft(iseg)
490                     icount = icount + 1
491                     nbidta(icount, igrd, ib_bdy) = ii
492                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
493                     nbrdta(icount, igrd, ib_bdy) = ir
494                  ENDDO
495               ENDDO
496               ! south
497               DO iseg = 1, nbdysegs
498                  DO ii = jpisdt(iseg), jpisft(iseg)
499                     icount = icount + 1
500                     nbidta(icount, igrd, ib_bdy) = ii
501                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
502                     nbrdta(icount, igrd, ib_bdy) = ir
503                  ENDDO
504               ENDDO
505            ENDDO
506
507         ELSE            ! Read global index arrays from boundary coordinates file.
508
509            DO igrd = 1, jpbgrd
510               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
511               DO ii = 1,nblendta(igrd,ib_bdy)
512                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
513               END DO
514               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
515               DO ii = 1,nblendta(igrd,ib_bdy)
516                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
517               END DO
518               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
519               DO ii = 1,nblendta(igrd,ib_bdy)
520                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
521               END DO
522
523               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
524               IF(lwp) WRITE(numout,*)
525               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
526               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
527               IF (ibr_max < nn_rimwidth(ib_bdy))   &
528                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
529
530            END DO
531            CALL iom_close( inum )
532
533         ENDIF
534
535      ENDDO 
536
537      ! Work out dimensions of boundary data on each processor
538      ! ------------------------------------------------------
539     
540      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
541      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
542      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
543      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
544
545      DO ib_bdy = 1, nb_bdy
546         DO igrd = 1, jpbgrd
547            icount  = 0
548            icountr = 0
549            idx_bdy(ib_bdy)%nblen(igrd)    = 0
550            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
551            DO ib = 1, nblendta(igrd,ib_bdy)
552               ! check that data is in correct order in file
553               ibm1 = MAX(1,ib-1)
554               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
555                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
556                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
557                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
558                  ENDIF   
559               ENDIF
560               ! check if point is in local domain
561               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
562                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
563                  !
564                  icount = icount  + 1
565                  !
566                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
567               ENDIF
568            ENDDO
569            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
570            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc       
571         ENDDO  ! igrd
572
573         ! Allocate index arrays for this boundary set
574         !--------------------------------------------
575         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
576         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
577         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
578         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
579         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
580         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
581         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
582         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
583
584         ! Dispatch mapping indices and discrete distances on each processor
585         ! -----------------------------------------------------------------
586
587         DO igrd = 1, jpbgrd
588            icount  = 0
589            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
590            DO ir=1, nn_rimwidth(ib_bdy)
591               DO ib = 1, nblendta(igrd,ib_bdy)
592                  ! check if point is in local domain and equals ir
593                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
594                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
595                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
596                     !
597                     icount = icount  + 1
598                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
599                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
600                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
601                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
602                  ENDIF
603               ENDDO
604            ENDDO
605         ENDDO 
606
607         ! Compute rim weights for FRS scheme
608         ! ----------------------------------
609         DO igrd = 1, jpbgrd
610            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
611               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
612               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
613!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
614!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
615            END DO
616         END DO
617
618      ENDDO
619
620      ! ------------------------------------------------------
621      ! Initialise masks and find normal/tangential directions
622      ! ------------------------------------------------------
623
624      ! Read global 2D mask at T-points: bdytmask
625      ! -----------------------------------------
626      ! bdytmask = 1  on the computational domain AND on open boundaries
627      !          = 0  elsewhere   
628 
629      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
630         zmask(         :                ,:) = 0.e0
631         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0         
632      ELSE IF( ln_mask_file ) THEN
633         CALL iom_open( cn_mask_file, inum )
634         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
635         CALL iom_close( inum )
636      ELSE
637         zmask(:,:) = 1.e0
638      ENDIF
639
640      DO ij = 1, nlcj      ! Save mask over local domain     
641         DO ii = 1, nlci
642            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
643         END DO
644      END DO
645
646      ! Derive mask on U and V grid from mask on T grid
647      bdyumask(:,:) = 0.e0
648      bdyvmask(:,:) = 0.e0
649      DO ij=1, jpjm1
650         DO ii=1, jpim1
651            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
652            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
653         END DO
654      END DO
655      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
656
657
658      ! Mask corrections
659      ! ----------------
660      DO ik = 1, jpkm1
661         DO ij = 1, jpj
662            DO ii = 1, jpi
663               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
664               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
665               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
666               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
667            END DO     
668         END DO
669      END DO
670
671      DO ik = 1, jpkm1
672         DO ij = 2, jpjm1
673            DO ii = 2, jpim1
674               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
675                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
676            END DO     
677         END DO
678      END DO
679
680      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
681      bdytmask(:,:) = tmask(:,:,1)
682
683      ! bdy masks and bmask are now set to zero on boundary points:
684      igrd = 1       ! In the free surface case, bmask is at T-points
685      DO ib_bdy = 1, nb_bdy
686        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
687          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
688        ENDDO
689      ENDDO
690      !
691      igrd = 1
692      DO ib_bdy = 1, nb_bdy
693        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
694          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
695        ENDDO
696      ENDDO
697      !
698      igrd = 2
699      DO ib_bdy = 1, nb_bdy
700        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
701          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
702        ENDDO
703      ENDDO
704      !
705      igrd = 3
706      DO ib_bdy = 1, nb_bdy
707        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
708          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
709        ENDDO
710      ENDDO
711
712      ! Lateral boundary conditions
713      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
714      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
715
716      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
717
718         idx_bdy(ib_bdy)%flagu(:) = 0.e0
719         idx_bdy(ib_bdy)%flagv(:) = 0.e0
720         icount = 0 
721
722         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
723         !flagu =  0 : u is tangential
724         !flagu =  1 : u is normal to the boundary and is direction is inward
725 
726         igrd = 2      ! u-component
727         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
728            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
729            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
730            zefl = bdytmask(nbi  ,nbj)
731            zwfl = bdytmask(nbi+1,nbj)
732            IF( zefl + zwfl == 2 ) THEN
733               icount = icount + 1
734            ELSE
735               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
736            ENDIF
737         END DO
738
739         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
740         !flagv =  0 : u is tangential
741         !flagv =  1 : u is normal to the boundary and is direction is inward
742
743         igrd = 3      ! v-component
744         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
745            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
746            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
747            znfl = bdytmask(nbi,nbj  )
748            zsfl = bdytmask(nbi,nbj+1)
749            IF( znfl + zsfl == 2 ) THEN
750               icount = icount + 1
751            ELSE
752               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
753            END IF
754         END DO
755 
756         IF( icount /= 0 ) THEN
757            IF(lwp) WRITE(numout,*)
758            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
759               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
760            IF(lwp) WRITE(numout,*) ' ========== '
761            IF(lwp) WRITE(numout,*)
762            nstop = nstop + 1
763         ENDIF
764   
765      ENDDO
766
767      ! Compute total lateral surface for volume correction:
768      ! ----------------------------------------------------
769      bdysurftot = 0.e0 
770      IF( ln_vol ) THEN 
771         igrd = 2      ! Lateral surface at U-points
772         DO ib_bdy = 1, nb_bdy
773            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
774               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
775               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
776               flagu => idx_bdy(ib_bdy)%flagu(ib)
777               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
778                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
779                  &                    * tmask_i(nbi  , nbj)                           &
780                  &                    * tmask_i(nbi+1, nbj)                   
781            ENDDO
782         ENDDO
783
784         igrd=3 ! Add lateral surface at V-points
785         DO ib_bdy = 1, nb_bdy
786            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
787               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
788               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
789               flagv => idx_bdy(ib_bdy)%flagv(ib)
790               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
791                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
792                  &                    * tmask_i(nbi, nbj  )                           &
793                  &                    * tmask_i(nbi, nbj+1)
794            ENDDO
795         ENDDO
796         !
797         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
798      END IF   
799      !
800      ! Tidy up
801      !--------
802      DEALLOCATE(nbidta, nbjdta, nbrdta)
803
804      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
805
806   END SUBROUTINE bdy_init
807
808#else
809   !!---------------------------------------------------------------------------------
810   !!   Dummy module                                   NO open boundaries
811   !!---------------------------------------------------------------------------------
812CONTAINS
813   SUBROUTINE bdy_init      ! Dummy routine
814   END SUBROUTINE bdy_init
815#endif
816
817   !!=================================================================================
818END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.