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.
obcini.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90 @ 2814

Last change on this file since 2814 was 2814, checked in by davestorkey, 13 years ago
  1. Implement tidal harmonics forcing (UKMO version) in new structure.
  2. Other bug fixes and updates.
  • Property svn:keywords set to Id
File size: 23.8 KB
Line 
1MODULE obcini
2   !!======================================================================
3   !!                       ***  MODULE  obcini  ***
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, J. Chanut) OBC-BDY merge
13   !!                 !  --- Renamed bdyini.F90 -> obcini.F90 ---
14   !!----------------------------------------------------------------------
15#if defined key_obc
16   !!----------------------------------------------------------------------
17   !!   'key_obc'                     Unstructured Open Boundary Conditions
18   !!----------------------------------------------------------------------
19   !!   obc_init       : Initialization of unstructured open boundaries
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers variables
22   USE dom_oce         ! ocean space and time domain
23   USE obc_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   obc_init   ! routine called by opa.F90
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 obc_init
42      !!----------------------------------------------------------------------
43      !!                 ***  ROUTINE obc_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   :  obc_init.nc, input file for unstructured open boundaries
52      !!----------------------------------------------------------------------     
53      INTEGER  ::   ib_obc, ii, ij, ik, igrd, ib, ir   ! dummy loop indices
54      INTEGER  ::   icount, icountr, ibr_max, ilen1    ! local integers
55      INTEGER  ::   iw, ie, is, in, inum, id_dummy     !   -       -
56      INTEGER  ::   igrd_start, igrd_end, jpbdta       !   -       -
57      INTEGER, POINTER  ::  nbi, nbj, nbr              ! short cuts
58      REAL   , POINTER  ::  flagu, flagv               !    -   -
59      REAL(wp) ::   zefl, zwfl, znfl, zsfl             ! local scalars
60      INTEGER, DIMENSION (2)                ::   kdimsz
61      INTEGER, DIMENSION(jpbgrd,jp_obc)       ::   nblendta         ! Length of index arrays
62      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of obc dta
63      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
64      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
65      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
66      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
67      !!
68      NAMELIST/namobc/ nb_obc, ln_coords_file, cn_coords_file,             &
69         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn3d,     &
70         &             nn_tra,                                             &
71#if defined key_lim2
72         &             nn_ice_lim2,                                        &
73#endif
74         &             nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl,    &
75         &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             &
76         &             nn_dmp3d_in, nn_dmp3d_out
77      !!----------------------------------------------------------------------
78
79      IF( obc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'obc_init : unable to allocate oce arrays' )
80
81      IF(lwp) WRITE(numout,*)
82      IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries'
83      IF(lwp) WRITE(numout,*) '~~~~~~~~'
84      !
85
86      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
87         &                               ' and general open boundary condition are not compatible' )
88
89      cgrid= (/'T','U','V'/)
90
91      ! -----------------------------------------
92      ! Initialise and read namelist parameters
93      ! -----------------------------------------
94
95      nb_obc            = 0
96      ln_coords_file(:) = .false.
97      cn_coords_file(:) = ''
98      ln_mask_file      = .false.
99      cn_mask_file(:)   = ''
100      nn_dyn2d(:)       = 0
101      nn_dyn3d(:)       = 0
102      nn_tra(:)         = 0
103#if defined key_lim2
104      nn_ice_lim2(:)    = 0
105#endif
106      ln_vol            = .false.
107      ln_clim(:)        = .false.
108      nn_dtactl(:)      = -1  ! uninitialised flag
109      nn_tides(:)       =  0  ! default to no tidal forcing
110      nn_volctl         = -1  ! uninitialised flag
111      nn_rimwidth(:)    = -1  ! uninitialised flag
112      nn_dmp2d_in(:)    = -1  ! uninitialised flag
113      nn_dmp2d_out(:)   = -1  ! uninitialised flag
114      nn_dmp3d_in(:)    = -1  ! uninitialised flag
115      nn_dmp3d_out(:)   = -1  ! uninitialised flag
116
117      REWIND( numnam )                   
118      READ  ( numnam, namobc )
119
120      ! -----------------------------------------
121      ! Check and write out namelist parameters
122      ! -----------------------------------------
123
124      !                                   ! control prints
125      IF(lwp) WRITE(numout,*) '         namobc'
126
127      IF( nb_obc .eq. 0 ) THEN
128        IF(lwp) WRITE(numout,*) 'nb_obc = 0, NO OPEN BOUNDARIES APPLIED.'
129      ELSE
130        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_obc
131      ENDIF
132
133      DO ib_obc = 1,nb_obc
134        IF(lwp) WRITE(numout,*) ' ' 
135        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------' 
136
137        !                                         ! check type of data used (nn_dtactl value)
138        SELECT CASE( nn_dtactl(ib_obc) )                   !
139          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for obc data'       
140          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
141          CASE DEFAULT   ;   CALL ctl_stop( 'nn_dtactl must be 0 or 1' )
142        END SELECT
143        IF(lwp) WRITE(numout,*)
144
145        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
146        SELECT CASE( nn_dyn2d(ib_obc) )                 
147          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
148          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
149          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
150          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
151        END SELECT
152        IF(lwp) WRITE(numout,*)
153
154        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
155        SELECT CASE( nn_dyn3d(ib_obc) )                 
156          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
157          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
158          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
159        END SELECT
160        IF(lwp) WRITE(numout,*)
161
162        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
163        SELECT CASE( nn_tra(ib_obc) )                 
164          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
165          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
166          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
167        END SELECT
168        IF(lwp) WRITE(numout,*)
169
170#if defined key_lim2
171        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
172        SELECT CASE( nn_tra(ib_obc) )                 
173          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
174          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
175          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
176        END SELECT
177        IF(lwp) WRITE(numout,*)
178#endif
179
180        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth
181        IF(lwp) WRITE(numout,*)
182
183        SELECT CASE( nn_tides(ib_obc) ) 
184        CASE(0)
185          IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing'
186          IF(lwp) WRITE(numout,*)
187        CASE(1)
188          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution'
189          IF(lwp) WRITE(numout,*)
190        CASE(2)
191          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions'
192          IF(lwp) WRITE(numout,*)
193        CASE DEFAULT
194          CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' )
195        END SELECT
196
197      ENDDO
198
199     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
200       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
201       IF(lwp) WRITE(numout,*)
202       SELECT CASE ( nn_volctl )
203         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
204         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
205         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
206       END SELECT
207       IF(lwp) WRITE(numout,*)
208     ELSE
209       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
210       IF(lwp) WRITE(numout,*)
211     ENDIF
212
213      ! -------------------------------------------------
214      ! Initialise indices arrays for open boundaries
215      ! -------------------------------------------------
216
217      ! Work out global dimensions of boundary data
218      ! ---------------------------------------------
219      DO ib_obc = 1, nb_obc
220
221         jpbdta = 1
222         IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters
223 
224
225           !! 1. Read parameters from namelist
226           !! 2. Work out global size of boundary data arrays nblendta and jpbdta
227
228
229         ELSE            ! Read size of arrays in boundary coordinates file.
230
231
232            CALL iom_open( cn_coords_file(ib_obc), inum )
233            jpbdta = 1
234            DO igrd = 1, jpbgrd
235               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
236               nblendta(igrd,ib_obc) = kdimsz(1)
237               jpbdta = MAX(jpbdta, kdimsz(1))
238            ENDDO
239
240         ENDIF
241
242      ENDDO
243
244      ! Allocate arrays
245      !---------------
246      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_obc), nbjdta(jpbdta, jpbgrd, nb_obc),    &
247         &      nbrdta(jpbdta, jpbgrd, nb_obc) )
248
249      ALLOCATE( dta_global(jpbdta, 1, jpk) )
250
251      ! Calculate global boundary index arrays or read in from file
252      !------------------------------------------------------------
253      DO ib_obc = 1, nb_obc
254
255         IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters
256
257           !! Calculate global index arrays from namelist parameters
258
259         ELSE            ! Read global index arrays from boundary coordinates file.
260
261            DO igrd = 1, jpbgrd
262               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
263               DO ii = 1,nblendta(igrd,ib_obc)
264                  nbidta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
265               END DO
266               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
267               DO ii = 1,nblendta(igrd,ib_obc)
268                  nbjdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
269               END DO
270               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
271               DO ii = 1,nblendta(igrd,ib_obc)
272                  nbrdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
273               END DO
274
275               ibr_max = MAXVAL( nbrdta(:,igrd,ib_obc) )
276               IF(lwp) WRITE(numout,*)
277               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
278               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_obc)
279               IF (ibr_max < nn_rimwidth(ib_obc))   &
280                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_obc) )
281
282            END DO
283            CALL iom_close( inum )
284
285         ENDIF
286
287      ENDDO 
288
289      ! Work out dimensions of boundary data on each processor
290      ! ------------------------------------------------------
291     
292      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
293      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
294      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
295      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
296
297      DO ib_obc = 1, nb_obc
298         DO igrd = 1, jpbgrd
299            icount  = 0
300            icountr = 0
301            idx_obc(ib_obc)%nblen(igrd)    = 0
302            idx_obc(ib_obc)%nblenrim(igrd) = 0
303            DO ib = 1, nblendta(igrd,ib_obc)
304               ! check if point is in local domain
305               IF(  nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND.   &
306                  & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in       ) THEN
307                  !
308                  icount = icount  + 1
309                  !
310                  IF( nbrdta(ib,igrd,ib_obc) == 1 )   icountr = icountr+1
311               ENDIF
312            ENDDO
313            idx_obc(ib_obc)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
314            idx_obc(ib_obc)%nblen   (igrd) = icount  !: length of boundary data on each proc       
315         ENDDO  ! igrd
316
317         ! Allocate index arrays for this boundary set
318         !--------------------------------------------
319         ilen1 = MAXVAL(idx_obc(ib_obc)%nblen(:))
320         ALLOCATE( idx_obc(ib_obc)%nbi(ilen1,jpbgrd) )
321         ALLOCATE( idx_obc(ib_obc)%nbj(ilen1,jpbgrd) )
322         ALLOCATE( idx_obc(ib_obc)%nbr(ilen1,jpbgrd) )
323         ALLOCATE( idx_obc(ib_obc)%nbmap(ilen1,jpbgrd) )
324         ALLOCATE( idx_obc(ib_obc)%nbw(ilen1,jpbgrd) )
325         ALLOCATE( idx_obc(ib_obc)%flagu(ilen1) )
326         ALLOCATE( idx_obc(ib_obc)%flagv(ilen1) )
327
328         ! Dispatch mapping indices and discrete distances on each processor
329         ! -----------------------------------------------------------------
330
331         DO igrd = 1, jpbgrd
332            icount  = 0
333            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
334            DO ir=1, nn_rimwidth(ib_obc)
335               DO ib = 1, nblendta(igrd,ib_obc)
336                  ! check if point is in local domain and equals ir
337                  IF(  nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND.   &
338                     & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in .AND.   &
339                     & nbrdta(ib,igrd,ib_obc) == ir  ) THEN
340                     !
341                     icount = icount  + 1
342                     idx_obc(ib_obc)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_obc)- mig(1)+1
343                     idx_obc(ib_obc)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_obc)- mjg(1)+1
344                     idx_obc(ib_obc)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_obc)
345                     idx_obc(ib_obc)%nbmap(icount,igrd) = ib
346                  ENDIF
347               ENDDO
348            ENDDO
349         ENDDO 
350
351         ! Compute rim weights for FRS scheme
352         ! ----------------------------------
353         DO igrd = 1, jpbgrd
354            DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
355               nbr => idx_obc(ib_obc)%nbr(ib,igrd)
356               idx_obc(ib_obc)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
357!              idx_obc(ib_obc)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
358!              idx_obc(ib_obc)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
359            END DO
360         END DO
361
362      ENDDO
363
364      ! ------------------------------------------------------
365      ! Initialise masks and find normal/tangential directions
366      ! ------------------------------------------------------
367
368      ! Read global 2D mask at T-points: obctmask
369      ! -----------------------------------------
370      ! obctmask = 1  on the computational domain AND on open boundaries
371      !          = 0  elsewhere   
372 
373      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
374         zmask(         :                ,:) = 0.e0
375         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0         
376      ELSE IF( ln_mask_file ) THEN
377         CALL iom_open( cn_mask_file, inum )
378         CALL iom_get ( inum, jpdom_data, 'obc_msk', zmask(:,:) )
379         CALL iom_close( inum )
380      ELSE
381         zmask(:,:) = 1.e0
382      ENDIF
383
384      DO ij = 1, nlcj      ! Save mask over local domain     
385         DO ii = 1, nlci
386            obctmask(ii,ij) = zmask( mig(ii), mjg(ij) )
387         END DO
388      END DO
389
390      ! Derive mask on U and V grid from mask on T grid
391      obcumask(:,:) = 0.e0
392      obcvmask(:,:) = 0.e0
393      DO ij=1, jpjm1
394         DO ii=1, jpim1
395            obcumask(ii,ij)=obctmask(ii,ij)*obctmask(ii+1, ij )
396            obcvmask(ii,ij)=obctmask(ii,ij)*obctmask(ii  ,ij+1) 
397         END DO
398      END DO
399      CALL lbc_lnk( obcumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( obcvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
400
401
402      ! Mask corrections
403      ! ----------------
404      DO ik = 1, jpkm1
405         DO ij = 1, jpj
406            DO ii = 1, jpi
407               tmask(ii,ij,ik) = tmask(ii,ij,ik) * obctmask(ii,ij)
408               umask(ii,ij,ik) = umask(ii,ij,ik) * obcumask(ii,ij)
409               vmask(ii,ij,ik) = vmask(ii,ij,ik) * obcvmask(ii,ij)
410               bmask(ii,ij)    = bmask(ii,ij)    * obctmask(ii,ij)
411            END DO     
412         END DO
413      END DO
414
415      DO ik = 1, jpkm1
416         DO ij = 2, jpjm1
417            DO ii = 2, jpim1
418               fmask(ii,ij,ik) = fmask(ii,ij,ik) * obctmask(ii,ij  ) * obctmask(ii+1,ij  )   &
419                  &                              * obctmask(ii,ij+1) * obctmask(ii+1,ij+1)
420            END DO     
421         END DO
422      END DO
423
424      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
425      obctmask(:,:) = tmask(:,:,1)
426
427      ! obc masks and bmask are now set to zero on boundary points:
428      igrd = 1       ! In the free surface case, bmask is at T-points
429      DO ib_obc = 1, nb_obc
430        DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)     
431          bmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
432        ENDDO
433      ENDDO
434      !
435      igrd = 1
436      DO ib_obc = 1, nb_obc
437        DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)     
438          obctmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
439        ENDDO
440      ENDDO
441      !
442      igrd = 2
443      DO ib_obc = 1, nb_obc
444        DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
445          obcumask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
446        ENDDO
447      ENDDO
448      !
449      igrd = 3
450      DO ib_obc = 1, nb_obc
451        DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
452          obcvmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
453        ENDDO
454      ENDDO
455
456      ! Lateral boundary conditions
457      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( obctmask(:,:), 'T', 1. )
458      CALL lbc_lnk( obcumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( obcvmask(:,:), 'V', 1. )
459
460      DO ib_obc = 1, nb_obc       ! Indices and directions of rim velocity components
461
462         idx_obc(ib_obc)%flagu(:) = 0.e0
463         idx_obc(ib_obc)%flagv(:) = 0.e0
464         icount = 0 
465
466         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
467         !flagu =  0 : u is tangential
468         !flagu =  1 : u is normal to the boundary and is direction is inward
469 
470         igrd = 2      ! u-component
471         DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 
472            nbi => idx_obc(ib_obc)%nbi(ib,igrd)
473            nbj => idx_obc(ib_obc)%nbj(ib,igrd)
474            zefl = obctmask(nbi  ,nbj)
475            zwfl = obctmask(nbi+1,nbj)
476            IF( zefl + zwfl == 2 ) THEN
477               icount = icount + 1
478            ELSE
479               idx_obc(ib_obc)%flagu(ib)=-zefl+zwfl
480            ENDIF
481         END DO
482
483         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
484         !flagv =  0 : u is tangential
485         !flagv =  1 : u is normal to the boundary and is direction is inward
486
487         igrd = 3      ! v-component
488         DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 
489            nbi => idx_obc(ib_obc)%nbi(ib,igrd)
490            nbj => idx_obc(ib_obc)%nbj(ib,igrd)
491            znfl = obctmask(nbi,nbj  )
492            zsfl = obctmask(nbi,nbj+1)
493            IF( znfl + zsfl == 2 ) THEN
494               icount = icount + 1
495            ELSE
496               idx_obc(ib_obc)%flagv(ib) = -znfl + zsfl
497            END IF
498         END DO
499 
500         IF( icount /= 0 ) THEN
501            IF(lwp) WRITE(numout,*)
502            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
503               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_obc
504            IF(lwp) WRITE(numout,*) ' ========== '
505            IF(lwp) WRITE(numout,*)
506            nstop = nstop + 1
507         ENDIF
508   
509      ENDDO
510
511      ! Compute total lateral surface for volume correction:
512      ! ----------------------------------------------------
513      obcsurftot = 0.e0 
514      IF( ln_vol ) THEN 
515         igrd = 2      ! Lateral surface at U-points
516         DO ib_obc = 1, nb_obc
517            DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
518               nbi => idx_obc(ib_obc)%nbi(ib,igrd)
519               nbj => idx_obc(ib_obc)%nbi(ib,igrd)
520               flagu => idx_obc(ib_obc)%flagu(ib)
521               obcsurftot = obcsurftot + hu     (nbi  , nbj)                           &
522                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
523                  &                    * tmask_i(nbi  , nbj)                           &
524                  &                    * tmask_i(nbi+1, nbj)                   
525            ENDDO
526         ENDDO
527
528         igrd=3 ! Add lateral surface at V-points
529         DO ib_obc = 1, nb_obc
530            DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
531               nbi => idx_obc(ib_obc)%nbi(ib,igrd)
532               nbj => idx_obc(ib_obc)%nbi(ib,igrd)
533               flagv => idx_obc(ib_obc)%flagv(ib)
534               obcsurftot = obcsurftot + hv     (nbi, nbj  )                           &
535                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
536                  &                    * tmask_i(nbi, nbj  )                           &
537                  &                    * tmask_i(nbi, nbj+1)
538            ENDDO
539         ENDDO
540         !
541         IF( lk_mpp )   CALL mpp_sum( obcsurftot )      ! sum over the global domain
542      END IF   
543      !
544      ! Tidy up
545      !--------
546      DEALLOCATE(nbidta, nbjdta, nbrdta)
547
548   END SUBROUTINE obc_init
549
550#else
551   !!---------------------------------------------------------------------------------
552   !!   Dummy module                                   NO open boundaries
553   !!---------------------------------------------------------------------------------
554CONTAINS
555   SUBROUTINE obc_init      ! Dummy routine
556   END SUBROUTINE obc_init
557#endif
558
559   !!=================================================================================
560END MODULE obcini
Note: See TracBrowser for help on using the repository browser.