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 @ 911

Last change on this file since 911 was 911, checked in by ctlod, 16 years ago

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