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.
domzgr.F90 in trunk/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMO/OPA_SRC/DOM/domzgr.F90 @ 450

Last change on this file since 450 was 450, checked in by opalod, 18 years ago

nemo_v1_bugfix_042:RB: - initialization of bathy_level for partial steps

  • set eiv coeff. to 0 when aeiv0 = 0
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.4 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_zgr     : defined the ocean vertical coordinate system
9   !!       zgr_bat      : bathymetry fields (levels and meters)
10   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
11   !!       zgr_bat_ctl  : check the bathymetry files
12   !!       zgr_z        : reference z-coordinate
13   !!       zgr_zps      : z-coordinate with partial steps
14   !!       zgr_s        : s-coordinate
15   !!---------------------------------------------------------------------
16   !! * Modules used
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distributed memory computing library
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE closea
23   USE solisl
24   USE ini1d           ! initialization of the 1D configuration
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC dom_zgr        ! called by dom_init.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!   OPA 9.0 , LOCEAN-IPSL (2005)
37   !! $Header$
38   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
39   !!----------------------------------------------------------------------
40
41CONTAINS       
42
43   SUBROUTINE dom_zgr
44      !!----------------------------------------------------------------------
45      !!                ***  ROUTINE dom_zgr  ***
46      !!                   
47      !! ** Purpose :  set the depth of model levels and the resulting
48      !!      vertical scale factors.
49      !!
50      !! ** Method  reference vertical coordinate
51      !!        Z-coordinates : The depth of model levels is defined
52      !!      from an analytical function the derivative of which gives
53      !!      the vertical scale factors.
54      !!      both depth and scale factors only depend on k (1d arrays).
55      !!              w-level: gdepw  = fsdep(k)
56      !!                       e3w(k) = dk(fsdep)(k)     = fse3(k)
57      !!              t-level: gdept  = fsdep(k+0.5)
58      !!                       e3t(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
59      !!
60      !! ** Action : - gdept, gdepw : depth of T- and W-point (m)
61      !!             -  e3t, e3w    : scale factors at T- and W-levels (m)
62      !!
63      !! Reference :
64      !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
65      !!
66      !! History :
67      !!   9.0  !  03-08  (G. Madec)  original code
68      !!----------------------------------------------------------------------
69      INTEGER ::   ioptio = 0      ! temporary integer
70      !!----------------------------------------------------------------------
71
72      ! Check Vertical coordinate options
73      ! ---------------------------------
74      ioptio = 0
75      IF( lk_sco ) THEN
76         IF(lwp) WRITE(numout,*)
77         IF(lwp) WRITE(numout,*) 'dom_zgr : s-coordinate'
78         IF(lwp) WRITE(numout,*) '~~~~~~~'
79         ioptio = ioptio + 1
80      ENDIF
81      IF( lk_zps ) THEN
82         IF(lwp) WRITE(numout,*)
83         IF(lwp) WRITE(numout,*) 'dom_zgr : z-coordinate with partial steps'
84         IF(lwp) WRITE(numout,*) '~~~~~~~'
85         ioptio = ioptio + 1
86      ENDIF
87      IF( ioptio == 0 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'dom_zgr : z-coordinate'
90         IF(lwp) WRITE(numout,*) '~~~~~~~'
91      ENDIF
92
93      IF ( ioptio > 1 ) THEN
94          IF(lwp) WRITE(numout,cform_err)
95          IF(lwp) WRITE(numout,*) ' several vertical coordinate options used'
96          nstop = nstop + 1
97      ENDIF
98
99      ! Build the vertical coordinate system
100      ! ------------------------------------
101
102      CALL zgr_z                       ! Reference z-coordinate system
103
104      CALL zgr_bat                     ! Bathymetry fields (levels and meters)
105
106      CALL zgr_zps                     ! Partial step z-coordinate
107
108      CALL zgr_s                       ! s-coordinate
109
110   END SUBROUTINE dom_zgr
111
112
113   SUBROUTINE zgr_z
114      !!----------------------------------------------------------------------
115      !!                   ***  ROUTINE zgr_z  ***
116      !!                   
117      !! ** Purpose :   set the depth of model levels and the resulting
118      !!      vertical scale factors.
119      !!
120      !! ** Method  :   z-coordinate system (use in all type of coordinate)
121      !!        The depth of model levels is defined from an analytical
122      !!      function the derivative of which gives the scale factors.
123      !!        both depth and scale factors only depend on k (1d arrays).
124      !!              w-level: gdepw  = fsdep(k)
125      !!                       e3w(k) = dk(fsdep)(k)     = fse3(k)
126      !!              t-level: gdept  = fsdep(k+0.5)
127      !!                       e3t(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
128      !!
129      !! ** Action  : - gdept, gdepw : depth of T- and W-point (m)
130      !!              -  e3t, e3w    : scale factors at T- and W-levels (m)
131      !!
132      !! Reference :
133      !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
134      !!
135      !! History :
136      !!   9.0  !  03-08  (G. Madec)  F90: Free form and module
137      !!----------------------------------------------------------------------
138      !! * Local declarations
139      INTEGER  ::   jk                     ! dummy loop indices
140      REAL(wp) ::   zt, zw                 ! temporary scalars
141      REAL(wp) ::   &
142      zsur , za0, za1, zkth, zacr,      &  ! Values set from parameters in
143      zdzmin, zhmax                        ! par_CONFIG_Rxx.h90
144      !!----------------------------------------------------------------------
145
146      ! Set variables from parameters
147      ! ------------------------------
148       zkth = ppkth       ;   zacr = ppacr
149       zdzmin = ppdzmin   ;   zhmax = pphmax
150
151      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
152      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
153      !
154       IF(  ppa1  == pp_to_be_computed  .AND.  &
155         &  ppa0  == pp_to_be_computed  .AND.  &
156         &  ppsur == pp_to_be_computed           ) THEN
157         za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) )          &
158             / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) &
159             &                         *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
160             &                             - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
161
162         za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
163         zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
164
165       ELSE
166         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
167       ENDIF
168
169
170      ! Parameter print
171      ! ---------------
172      IF(lwp) THEN
173         WRITE(numout,*)
174         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
175         WRITE(numout,*) '    ~~~~~~~'
176         IF (  ppkth == 0. ) THEN             
177              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
178              WRITE(numout,*) '            Total depth    :', zhmax
179              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
180         ELSE
181            IF ( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
182               WRITE(numout,*) '         zsur, za0, za1 computed from '
183               WRITE(numout,*) '                 zdzmin = ', zdzmin
184               WRITE(numout,*) '                 zhmax  = ', zhmax
185            ENDIF
186            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
187            WRITE(numout,*) '                 zsur = ', zsur
188            WRITE(numout,*) '                 za0  = ', za0
189            WRITE(numout,*) '                 za1  = ', za1
190            WRITE(numout,*) '                 zkth = ', zkth
191            WRITE(numout,*) '                 zacr = ', zacr
192         ENDIF
193      ENDIF
194
195
196      ! Reference z-coordinate (depth - scale factor at T- and W-points)
197      ! ======================
198      IF (  ppkth == 0. ) THEN            !  uniform vertical grid       
199
200         za1 = zhmax/FLOAT(jpk-1) 
201         DO jk = 1, jpk
202            zw = FLOAT( jk )
203            zt = FLOAT( jk ) + 0.5
204            gdepw(jk) = ( zw - 1 ) * za1
205            gdept(jk) = ( zt - 1 ) * za1
206            e3w  (jk) =  za1
207            e3t  (jk) =  za1
208         END DO
209
210      ELSE
211
212         DO jk = 1, jpk
213            zw = FLOAT( jk )
214            zt = FLOAT( jk ) + 0.5
215            gdepw(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) )  )
216            gdept(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) )  )
217            e3w  (jk) =          za0      + za1        * TANH(      (zw-zkth)/zacr   )
218            e3t  (jk) =          za0      + za1        * TANH(      (zt-zkth)/zacr   )
219         END DO
220         gdepw(1) = 0.e0   ! force first w-level to be exactly at zero
221
222      ENDIF
223
224      ! Control and  print
225      ! ==================
226
227      IF(lwp) THEN
228         WRITE(numout,*)
229         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
230         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
231         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept(jk), gdepw(jk), e3t(jk), e3w(jk), jk = 1, jpk )
232      ENDIF
233
234      DO jk = 1, jpk
235         IF( e3w(jk) <= 0. .OR. e3t(jk) <= 0. ) THEN
236            IF(lwp) WRITE(numout,cform_err)
237            IF(lwp) WRITE(numout,*) ' e3w or e3t =< 0 '
238            nstop = nstop + 1
239         ENDIF
240         IF( gdepw(jk) < 0. .OR. gdept(jk) < 0.) THEN
241            IF(lwp) WRITE(numout,cform_err)
242            IF(lwp) WRITE(numout,*) ' gdepw or gdept < 0 '
243            nstop = nstop + 1
244         ENDIF
245      END DO
246
247   END SUBROUTINE zgr_z
248
249
250   SUBROUTINE zgr_bat
251      !!----------------------------------------------------------------------
252      !!                    ***  ROUTINE zgr_bat  ***
253      !!
254      !! ** Purpose :   set bathymetry both in levels and meters
255      !!
256      !! ** Method  :   read or define mbathy and bathy arrays
257      !!       * level bathymetry:
258      !!      The ocean basin geometry is given by a two-dimensional array,
259      !!      mbathy, which is defined as follow :
260      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
261      !!                              at t-point (ji,jj).
262      !!                            = 0  over the continental t-point.
263      !!                            = -n over the nth island t-point.
264      !!      The array mbathy is checked to verified its consistency with
265      !!      model option. in particular:
266      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
267      !!                  along closed boundary.
268      !!            mbathy must be cyclic IF jperio=1.
269      !!            mbathy must be lower or equal to jpk-1.
270      !!            isolated ocean grid points are suppressed from mbathy
271      !!                  since they are only connected to remaining
272      !!                  ocean through vertical diffusion.
273      !!      ntopo=-1 :   rectangular channel or bassin with a bump
274      !!      ntopo= 0 :   flat rectangular channel or basin
275      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
276      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
277      !!      C A U T I O N : mbathy will be modified during the initializa-
278      !!      tion phase to become the number of non-zero w-levels of a water
279      !!      column, with a minimum value of 1.
280      !!
281      !! ** Action  : - mbathy: level bathymetry (in level index)
282      !!              - bathy : meter bathymetry (in meters)
283      !!
284      !! History :
285      !!   9.0  !  03-08  (G. Madec)  Original code
286      !!----------------------------------------------------------------------
287      !! * Modules used
288      USE ioipsl
289
290      !! * Local declarations
291      CHARACTER (len=18) ::   clname    ! temporary characters
292      LOGICAL ::   llbon                ! check the existence of bathy files
293      INTEGER ::   ji, jj, jl, jk       ! dummy loop indices
294      INTEGER ::   inum = 11            ! temporary logical unit
295      INTEGER  ::   &
296         ipi, ipj, ipk,              &  ! temporary integers
297         itime,                      &  !    "          "
298         ii_bump, ij_bump               ! bump center position
299      INTEGER, DIMENSION (1) ::   istep
300      INTEGER , DIMENSION(jpidta,jpjdta) ::   &
301         idta                           ! global domain integer data
302      REAL(wp) ::   &
303         r_bump, h_bump, h_oce,      &  ! bump characteristics
304         zi, zj, zdate0, zdt            ! temporary scalars
305      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
306         zlamt, zphit,               &  ! temporary workspace (NetCDF read)
307         zdta                           ! global domain scalar data
308      REAL(wp), DIMENSION(jpk) ::   &
309         zdept                          ! temporary workspace (NetCDF read)
310      !!----------------------------------------------------------------------
311
312      IF(lwp) WRITE(numout,*)
313      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
314      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
315
316      ! ========================================
317      ! global domain level and meter bathymetry (idta,zdta)
318      ! ========================================
319      !                                               ! =============== !
320      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          ! defined by hand !
321         !                                            ! =============== !
322
323         IF( ntopo == 0 ) THEN                        ! flat basin
324
325            IF(lwp) WRITE(numout,*)
326            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
327
328            idta(:,:) = jpkm1                            ! flat basin
329            zdta(:,:) = gdepw(jpk)
330
331         ELSE                                         ! bump
332            IF(lwp) WRITE(numout,*)
333            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
334
335            ii_bump = jpidta / 3 + 3       ! i-index of the bump center
336            ij_bump = jpjdta / 2           ! j-index of the bump center
337            r_bump  =    6              ! bump radius (index)       
338            h_bump  =  240.e0           ! bump height (meters)
339            h_oce   = gdepw(jpk)        ! background ocean depth (meters)
340            IF(lwp) WRITE(numout,*) '            bump characteristics: '
341            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
342            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
343            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
344            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
345            ! zdta :
346            DO jj = 1, jpjdta
347               DO ji = 1, jpidta
348                  zi = FLOAT( ji - ii_bump ) / r_bump     
349                  zj = FLOAT( jj - ij_bump ) / r_bump       
350                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
351               END DO
352            END DO
353            ! idta :
354            idta(:,:) = jpkm1
355            DO jk = 1, jpkm1
356               DO jj = 1, jpjdta
357                  DO ji = 1, jpidta
358                     IF( gdept(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept(jk+1) )   idta(ji,jj) = jk
359                  END DO
360               END DO
361            END DO
362         ENDIF
363
364         ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
365         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
366            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0
367            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0.e0
368         ELSEIF( jperio == 2 ) THEN
369            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
370            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0
371            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0
372            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0
373         ELSE
374            idta( :    , 1    ) = 0                 ;      zdta( :    , 1    ) =  0.e0
375            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0
376            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0
377            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0
378         ENDIF
379
380         !  EEL R5 configuration with east and west open boundaries.
381         !  Two rows of zeroes are needed at the south and north for OBCs
382         !  This is for compatibility with the rigid lid option.
383         
384         IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
385            idta( : , 2      ) = 0                 ;      zdta( : , 2      ) =  0.e0
386            idta( : ,jpjdta-1) = 0                 ;      zdta( : ,jpjdta-1) =  0.e0
387         ENDIF
388
389         !                                            ! =============== !
390      ELSEIF( ntopo == 1 ) THEN                       !   read in file  !
391         !                                            ! =============== !
392         IF( lk_zco ) THEN
393            clname = 'bathy_level.nc'                       ! Level bathymetry
394#if defined key_agrif
395            IF( .NOT. Agrif_Root() ) THEN
396               clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
397            ENDIF
398#endif         
399            INQUIRE( FILE=clname, EXIST=llbon )
400            IF( llbon ) THEN
401               IF(lwp) WRITE(numout,*)
402               IF(lwp) WRITE(numout,*) '         read level bathymetry in ', clname
403               IF(lwp) WRITE(numout,*)
404               itime = 1
405               ipi = jpidta
406               ipj = jpjdta
407               ipk = 1
408               zdt = rdt
409               CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE.,   &
410                              ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum )
411               CALL flinget( inum, 'Bathy_level', jpidta, jpjdta, 1,   &
412                             itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) )
413               idta(:,:) = zdta(:,:)
414               CALL flinclo( inum )
415
416            ELSE
417               IF(lwp) WRITE(numout,cform_err)
418               IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
419               nstop = nstop + 1
420            ENDIF   
421 
422         ELSEIF( lk_zps ) THEN
423            clname = 'bathy_meter.nc'                       ! meter bathymetry
424#if defined key_agrif
425            IF( .NOT. Agrif_Root() ) THEN
426               clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
427            ENDIF
428#endif       
429            INQUIRE( FILE=clname, EXIST=llbon )
430            IF( llbon ) THEN
431               IF(lwp) WRITE(numout,*)
432               IF(lwp) WRITE(numout,*) '         read meter bathymetry in ', clname
433               IF(lwp) WRITE(numout,*)
434               itime = 1
435               ipi = jpidta
436               ipj = jpjdta
437               ipk = 1
438               zdt = rdt
439               CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE.,   &   
440                              ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum )
441               CALL flinget( inum, 'Bathymetry', jpidta, jpjdta, 1,   &
442                             itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 
443               CALL flinclo( inum )
444               idta(:,:) = jpkm1      ! initialisation
445            ELSE
446               IF(lwp) WRITE(numout,cform_err)       
447               IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
448               nstop = nstop + 1
449            ENDIF
450         ENDIF
451         !                                            ! =============== !
452      ELSE                                            !      error      !
453         !                                            ! =============== !
454         IF(lwp) WRITE(numout,cform_err)
455         IF(lwp) WRITE(numout,*) '          parameter , ntopo = ', ntopo
456         nstop = nstop + 1
457      ENDIF
458
459
460      ! =======================================
461      ! local domain level and meter bathymetry (mbathy,bathy)
462      ! =======================================
463
464      mbathy(:,:) = 0                                 ! set to zero extra halo points
465      bathy (:,:) = 0.e0                              ! (require for mpp case)
466
467      DO jj = 1, nlcj                                 ! interior values
468         DO ji = 1, nlci
469            mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
470            bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
471         END DO
472      END DO
473
474      ! =======================
475      ! NO closed seas or lakes
476      ! =======================
477
478      IF( nclosea == 0 ) THEN
479         DO jl = 1, jpncs
480            DO jj = ncsj1(jl), ncsj2(jl)
481               DO ji = ncsi1(jl), ncsi2(jl)
482                  mbathy(ji,jj) = 0                   ! suppress closed seas
483                  bathy (ji,jj) = 0.e0                ! and lakes
484               END DO
485            END DO
486         END DO
487      ENDIF
488
489      ! ===========
490      ! Zoom domain
491      ! ===========
492
493      IF( lzoom )   CALL zgr_bat_zoom
494
495      ! ================
496      ! Bathymetry check
497      ! ================
498
499      CALL zgr_bat_ctl
500
501   END SUBROUTINE zgr_bat
502
503
504   SUBROUTINE zgr_bat_zoom
505      !!----------------------------------------------------------------------
506      !!                    ***  ROUTINE zgr_bat_zoom  ***
507      !!
508      !! ** Purpose : - Close zoom domain boundary if necessary
509      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
510      !!
511      !! ** Method  :
512      !!
513      !! ** Action  : - update mbathy: level bathymetry (in level index)
514      !!
515      !! History :
516      !!   9.0  !  03-08  (G. Madec)  Original code
517      !!----------------------------------------------------------------------
518      !! * Local variables
519      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
520      !!----------------------------------------------------------------------
521
522      IF(lwp) WRITE(numout,*)
523      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
524      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
525
526      ! Zoom domain
527      ! ===========
528
529      ! Forced closed boundary if required
530      IF( lzoom_w )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0
531      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
532      IF( lzoom_e )   mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
533      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
534
535      ! Configuration specific domain modifications
536      ! (here, ORCA arctic configuration: suppress Med Sea)
537      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
538         SELECT CASE ( jp_cfg )
539         !                                        ! =======================
540         CASE ( 2 )                               !  ORCA_R2 configuration
541            !                                     ! =======================
542            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
543
544            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
545            ij0 =  98   ;   ij1 = 110
546            !                                     ! =======================
547         CASE ( 05 )                              !  ORCA_R05 configuration
548            !                                     ! =======================
549            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
550 
551            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
552            ij0 = 314   ;   ij1 = 370 
553
554         END SELECT
555         !
556         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
557         !
558      ENDIF
559
560   END SUBROUTINE zgr_bat_zoom
561
562
563   SUBROUTINE zgr_bat_ctl
564      !!----------------------------------------------------------------------
565      !!                    ***  ROUTINE zgr_bat_ctl  ***
566      !!
567      !! ** Purpose :   check the bathymetry in levels
568      !!
569      !! ** Method  :   The array mbathy is checked to verified its consistency
570      !!      with the model options. in particular:
571      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
572      !!                  along closed boundary.
573      !!            mbathy must be cyclic IF jperio=1.
574      !!            mbathy must be lower or equal to jpk-1.
575      !!            isolated ocean grid points are suppressed from mbathy
576      !!                  since they are only connected to remaining
577      !!                  ocean through vertical diffusion.
578      !!      C A U T I O N : mbathy will be modified during the initializa-
579      !!      tion phase to become the number of non-zero w-levels of a water
580      !!      column, with a minimum value of 1.
581      !!
582      !! ** Action  : - update mbathy: level bathymetry (in level index)
583      !!              - update bathy : meter bathymetry (in meters)
584      !!
585      !! History :
586      !!   9.0  !  03-08  (G. Madec)  Original code
587      !!----------------------------------------------------------------------
588      !! * Local declarations
589      INTEGER ::   ji, jj, jl           ! dummy loop indices
590      INTEGER ::   &
591         icompt, ibtest, ikmax          ! temporary integers
592      REAL(wp), DIMENSION(jpi,jpj) ::   &
593         zbathy                         ! temporary workspace
594      !!----------------------------------------------------------------------
595
596      IF(lwp) WRITE(numout,*)
597      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
598      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
599
600      ! ================
601      ! Bathymetry check
602      ! ================
603
604      IF( .NOT. lk_cfg_1d )   THEN
605
606         ! Suppress isolated ocean grid points
607
608         IF(lwp) WRITE(numout,*)
609         IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
610         IF(lwp) WRITE(numout,*)'                   -----------------------------------'
611
612         icompt = 0
613
614         DO jl = 1, 2
615
616            IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
617               mbathy( 1 ,:) = mbathy(jpim1,:)
618               mbathy(jpi,:) = mbathy(  2  ,:)
619            ENDIF
620            DO jj = 2, jpjm1
621               DO ji = 2, jpim1
622                  ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   &
623                     mbathy(ji,jj-1),mbathy(ji,jj+1) )
624                  IF( ibtest < mbathy(ji,jj) ) THEN
625                     IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
626                        'grid-point (i,j) =  ',ji,jj,' is changed from ',   &
627                        mbathy(ji,jj),' to ', ibtest
628                     mbathy(ji,jj) = ibtest
629                     icompt = icompt + 1
630                  ENDIF
631               END DO
632            END DO
633
634         END DO
635         IF( icompt == 0 ) THEN
636            IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
637         ELSE
638            IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
639         ENDIF
640         IF( lk_mpp ) THEN
641            zbathy(:,:) = FLOAT( mbathy(:,:) )
642            CALL lbc_lnk( zbathy, 'T', 1. )
643            mbathy(:,:) = INT( zbathy(:,:) )
644         ENDIF
645
646         ! 3.2 East-west cyclic boundary conditions
647
648         IF( nperio == 0 ) THEN
649            IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   &
650               ' boundary: nperio = ', nperio
651            IF( lk_mpp ) THEN
652               IF( nbondi == -1 .OR. nbondi == 2 ) THEN
653                  IF( jperio /= 1 )   mbathy(1,:) = 0
654               ENDIF
655               IF( nbondi == 1 .OR. nbondi == 2 ) THEN
656                  IF( jperio /= 1 )   mbathy(nlci,:) = 0
657               ENDIF
658            ELSE
659               mbathy( 1 ,:) = 0
660               mbathy(jpi,:) = 0
661            ENDIF
662         ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
663            IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   &
664               ' on mbathy: nperio = ', nperio
665            mbathy( 1 ,:) = mbathy(jpim1,:)
666            mbathy(jpi,:) = mbathy(  2  ,:)
667         ELSEIF( nperio == 2 ) THEN
668            IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   &
669               ' on mbathy: nperio = ', nperio
670         ELSE
671            IF(lwp) WRITE(numout,*) '    e r r o r'
672            IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
673            !         STOP 'dom_mba'
674         ENDIF
675
676         ! Set to zero mbathy over islands if necessary  (lk_isl=F)
677         IF( .NOT. lk_isl ) THEN    ! No island
678            IF(lwp) WRITE(numout,*)
679            IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
680            IF(lwp) WRITE(numout,*) '         ----------------------------'
681
682            mbathy(:,:) = MAX( 0, mbathy(:,:) )
683
684            !  Boundary condition on mbathy
685            IF( .NOT.lk_mpp ) THEN 
686               
687               !!bug ???  y reflechir!
688               !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
689               
690               zbathy(:,:) = FLOAT( mbathy(:,:) )
691               CALL lbc_lnk( zbathy, 'T', 1. )
692               mbathy(:,:) = INT( zbathy(:,:) )
693            ENDIF
694
695         ENDIF
696
697      ENDIF
698
699      ! Number of ocean level inferior or equal to jpkm1
700
701      ikmax = 0
702      DO jj = 1, jpj
703         DO ji = 1, jpi
704            ikmax = MAX( ikmax, mbathy(ji,jj) )
705         END DO
706      END DO
707      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
708
709      IF( ikmax > jpkm1 ) THEN
710         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
711         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
712      ELSE IF( ikmax < jpkm1 ) THEN
713         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
714         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
715      ENDIF
716
717      IF( lwp .AND. nprint == 1 ) THEN
718         WRITE(numout,*)
719         WRITE(numout,*) ' bathymetric field '
720         WRITE(numout,*) ' ------------------'
721         WRITE(numout,*) ' number of non-zero T-levels '
722         CALL prihin( mbathy, jpi, jpj, 1, jpi,   &
723                      1     , 1  , jpj, 1, 3  ,   &
724                      numout )
725         WRITE(numout,*)
726      ENDIF
727
728   END SUBROUTINE zgr_bat_ctl
729
730
731#  include "domzgr_zps.h90"
732
733
734#  include "domzgr_s.h90"
735
736
737   !!======================================================================
738END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.