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

source: tags/nemo_v1_05/NEMO/OPA_SRC/DOM/domzgr.F90 @ 310

Last change on this file since 310 was 293, checked in by opalod, 19 years ago

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