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

source: trunk/NEMO/OPA_SRC/BDY/bdyini.F90 @ 1528

Last change on this file since 1528 was 1528, checked in by rblod, 15 years ago

Suppress rigid-lid option, see ticket #486

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 18.9 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   !!----------------------------------------------------------------------
11#if defined key_bdy
12   !!----------------------------------------------------------------------
13   !!   'key_bdy'                     Unstructured Open Boundary Conditions
14   !!----------------------------------------------------------------------
15   !!   bdy_init       : Initialization of unstructured open boundaries
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain
19   USE bdy_oce         ! unstructured open boundary conditions
20   USE bdytides        ! tides at open boundaries initialization (tide_init routine)
21   USE in_out_manager  ! I/O units
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp         ! for mpp_sum 
24   USE iom             ! I/O
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   bdy_init   ! routine called by opa.F90
30
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
33   !! $Id$
34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
35   !!---------------------------------------------------------------------------------
36
37CONTAINS
38   
39   SUBROUTINE bdy_init
40      !!----------------------------------------------------------------------
41      !!                 ***  ROUTINE bdy_init  ***
42      !!         
43      !! ** Purpose :   Initialization of the dynamics and tracer fields with
44      !!      unstructured open boundaries.
45      !!
46      !! ** Method  :  Read initialization arrays (mask, indices) to identify
47      !!               an unstructured open boundary
48      !!
49      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
50      !!
51      !!----------------------------------------------------------------------     
52      INTEGER ::   ii, ij, ik, igrd, ib, ir   ! dummy loop indices
53      INTEGER ::   icount, icountr
54      INTEGER ::   ib_len, ibr_max
55      INTEGER ::   iw, ie, is, in 
56      INTEGER ::   inum                 ! temporary logical unit
57      INTEGER ::   id_dummy             ! temporary integers
58      INTEGER ::   igrd_start, igrd_end ! start and end of loops on igrd
59      INTEGER, DIMENSION (2)             ::   kdimsz
60      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
61      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbrdta           ! Discrete distance from rim points
62      REAL(wp) :: zefl, zwfl, znfl, zsfl                       ! temporary scalars
63      REAL(wp) , DIMENSION(jpidta,jpjdta) ::   zmask           ! global domain mask
64      REAL(wp) , DIMENSION(jpbdta,1)      ::   zdta            ! temporary array
65      CHARACTER(LEN=80),DIMENSION(3)      ::   clfile
66      !!
67      NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V,          &
68         &            ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask,                &
69         &            ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs,                    &
70         &            nbdy_dta   , nb_rimwidth  , volbdy
71      !!----------------------------------------------------------------------
72
73      IF(lwp) WRITE(numout,*)
74      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries'
75      IF(lwp) WRITE(numout,*) '~~~~~~~~'
76      !
77      IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,',   &
78           ' and unstructured open boundary condition are not compatible' )
79
80#if defined key_obc
81      CALL ctl_stop( 'Straight open boundaries,',   &
82           ' and unstructured open boundaries are not compatible' )
83#endif
84
85      ! Read namelist parameters
86      ! ---------------------------
87      REWIND( numnam )
88      READ  ( numnam, nambdy )
89
90      ! control prints
91      IF(lwp) WRITE(numout,*) '         nambdy'
92
93      ! Check nbdy_dta value
94      IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta     
95      IF(lwp) WRITE(numout,*) ' '
96      SELECT CASE( nbdy_dta )
97      CASE( 0 )
98        IF(lwp) WRITE(numout,*) '         initial state used for bdy data'       
99      CASE( 1 )
100        IF(lwp) WRITE(numout,*) '         boundary data taken from file'
101      CASE DEFAULT
102        CALL ctl_stop( 'nbdy_dta must be 0 or 1' )
103      END SELECT
104
105      IF(lwp) WRITE(numout,*) ' '
106      IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth
107
108      IF(lwp) WRITE(numout,*) ' '
109      IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy
110
111      IF (ln_bdy_vol) THEN
112        SELECT CASE ( volbdy ) ! Check volbdy value
113        CASE( 1 )
114          IF(lwp) WRITE(numout,*) '         The total volume will be constant'
115        CASE( 0 )
116          IF(lwp) WRITE(numout,*) '         The total volume will vary according'
117          IF(lwp) WRITE(numout,*) '         to the surface E-P flux'
118        CASE DEFAULT
119          CALL ctl_stop( 'volbdy must be 0 or 1' )
120        END SELECT
121      ELSE
122        IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
123        IF(lwp) WRITE(numout,*) ' '
124      ENDIF
125
126      IF (ln_bdy_tides) THEN
127        IF(lwp) WRITE(numout,*) ' '
128        IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries'
129        IF(lwp) WRITE(numout,*) ' '
130      ENDIF
131
132      IF (ln_bdy_dyn_fla) THEN
133        IF(lwp) WRITE(numout,*) ' '
134        IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries'
135        IF(lwp) WRITE(numout,*) ' '
136      ENDIF
137
138      IF (ln_bdy_dyn_frs) THEN
139        IF(lwp) WRITE(numout,*) ' '
140        IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries'
141        IF(lwp) WRITE(numout,*) ' '
142      ENDIF
143
144      IF (ln_bdy_tra_frs) THEN
145        IF(lwp) WRITE(numout,*) ' '
146        IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries'
147        IF(lwp) WRITE(numout,*) ' '
148      ENDIF
149
150      ! Read tides namelist
151      ! ------------------------
152      IF( ln_bdy_tides )   CALL tide_init
153
154      ! Read arrays defining unstructured open boundaries
155      ! -------------------------------------------------
156
157      ! Read global 2D mask at T-points: bdytmask
158      ! *****************************************
159      ! bdytmask = 1  on the computational domain AND on open boundaries
160      !          = 0  elsewhere   
161 
162      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
163         zmask(         :                ,:) = 0.e0
164         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0         
165      ELSE IF ( ln_bdy_mask ) THEN
166         CALL iom_open( filbdy_mask, inum )
167         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
168         CALL iom_close( inum )
169      ELSE
170         zmask(:,:) = 1.e0
171      ENDIF
172
173      ! Save mask over local domain     
174      DO ij = 1, nlcj
175         DO ii = 1, nlci
176            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
177         END DO
178      END DO
179
180      ! Derive mask on U and V grid from mask on T grid
181      bdyumask(:,:) = 0.e0
182      bdyvmask(:,:) = 0.e0
183      DO ij=1, jpjm1
184         DO ii=1, jpim1
185            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
186            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
187         END DO
188      END DO
189
190      ! Lateral boundary conditions
191      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
192      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
193
194      ! Read discrete distance and mapping indices
195      ! ******************************************
196      nbidta(:,:) = 0.e0
197      nbjdta(:,:) = 0.e0
198      nbrdta(:,:) = 0.e0
199
200      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
201         icount = 0
202         ! Define west boundary (from ii=2 to ii=1+nb_rimwidth):
203         DO ir = 1, nb_rimwidth         
204            DO ij = 3, jpjglo-2
205               icount=icount+1
206               nbidta(icount,:) = ir + 1 + (jpizoom-1)
207               nbjdta(icount,:) = ij + (jpjzoom-1) 
208               nbrdta(icount,:) = ir
209            END DO
210         END DO
211
212         ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth):
213         DO ir=1,nb_rimwidth         
214            DO ij=3,jpjglo-2
215               icount=icount+1
216               nbidta(icount,:) = jpiglo-ir + (jpizoom-1)
217               nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points
218               nbjdta(icount,:) = ij + (jpjzoom-1)
219               nbrdta(icount,:) = ir
220            END DO
221         END DO
222           
223      ELSE            ! Read indices and distances in unstructured boundary data files
224
225         IF( ln_bdy_tides ) THEN 
226            ! Read tides input files for preference in case there are
227            ! no bdydata files.
228            clfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'
229            clfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'
230            clfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'
231         ELSE
232            clfile(1) = filbdy_data_T
233            clfile(2) = filbdy_data_U
234            clfile(3) = filbdy_data_V
235         ENDIF         
236
237         ! how many files are we to read in?
238         igrd_start = 1
239         igrd_end   = 3
240         IF(.NOT. ln_bdy_tides ) THEN
241            IF(.NOT. (ln_bdy_dyn_fla) .AND..NOT. (ln_bdy_tra_frs)) THEN
242               ! No T-grid file.
243               igrd_start = 2
244            ELSEIF ( .NOT. ln_bdy_dyn_frs .AND..NOT. ln_bdy_dyn_fla ) THEN
245               ! No U-grid or V-grid file.
246               igrd_end   = 1         
247            ENDIF
248         ENDIF
249
250         DO igrd = igrd_start, igrd_end
251            CALL iom_open( clfile(igrd), inum )
252            id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) 
253            WRITE(numout,*) 'kdimsz : ',kdimsz
254            ib_len = kdimsz(1)
255            IF( ib_len > jpbdta) CALL ctl_stop(          &
256                'Boundary data array in file too long.', &
257                'File :', TRIM(clfile(igrd)),            &
258                'increase parameter jpbdta.' )
259
260            CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) )
261            DO ii = 1,ib_len
262               nbidta(ii,igrd) = INT( zdta(ii,1) )
263            END DO
264            CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) )
265            DO ii = 1,ib_len
266              nbjdta(ii,igrd) = INT( zdta(ii,1) )
267            END DO
268            CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )
269            DO ii = 1,ib_len
270              nbrdta(ii,igrd) = INT( zdta(ii,1) )
271            END DO
272            CALL iom_close( inum )
273
274            ! Check that rimwidth in file is big enough:
275            ibr_max = MAXVAL( nbrdta(:,igrd) )
276            IF(lwp) WRITE(numout,*)
277            IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
278            IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth
279            IF (ibr_max < nb_rimwidth) CALL ctl_stop( &
280                'nb_rimwidth is larger than maximum rimwidth in file' )
281            !
282         END DO
283         !
284      ENDIF 
285
286      ! Dispatch mapping indices and discrete distances on each processor
287      ! *****************************************************************
288     
289      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
290      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
291      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
292      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
293
294      DO igrd = igrd_start, igrd_end
295        icount  = 0
296        icountr = 0
297        nblen(igrd) = 0
298        nblenrim(igrd) = 0
299        nblendta(igrd) = 0
300        DO ir=1, nb_rimwidth
301          DO ib = 1, jpbdta
302          ! check if point is in local domain and equals ir
303            IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.   &
304               & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   &
305               & nbrdta(ib,igrd) == ir  ) THEN
306               !
307               icount = icount  + 1
308               !
309               IF( ir == 1 )   icountr = icountr+1
310                  IF (icount > jpbdim) THEN
311                     IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'
312                     nstop = nstop + 1
313                  ELSE
314                     nbi(icount, igrd)  = nbidta(ib,igrd)- mig(1)+1
315                     nbj(icount, igrd)  = nbjdta(ib,igrd)- mjg(1)+1
316                     nbr(icount, igrd)  = nbrdta(ib,igrd)
317                     nbmap(icount,igrd) = ib
318                  ENDIF           
319               ENDIF
320            END DO
321         END DO
322         nblenrim(igrd) = icountr !: length of rim boundary data on each proc
323         nblen   (igrd) = icount  !: length of boundary data on each proc       
324      END DO 
325
326      ! Compute rim weights
327      ! -------------------
328      DO igrd = igrd_start, igrd_end
329         DO ib = 1, nblen(igrd)
330            ! tanh formulation
331            nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )
332            ! quadratic
333!           nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2
334            ! linear
335!           nbw(ib,igrd) =  FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth)
336         END DO
337      END DO 
338   
339      ! Mask corrections
340      ! ----------------
341      DO ik = 1, jpkm1
342         DO ij = 1, jpj
343            DO ii = 1, jpi
344               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
345               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
346               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
347               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
348            END DO     
349         END DO
350      END DO
351
352      DO ik = 1, jpkm1
353         DO ij = 2, jpjm1
354            DO ii = 2, jpim1
355               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
356                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
357            END DO     
358         END DO
359      END DO
360
361      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
362      bdytmask(:,:) = tmask(:,:,1)
363
364      ! bdy masks and bmask are now set to zero on boundary points:
365      igrd = 1       ! In the free surface case, bmask is at T-points
366      DO ib = 1, nblenrim(igrd)     
367        bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
368      END DO
369      !
370      igrd = 1
371      DO ib = 1, nblenrim(igrd)     
372        bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
373      END DO
374      !
375      igrd = 2
376      DO ib = 1, nblenrim(igrd)
377        bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
378      END DO
379      !
380      igrd = 3
381      DO ib = 1, nblenrim(igrd)
382        bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
383      END DO
384
385      ! Lateral boundary conditions
386      CALL lbc_lnk( fmask        , 'F', 1. )
387      CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
388      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
389      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
390
391      IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN      ! Indices and directions of rim velocity components
392         !
393         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
394         !flagu =  0 : u is tangential
395         !flagu =  1 : u is normal to the boundary and is direction is inward
396         icount = 0 
397         flagu(:) = 0.e0
398 
399         igrd = 2      ! u-component
400         DO ib = 1, nblenrim(igrd) 
401            zefl=bdytmask(nbi(ib,igrd)  , nbj(ib,igrd))
402            zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd))
403            IF( zefl + zwfl ==2 ) THEN
404               icount = icount +1
405            ELSE
406               flagu(ib)=-zefl+zwfl
407            ENDIF
408         END DO
409
410         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
411         !flagv =  0 : u is tangential
412         !flagv =  1 : u is normal to the boundary and is direction is inward
413         flagv(:) = 0.e0
414
415         igrd = 3      ! v-component
416         DO ib = 1, nblenrim(igrd) 
417            znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd))
418            zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1)
419            IF( znfl + zsfl ==2 ) THEN
420               icount = icount + 1
421            ELSE
422               flagv(ib) = -znfl + zsfl
423            END IF
424         END DO
425 
426         IF( icount /= 0 ) THEN
427            IF(lwp) WRITE(numout,*)
428            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
429               ' are not boundary points. Check nbi, nbj, indices.'
430            IF(lwp) WRITE(numout,*) ' ========== '
431            IF(lwp) WRITE(numout,*)
432            nstop = nstop + 1
433         ENDIF
434   
435      ENDIF
436
437      ! Compute total lateral surface for volume correction:
438      ! ----------------------------------------------------
439 
440      bdysurftot = 0.e0 
441      IF( ln_bdy_vol ) THEN 
442         igrd = 2      ! Lateral surface at U-points
443         DO ib = 1, nblenrim(igrd)
444            bdysurftot = bdysurftot + hu     (nbi(ib,igrd)  ,nbj(ib,igrd))                      &
445               &                    * e2u    (nbi(ib,igrd)  ,nbj(ib,igrd)) * ABS( flagu(ib) )   &
446               &                    * tmask_i(nbi(ib,igrd)  ,nbj(ib,igrd))                      &
447               &                    * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd))                   
448         END DO
449
450         igrd=3 ! Add lateral surface at V-points
451         DO ib = 1, nblenrim(igrd)
452            bdysurftot = bdysurftot + hv     (nbi(ib,igrd),nbj(ib,igrd)  )                      &
453               &                    * e1v    (nbi(ib,igrd),nbj(ib,igrd)  ) * ABS( flagv(ib) )   &
454               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)  )                      &
455               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1)
456         END DO
457
458         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
459      END IF   
460
461      ! Initialise bdy data arrays
462      ! --------------------------
463      tbdy(:,:) = 0.e0
464      sbdy(:,:) = 0.e0
465      ubdy(:,:) = 0.e0
466      vbdy(:,:) = 0.e0
467      sshbdy(:) = 0.e0
468      ubtbdy(:) = 0.e0
469      vbtbdy(:) = 0.e0
470
471      ! Read in tidal constituents and adjust for model start time
472      ! ----------------------------------------------------------
473      IF( ln_bdy_tides )   CALL tide_data
474      !
475   END SUBROUTINE bdy_init
476
477#else
478   !!---------------------------------------------------------------------------------
479   !!   Dummy module                                   NO unstructured open boundaries
480   !!---------------------------------------------------------------------------------
481CONTAINS
482   SUBROUTINE bdy_init      ! Dummy routine
483   END SUBROUTINE bdy_init
484#endif
485
486   !!=================================================================================
487END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.