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

source: tags/nemo_dev_x2/NEMO/OPA_SRC/DOM/domzgr.F90 @ 630

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

CT : BUGFIX024 : # Remove the commentary '!CT' for EEL5 configuration

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