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

source: tags/nemo_dev_x5/NEMO/OPA_SRC/DOM/domzgr.F90 @ 1792

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