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

Last change on this file since 136 was 136, checked in by opalod, 20 years ago

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