source: trunk/NEMO/OPA_SRC/DOM/domzgr.F90 @ 392

Last change on this file since 392 was 392, checked in by opalod, 15 years ago

RB:nemo_v1_update_038: first integration of Agrif :

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