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

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

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