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 tags/nemo_dev_x1/NEMO/OPA_SRC/DOM – NEMO

source: tags/nemo_dev_x1/NEMO/OPA_SRC/DOM/domzgr.F90 @ 46

Last change on this file since 46 was 46, checked in by cvs2svn, 20 years ago

This commit was manufactured by cvs2svn to create tag 'nemo_dev_x1'.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.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
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' 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) ::  &
269         clexp ,                     &  ! temporary characters
270         clname                         !    "           "
271      LOGICAL ::   llbon                ! check the existence of bathy files
272      INTEGER ::   ji, jj, jl, jn, jk   ! dummy loop indices
273      INTEGER ::   inum = 11            ! temporary logical unit
274      INTEGER ::   &
275         iim, ijm, ifreq,            &  ! temporary integers
276         il1, il2, ij, ii               ! (user to read "bathy_level" file
277      INTEGER  ::   &
278         ipi, ipj, ipk,              &  ! temporary integers
279         itime,                      &  !    "          "
280         ii_bump, ij_bump               ! bump center position
281      INTEGER, DIMENSION (1) ::   istep
282      INTEGER , DIMENSION(jpidta,jpjdta) ::   &
283         idta                           ! global domain integer data
284      REAL(wp) ::   &
285         r_bump, h_bump, h_oce,      &  ! bump characteristics
286         zi, zj, zdate0                 ! temporary scalars
287      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
288         zlamt, zphit,               &  ! temporary workspace (NetCDF read)
289         zdta                           ! global domain scalar data
290      REAL(wp), DIMENSION(jpk) ::   &
291         zdept                          ! temporary workspace (NetCDF read)
292      !!----------------------------------------------------------------------
293
294      IF(lwp) WRITE(numout,*)
295      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
296      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
297
298      ! ========================================
299      ! global domain level and meter bathymetry (idta,zdta)
300      ! ========================================
301      !                                               ! =============== !
302      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          ! defined by hand !
303         !                                            ! =============== !
304
305         IF( ntopo == 0 ) THEN                        ! flat basin
306
307            IF(lwp) WRITE(numout,*)
308            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
309
310            idta(:,:) = jpkm1                            ! flat basin
311            zdta(:,:) = gdepw(jpk)
312
313         ELSE                                         ! bump
314            IF(lwp) WRITE(numout,*)
315            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
316
317            ii_bump = jpidta / 3 + 3       ! i-index of the bump center
318            ij_bump = jpjdta / 2           ! j-index of the bump center
319            r_bump  =    6              ! bump radius (index)       
320            h_bump  =  240.e0           ! bump height (meters)
321            h_oce   = gdepw(jpk)        ! background ocean depth (meters)
322            IF(lwp) WRITE(numout,*) '            bump characteristics: '
323            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
324            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
325            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
326            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
327            ! zdta :
328            DO jj = 1, jpjdta
329               DO ji = 1, jpidta
330                  zi = FLOAT( ji - ii_bump ) / r_bump     
331                  zj = FLOAT( jj - ij_bump ) / r_bump       
332                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
333               END DO
334            END DO
335            ! idta :
336            idta(:,:) = jpkm1
337            DO jk = 1, jpkm1
338               DO jj = 1, jpjdta
339                  DO ji = 1, jpidta
340                     IF( gdept(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept(jk+1) )   idta(ji,jj) = jk
341                  END DO
342               END DO
343            END DO
344         ENDIF
345
346         ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
347         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
348            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0
349            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0.e0
350         ELSEIF( jperio == 2 ) THEN
351            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
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         ELSE
356            idta( :    , 1    ) = 0                 ;      zdta( :    , 1    ) =  0.e0
357            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0
358            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0
359            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0
360         ENDIF
361
362         !  EEL R5 configuration with east and west open boundaries.
363         !  Two rows of zeroes are needed at the south and north for OBCs
364         
365         IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
366            idta( : , 2      ) = 0                 ;      zdta( : , 2      ) =  0.e0
367!!CT            idta( : ,jpjdta-1) = 0                 ;      zdta( : ,jpjdta-1) =  0.e0
368         ENDIF
369
370         !                                            ! =============== !
371      ELSEIF( ntopo == 1 ) THEN                       !   read in file  !
372         !                                            ! =============== !
373
374         clname = 'bathy_level'                          ! level bathymetry
375         INQUIRE( FILE=clname, EXIST=llbon )
376         IF( llbon ) THEN
377            IF(lwp) WRITE(numout,*)
378            IF(lwp) WRITE(numout,*) '         read level bathymetry in ', clname
379            IF(lwp) WRITE(numout,*)
380!!bugs      OPEN( UNIT=inum, FILE=clname, FORM='FORMATTED', ACCESS='SEQUENTIAL', RECL=1 )
381            OPEN( UNIT=inum, FILE=clname, FORM='FORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD' )
382
383         REWIND(inum)
384         READ  (inum,9101) clexp, iim, ijm
385         READ  (inum,'(/)')
386         ifreq = 40
387         il1   = 1
388         DO jn = 1, jpidta/ifreq+1
389            READ(inum,'(/)')
390            il2 = MIN( jpidta, il1+ifreq-1 )
391            READ(inum,9201) ( ii, ji = il1, il2, 5 )
392            READ(inum,'(/)')
393            DO jj = jpjdta, 1, -1
394               READ(inum,9202) ij, ( idta(ji,jj), ji = il1, il2 )
395            END DO
396            il1 = il1 + ifreq
397         END DO
398            CLOSE(inum)
3999101     FORMAT(1x,a15,2i8)
4009201     FORMAT(3x,13(i3,12x))
4019202     FORMAT(i3,41i3)
402
403         ELSE
404            IF(lwp) WRITE(numout,cform_err)
405            IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
406            nstop = nstop + 1
407         ENDIF     
408
409         clname = 'bathy_meter.nc'                       ! meter bathymetry
410         INQUIRE( FILE=clname, EXIST=llbon )
411         IF( llbon ) THEN
412            IF(lwp) WRITE(numout,*)
413            IF(lwp) WRITE(numout,*) '         read meter bathymetry in ', clname
414            IF(lwp) WRITE(numout,*)
415            itime = 1
416            ipi = jpi
417            ipj = jpj
418            ipk = 1
419            CALL flinopen( clname, jpizoom, jpi, jpjzoom, jpj, .FALSE.,   &   
420                           ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, rdt, inum )
421            CALL flinget( inum, 'Bathymetry', jpidta, jpjdta, 1,   &
422                          itime, 1, 1, jpizoom, jpi, jpjzoom, jpj, zdta(:,:) ) 
423            CALL flinclo( inum )
424         ELSE
425            IF( lk_zps .OR. lk_sco ) THEN
426               IF(lwp) WRITE(numout,cform_err)       
427               IF(lwp) WRITE(numout,*)'    zgr_bat : unable to read the file', clname
428               nstop = nstop + 1
429            ELSE
430               zdta(:,:) = 0.e0
431               IF(lwp) WRITE(numout,*)'    zgr_bat : bathy_meter not found, but not used, bathy array set to zero'
432            ENDIF
433         ENDIF
434         !                                            ! =============== !
435      ELSE                                            !      error      !
436         !                                            ! =============== !
437         IF(lwp) WRITE(numout,cform_err)
438         IF(lwp) WRITE(numout,*) '          parameter , ntopo = ', ntopo
439         nstop = nstop + 1
440      ENDIF
441
442
443      ! =======================================
444      ! local domain level and meter bathymetry (mbathy,bathy)
445      ! =======================================
446
447      mbathy(:,:) = 0                                 ! set to zero extra halo points
448      bathy (:,:) = 0.e0                              ! (require for mpp case)
449
450      DO jj = 1, nlcj                                 ! interior values
451         DO ji = 1, nlci
452            mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
453            bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
454         END DO
455      END DO
456
457      ! =======================
458      ! NO closed seas or lakes
459      ! =======================
460
461      IF( .NOT. lclosea ) THEN
462         DO jl = 1, jpncs
463            DO jj = ncsj1(jl), ncsj2(jl)
464               DO ji = ncsi1(jl), ncsi2(jl)
465                  mbathy(ji,jj) = 0                   ! suppress closed seas
466                  bathy (ji,jj) = 0.e0                ! and lakes
467               END DO
468            END DO
469         END DO
470      ENDIF
471
472      ! ===========
473      ! Zoom domain
474      ! ===========
475
476      IF( lzoom )   CALL zgr_bat_zoom
477
478      ! ================
479      ! Bathymetry check
480      ! ================
481
482      CALL zgr_bat_ctl
483
484   END SUBROUTINE zgr_bat
485
486
487   SUBROUTINE zgr_bat_zoom
488      !!----------------------------------------------------------------------
489      !!                    ***  ROUTINE zgr_bat_zoom  ***
490      !!
491      !! ** Purpose : - Close zoom domain boundary if necessary
492      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
493      !!
494      !! ** Method  :
495      !!
496      !! ** Action  : - update mbathy: level bathymetry (in level index)
497      !!
498      !! History :
499      !!   9.0  !  03-08  (G. Madec)  Original code
500      !!----------------------------------------------------------------------
501      !! * Local variables
502      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
503      !!----------------------------------------------------------------------
504
505      IF(lwp) WRITE(numout,*)
506      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
507      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
508
509      ! Zoom domain
510      ! ===========
511
512      ! Forced closed boundary if required
513      IF( lzoom_e )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0
514      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
515      IF( lzoom_w )   mbathy( mi0(jpidta ):mi1(jpidta ) , :  ) = 0
516      IF( lzoom_n )   mbathy(  : , mj0(jpjdta ):mj1(jpjdta ) ) = 0
517
518      ! Configuration specific domain modifications
519      ! (here, ORCA arctic configuration: suppress Med Sea)
520      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
521         SELECT CASE ( jp_cfg )
522         !                                        ! =======================
523         CASE ( 2 )                               !  ORCA_R2 configuration
524            !                                     ! =======================
525            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
526
527            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
528            ij0 =  98   ;   ij1 = 110
529            !                                     ! =======================
530         CASE ( 05 )                              !  ORCA_R05 configuration
531            !                                     ! =======================
532            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
533 
534            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
535            ij0 = 314   ;   ij1 = 370 
536
537         END SELECT
538         !
539         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
540         !
541      ENDIF
542
543   END SUBROUTINE zgr_bat_zoom
544
545
546   SUBROUTINE zgr_bat_ctl
547      !!----------------------------------------------------------------------
548      !!                    ***  ROUTINE zgr_bat_ctl  ***
549      !!
550      !! ** Purpose :   check the bathymetry in levels
551      !!
552      !! ** Method  :   The array mbathy is checked to verified its consistency
553      !!      with the model options. in particular:
554      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
555      !!                  along closed boundary.
556      !!            mbathy must be cyclic IF jperio=1.
557      !!            mbathy must be lower or equal to jpk-1.
558      !!            isolated ocean grid points are suppressed from mbathy
559      !!                  since they are only connected to remaining
560      !!                  ocean through vertical diffusion.
561      !!      C A U T I O N : mbathy will be modified during the initializa-
562      !!      tion phase to become the number of non-zero w-levels of a water
563      !!      column, with a minimum value of 1.
564      !!
565      !! ** Action  : - update mbathy: level bathymetry (in level index)
566      !!              - update bathy : meter bathymetry (in meters)
567      !!
568      !! History :
569      !!   9.0  !  03-08  (G. Madec)  Original code
570      !!----------------------------------------------------------------------
571      !! * Local declarations
572      INTEGER ::   ji, jj, jl           ! dummy loop indices
573      INTEGER ::   &
574         icompt, ibtest, ikmax          ! temporary integers
575      REAL(wp), DIMENSION(jpi,jpj) ::   &
576         zbathy                         ! temporary workspace
577      !!----------------------------------------------------------------------
578
579      IF(lwp) WRITE(numout,*)
580      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
581      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
582
583      ! ================
584      ! Bathymetry check
585      ! ================
586
587      ! Suppress isolated ocean grid points
588
589      IF(lwp) WRITE(numout,*)
590      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
591      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
592
593      icompt = 0
594
595      DO jl = 1, 2
596
597         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
598            mbathy( 1 ,:) = mbathy(jpim1,:)
599            mbathy(jpi,:) = mbathy(  2  ,:)
600         ENDIF
601         DO jj = 2, jpjm1
602            DO ji = 2, jpim1
603               ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   &
604                  mbathy(ji,jj-1),mbathy(ji,jj+1) )
605               IF( ibtest < mbathy(ji,jj) ) THEN
606                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
607                     'grid-point (i,j) =  ',ji,jj,' is changed from ',   &
608                     mbathy(ji,jj),' to ', ibtest
609                  mbathy(ji,jj) = ibtest
610                  icompt = icompt + 1
611               ENDIF
612            END DO
613         END DO
614
615      END DO
616      IF( icompt == 0 ) THEN
617         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
618      ELSE
619         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
620      ENDIF
621      IF( lk_mpp ) THEN
622         zbathy(:,:) = FLOAT( mbathy(:,:) )
623         CALL lbc_lnk( zbathy, 'T', 1. )
624         mbathy(:,:) = INT( zbathy(:,:) )
625      ENDIF
626
627      ! 3.2 East-west cyclic boundary conditions
628
629      IF( nperio == 0 ) THEN
630         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   &
631            ' boundary: nperio = ', nperio
632         IF( lk_mpp ) THEN
633            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
634               IF( jperio /= 1 )   mbathy(1,:) = 0
635            ENDIF
636            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
637               IF( jperio /= 1 )   mbathy(nlci,:) = 0
638            ENDIF
639         ELSE
640            mbathy( 1 ,:) = 0
641            mbathy(jpi,:) = 0
642         ENDIF
643      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
644         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   &
645            ' on mbathy: nperio = ', nperio
646         mbathy( 1 ,:) = mbathy(jpim1,:)
647         mbathy(jpi,:) = mbathy(  2  ,:)
648      ELSEIF( nperio == 2 ) THEN
649         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   &
650            ' on mbathy: nperio = ', nperio
651      ELSE
652         IF(lwp) WRITE(numout,*) '    e r r o r'
653         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
654         !         STOP 'dom_mba'
655      ENDIF
656
657      ! Set to zero mbathy over islands if necessary  (l_isl=F)
658      IF( .NOT. l_isl ) THEN    ! No island
659         IF(lwp) WRITE(numout,*)
660         IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
661         IF(lwp) WRITE(numout,*) '         ----------------------------'
662
663         mbathy(:,:) = MAX( 0, mbathy(:,:) )
664
665         !  Boundary condition on mbathy
666         IF( .NOT.lk_mpp ) THEN 
667            !!bug ???  y reflechir!
668            !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
669            zbathy(:,:) = FLOAT( mbathy(:,:) )
670            CALL lbc_lnk( zbathy, 'T', 1. )
671            mbathy(:,:) = INT( zbathy(:,:) )
672         ENDIF
673
674      ENDIF
675
676      ! Number of ocean level inferior or equal to jpkm1
677
678      ikmax = 0
679      DO jj = 1, jpj
680         DO ji = 1, jpi
681            ikmax = MAX( ikmax, mbathy(ji,jj) )
682         END DO
683      END DO
684      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
685
686      IF( ikmax > jpkm1 ) THEN
687         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
688         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
689      ELSE IF( ikmax < jpkm1 ) THEN
690         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
691         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
692      ENDIF
693
694      IF( lwp .AND. nprint == 1 ) THEN
695         WRITE(numout,*)
696         WRITE(numout,*) ' bathymetric field '
697         WRITE(numout,*) ' ------------------'
698         WRITE(numout,*) ' number of non-zero T-levels '
699         CALL prihin( mbathy, jpi, jpj, 1, jpi,   &
700                      1     , 1  , jpj, 1, 3  ,   &
701                      numout )
702         WRITE(numout,*)
703      ENDIF
704
705   END SUBROUTINE zgr_bat_ctl
706
707
708#  include "domzgr_zps.h90"
709
710
711#  include "domzgr_s.h90"
712
713
714   !!======================================================================
715END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.