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

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

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

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