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

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

nemo_v1_bugfix_019 : CT : allow to read either the levels (lk_zco) either the meters (lk_zps) bathymetry
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.0 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=15) ::   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            INQUIRE( FILE=clname, EXIST=llbon )
395            IF( llbon ) THEN
396               IF(lwp) WRITE(numout,*)
397               IF(lwp) WRITE(numout,*) '         read level bathymetry in ', clname
398               IF(lwp) WRITE(numout,*)
399               itime = 1
400               ipi = jpidta
401               ipj = jpjdta
402               ipk = 1
403               zdt = rdt
404               CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE.,   &
405                              ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum )
406               CALL flinget( inum, 'Bathy_level', jpidta, jpjdta, 1,   &
407                             itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) )
408               idta(:,:) = zdta(:,:)
409               CALL flinclo( inum )
410
411            ELSE
412               IF(lwp) WRITE(numout,cform_err)
413               IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
414               nstop = nstop + 1
415            ENDIF   
416 
417         ELSEIF( lk_zps ) THEN
418            clname = 'bathy_meter.nc'                       ! meter bathymetry
419            INQUIRE( FILE=clname, EXIST=llbon )
420            IF( llbon ) THEN
421               IF(lwp) WRITE(numout,*)
422               IF(lwp) WRITE(numout,*) '         read meter bathymetry in ', clname
423               IF(lwp) WRITE(numout,*)
424               itime = 1
425               ipi = jpidta
426               ipj = jpjdta
427               ipk = 1
428               zdt = rdt
429               CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE.,   &   
430                              ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum )
431               CALL flinget( inum, 'Bathymetry', jpidta, jpjdta, 1,   &
432                             itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 
433               CALL flinclo( inum )
434            ELSE
435               IF(lwp) WRITE(numout,cform_err)       
436               IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
437               nstop = nstop + 1
438            ENDIF
439         ENDIF
440         !                                            ! =============== !
441      ELSE                                            !      error      !
442         !                                            ! =============== !
443         IF(lwp) WRITE(numout,cform_err)
444         IF(lwp) WRITE(numout,*) '          parameter , ntopo = ', ntopo
445         nstop = nstop + 1
446      ENDIF
447
448
449      ! =======================================
450      ! local domain level and meter bathymetry (mbathy,bathy)
451      ! =======================================
452
453      mbathy(:,:) = 0                                 ! set to zero extra halo points
454      bathy (:,:) = 0.e0                              ! (require for mpp case)
455
456      DO jj = 1, nlcj                                 ! interior values
457         DO ji = 1, nlci
458            mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
459            bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
460         END DO
461      END DO
462
463      ! =======================
464      ! NO closed seas or lakes
465      ! =======================
466
467      IF( nclosea == 0 ) THEN
468         DO jl = 1, jpncs
469            DO jj = ncsj1(jl), ncsj2(jl)
470               DO ji = ncsi1(jl), ncsi2(jl)
471                  mbathy(ji,jj) = 0                   ! suppress closed seas
472                  bathy (ji,jj) = 0.e0                ! and lakes
473               END DO
474            END DO
475         END DO
476      ENDIF
477
478      ! ===========
479      ! Zoom domain
480      ! ===========
481
482      IF( lzoom )   CALL zgr_bat_zoom
483
484      ! ================
485      ! Bathymetry check
486      ! ================
487
488      CALL zgr_bat_ctl
489
490   END SUBROUTINE zgr_bat
491
492
493   SUBROUTINE zgr_bat_zoom
494      !!----------------------------------------------------------------------
495      !!                    ***  ROUTINE zgr_bat_zoom  ***
496      !!
497      !! ** Purpose : - Close zoom domain boundary if necessary
498      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
499      !!
500      !! ** Method  :
501      !!
502      !! ** Action  : - update mbathy: level bathymetry (in level index)
503      !!
504      !! History :
505      !!   9.0  !  03-08  (G. Madec)  Original code
506      !!----------------------------------------------------------------------
507      !! * Local variables
508      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
509      !!----------------------------------------------------------------------
510
511      IF(lwp) WRITE(numout,*)
512      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
513      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
514
515      ! Zoom domain
516      ! ===========
517
518      ! Forced closed boundary if required
519      IF( lzoom_w )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0
520      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
521      IF( lzoom_e )   mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
522      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
523
524      ! Configuration specific domain modifications
525      ! (here, ORCA arctic configuration: suppress Med Sea)
526      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
527         SELECT CASE ( jp_cfg )
528         !                                        ! =======================
529         CASE ( 2 )                               !  ORCA_R2 configuration
530            !                                     ! =======================
531            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
532
533            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
534            ij0 =  98   ;   ij1 = 110
535            !                                     ! =======================
536         CASE ( 05 )                              !  ORCA_R05 configuration
537            !                                     ! =======================
538            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
539 
540            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
541            ij0 = 314   ;   ij1 = 370 
542
543         END SELECT
544         !
545         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
546         !
547      ENDIF
548
549   END SUBROUTINE zgr_bat_zoom
550
551
552   SUBROUTINE zgr_bat_ctl
553      !!----------------------------------------------------------------------
554      !!                    ***  ROUTINE zgr_bat_ctl  ***
555      !!
556      !! ** Purpose :   check the bathymetry in levels
557      !!
558      !! ** Method  :   The array mbathy is checked to verified its consistency
559      !!      with the model options. in particular:
560      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
561      !!                  along closed boundary.
562      !!            mbathy must be cyclic IF jperio=1.
563      !!            mbathy must be lower or equal to jpk-1.
564      !!            isolated ocean grid points are suppressed from mbathy
565      !!                  since they are only connected to remaining
566      !!                  ocean through vertical diffusion.
567      !!      C A U T I O N : mbathy will be modified during the initializa-
568      !!      tion phase to become the number of non-zero w-levels of a water
569      !!      column, with a minimum value of 1.
570      !!
571      !! ** Action  : - update mbathy: level bathymetry (in level index)
572      !!              - update bathy : meter bathymetry (in meters)
573      !!
574      !! History :
575      !!   9.0  !  03-08  (G. Madec)  Original code
576      !!----------------------------------------------------------------------
577      !! * Local declarations
578      INTEGER ::   ji, jj, jl           ! dummy loop indices
579      INTEGER ::   &
580         icompt, ibtest, ikmax          ! temporary integers
581      REAL(wp), DIMENSION(jpi,jpj) ::   &
582         zbathy                         ! temporary workspace
583      !!----------------------------------------------------------------------
584
585      IF(lwp) WRITE(numout,*)
586      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
587      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
588
589      ! ================
590      ! Bathymetry check
591      ! ================
592
593      IF( .NOT. lk_cfg_1d )   THEN
594
595         ! Suppress isolated ocean grid points
596
597         IF(lwp) WRITE(numout,*)
598         IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
599         IF(lwp) WRITE(numout,*)'                   -----------------------------------'
600
601         icompt = 0
602
603         DO jl = 1, 2
604
605            IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
606               mbathy( 1 ,:) = mbathy(jpim1,:)
607               mbathy(jpi,:) = mbathy(  2  ,:)
608            ENDIF
609            DO jj = 2, jpjm1
610               DO ji = 2, jpim1
611                  ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   &
612                     mbathy(ji,jj-1),mbathy(ji,jj+1) )
613                  IF( ibtest < mbathy(ji,jj) ) THEN
614                     IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
615                        'grid-point (i,j) =  ',ji,jj,' is changed from ',   &
616                        mbathy(ji,jj),' to ', ibtest
617                     mbathy(ji,jj) = ibtest
618                     icompt = icompt + 1
619                  ENDIF
620               END DO
621            END DO
622
623         END DO
624         IF( icompt == 0 ) THEN
625            IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
626         ELSE
627            IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
628         ENDIF
629         IF( lk_mpp ) THEN
630            zbathy(:,:) = FLOAT( mbathy(:,:) )
631            CALL lbc_lnk( zbathy, 'T', 1. )
632            mbathy(:,:) = INT( zbathy(:,:) )
633         ENDIF
634
635         ! 3.2 East-west cyclic boundary conditions
636
637         IF( nperio == 0 ) THEN
638            IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   &
639               ' boundary: nperio = ', nperio
640            IF( lk_mpp ) THEN
641               IF( nbondi == -1 .OR. nbondi == 2 ) THEN
642                  IF( jperio /= 1 )   mbathy(1,:) = 0
643               ENDIF
644               IF( nbondi == 1 .OR. nbondi == 2 ) THEN
645                  IF( jperio /= 1 )   mbathy(nlci,:) = 0
646               ENDIF
647            ELSE
648               mbathy( 1 ,:) = 0
649               mbathy(jpi,:) = 0
650            ENDIF
651         ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
652            IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   &
653               ' on mbathy: nperio = ', nperio
654            mbathy( 1 ,:) = mbathy(jpim1,:)
655            mbathy(jpi,:) = mbathy(  2  ,:)
656         ELSEIF( nperio == 2 ) THEN
657            IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   &
658               ' on mbathy: nperio = ', nperio
659         ELSE
660            IF(lwp) WRITE(numout,*) '    e r r o r'
661            IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
662            !         STOP 'dom_mba'
663         ENDIF
664
665         ! Set to zero mbathy over islands if necessary  (lk_isl=F)
666         IF( .NOT. lk_isl ) THEN    ! No island
667            IF(lwp) WRITE(numout,*)
668            IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
669            IF(lwp) WRITE(numout,*) '         ----------------------------'
670
671            mbathy(:,:) = MAX( 0, mbathy(:,:) )
672
673            !  Boundary condition on mbathy
674            IF( .NOT.lk_mpp ) THEN 
675               !!bug ???  y reflechir!
676               !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
677               zbathy(:,:) = FLOAT( mbathy(:,:) )
678               CALL lbc_lnk( zbathy, 'T', 1. )
679               mbathy(:,:) = INT( zbathy(:,:) )
680            ENDIF
681
682         ENDIF
683
684      ENDIF
685
686      ! Number of ocean level inferior or equal to jpkm1
687
688      ikmax = 0
689      DO jj = 1, jpj
690         DO ji = 1, jpi
691            ikmax = MAX( ikmax, mbathy(ji,jj) )
692         END DO
693      END DO
694      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
695
696      IF( ikmax > jpkm1 ) THEN
697         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
698         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
699      ELSE IF( ikmax < jpkm1 ) THEN
700         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
701         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
702      ENDIF
703
704      IF( lwp .AND. nprint == 1 ) THEN
705         WRITE(numout,*)
706         WRITE(numout,*) ' bathymetric field '
707         WRITE(numout,*) ' ------------------'
708         WRITE(numout,*) ' number of non-zero T-levels '
709         CALL prihin( mbathy, jpi, jpj, 1, jpi,   &
710                      1     , 1  , jpj, 1, 3  ,   &
711                      numout )
712         WRITE(numout,*)
713      ENDIF
714
715   END SUBROUTINE zgr_bat_ctl
716
717
718#  include "domzgr_zps.h90"
719
720
721#  include "domzgr_s.h90"
722
723
724   !!======================================================================
725END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.