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 @ 3

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

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