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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2482

Last change on this file since 2482 was 2482, checked in by gm, 14 years ago

v3.3beta: #777 backward compatibility with the setting of the minimum ocean depth: add a namelist param
  • Property svn:keywords set to Id
File size: 79.7 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dom_zgr          : defined the ocean vertical coordinate system
21   !!       zgr_bat      : bathymetry fields (levels and meters)
22   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
23   !!       zgr_bat_ctl  : check the bathymetry files
24   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
25   !!       zgr_z        : reference z-coordinate
26   !!       zgr_zco      : z-coordinate
27   !!       zgr_zps      : z-coordinate with partial steps
28   !!       zgr_sco      : s-coordinate
29   !!       fssig        : sigma coordinate non-dimensional function
30   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
31   !!---------------------------------------------------------------------
32   USE oce               ! ocean variables
33   USE dom_oce           ! ocean domain
34   USE closea            ! closed seas
35   USE c1d               ! 1D vertical configuration
36   USE in_out_manager    ! I/O manager
37   USE iom               ! I/O library
38   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
39   USE lib_mpp           ! distributed memory computing library
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   dom_zgr      ! called by dom_init.F90
45
46   !                                       !!* Namelist namzgr_sco *
47   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m)
48   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
49   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20)
50   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1)
51   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1)
52   LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
53   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching
54   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
55   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates
56 
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS       
66
67   SUBROUTINE dom_zgr
68      !!----------------------------------------------------------------------
69      !!                ***  ROUTINE dom_zgr  ***
70      !!                   
71      !! ** Purpose :  set the depth of model levels and the resulting
72      !!      vertical scale factors.
73      !!
74      !! ** Method  : - reference 1D vertical coordinate (gdep._0, e3._0)
75      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
76      !!              - vertical coordinate (gdep., e3.) depending on the
77      !!                coordinate chosen :
78      !!                   ln_zco=T   z-coordinate   
79      !!                   ln_zps=T   z-coordinate with partial steps
80      !!                   ln_zco=T   s-coordinate
81      !!
82      !! ** Action  :   define gdep., e3., mbathy and bathy
83      !!----------------------------------------------------------------------
84      INTEGER ::   ioptio = 0   ! temporary integer
85      !
86      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
87      !!----------------------------------------------------------------------
88
89      REWIND( numnam )                 ! Read Namelist namzgr : vertical coordinate'
90      READ  ( numnam, namzgr )
91
92      IF(lwp) THEN                     ! Control print
93         WRITE(numout,*)
94         WRITE(numout,*) 'dom_zgr : vertical coordinate'
95         WRITE(numout,*) '~~~~~~~'
96         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
97         WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
98         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
99         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco
100      ENDIF
101
102      ioptio = 0                       ! Check Vertical coordinate options
103      IF( ln_zco ) ioptio = ioptio + 1
104      IF( ln_zps ) ioptio = ioptio + 1
105      IF( ln_sco ) ioptio = ioptio + 1
106      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
107      !
108      ! Build the vertical coordinate system
109      ! ------------------------------------
110                          CALL zgr_z            ! Reference z-coordinate system (always called)
111                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
112      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
113      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
114      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
115      !
116      ! final adjustment of mbathy & check
117      ! -----------------------------------
118      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
119      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isoated ocean points
120                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
121      !
122      !
123      IF( nprint == 1 .AND. lwp )   THEN
124         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
125         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
126            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
127         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
128            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
129            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
130            &                   ' w ',   MINVAL( fse3w(:,:,:) )
131
132         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
133            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
134         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
135            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
136            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
137            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
138      ENDIF
139      !
140   END SUBROUTINE dom_zgr
141
142
143   SUBROUTINE zgr_z
144      !!----------------------------------------------------------------------
145      !!                   ***  ROUTINE zgr_z  ***
146      !!                   
147      !! ** Purpose :   set the depth of model levels and the resulting
148      !!      vertical scale factors.
149      !!
150      !! ** Method  :   z-coordinate system (use in all type of coordinate)
151      !!        The depth of model levels is defined from an analytical
152      !!      function the derivative of which gives the scale factors.
153      !!        both depth and scale factors only depend on k (1d arrays).
154      !!              w-level: gdepw_0  = fsdep(k)
155      !!                       e3w_0(k) = dk(fsdep)(k)     = fse3(k)
156      !!              t-level: gdept_0  = fsdep(k+0.5)
157      !!                       e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
158      !!
159      !! ** Action  : - gdept_0, gdepw_0 : depth of T- and W-point (m)
160      !!              - e3t_0  , e3w_0   : scale factors at T- and W-levels (m)
161      !!
162      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
163      !!----------------------------------------------------------------------
164      INTEGER  ::   jk                     ! dummy loop indices
165      REAL(wp) ::   zt, zw                 ! temporary scalars
166      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
167      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
168      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
169      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
170      !!----------------------------------------------------------------------
171
172      ! Set variables from parameters
173      ! ------------------------------
174       zkth = ppkth       ;   zacr = ppacr
175       zdzmin = ppdzmin   ;   zhmax = pphmax
176       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
177
178      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
179      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
180      IF(   ppa1  == pp_to_be_computed  .AND.  &
181         &  ppa0  == pp_to_be_computed  .AND.  &
182         &  ppsur == pp_to_be_computed           ) THEN
183         !
184         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
185            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
186            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
187         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
188         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
189      ELSE
190         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
191         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
192      ENDIF
193
194      IF(lwp) THEN                         ! Parameter print
195         WRITE(numout,*)
196         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
197         WRITE(numout,*) '    ~~~~~~~'
198         IF(  ppkth == 0._wp ) THEN             
199              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
200              WRITE(numout,*) '            Total depth    :', zhmax
201              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
202         ELSE
203            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
204               WRITE(numout,*) '         zsur, za0, za1 computed from '
205               WRITE(numout,*) '                 zdzmin = ', zdzmin
206               WRITE(numout,*) '                 zhmax  = ', zhmax
207            ENDIF
208            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
209            WRITE(numout,*) '                 zsur = ', zsur
210            WRITE(numout,*) '                 za0  = ', za0
211            WRITE(numout,*) '                 za1  = ', za1
212            WRITE(numout,*) '                 zkth = ', zkth
213            WRITE(numout,*) '                 zacr = ', zacr
214            IF( ldbletanh ) THEN
215               WRITE(numout,*) ' (Double tanh    za2  = ', za2
216               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
217               WRITE(numout,*) '                 zacr2= ', zacr2
218            ENDIF
219         ENDIF
220      ENDIF
221
222
223      ! Reference z-coordinate (depth - scale factor at T- and W-points)
224      ! ======================
225      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid       
226         za1 = zhmax / FLOAT(jpk-1) 
227         DO jk = 1, jpk
228            zw = FLOAT( jk )
229            zt = FLOAT( jk ) + 0.5_wp
230            gdepw_0(jk) = ( zw - 1 ) * za1
231            gdept_0(jk) = ( zt - 1 ) * za1
232            e3w_0  (jk) =  za1
233            e3t_0  (jk) =  za1
234         END DO
235      ELSE                                ! Madec & Imbard 1996 function
236         IF( .NOT. ldbletanh ) THEN
237            DO jk = 1, jpk
238               zw = REAL( jk , wp )
239               zt = REAL( jk , wp ) + 0.5_wp
240               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
241               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
242               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
243               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
244            END DO
245         ELSE
246            DO jk = 1, jpk
247               zw = FLOAT( jk )
248               zt = FLOAT( jk ) + 0.5_wp
249               ! Double tanh function
250               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
251                  &                            + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
252               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
253                  &                            + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
254               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )    &
255                  &                            + za2        * TANH(       (zw-zkth2) / zacr2 )
256               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )    &
257                  &                            + za2        * TANH(       (zt-zkth2) / zacr2 )
258            END DO
259         ENDIF
260         gdepw_0(1) = 0._wp                    ! force first w-level to be exactly at zero
261      ENDIF
262
263!!gm BUG in s-coordinate this does not work!
264      ! deepest/shallowest W level Above/Below ~10m
265      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 )                    ! ref. depth with tolerance (10% of minimum layer thickness)
266      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m
267      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
268!!gm end bug
269
270      IF(lwp) THEN                        ! control print
271         WRITE(numout,*)
272         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
273         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
274         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
275      ENDIF
276      DO jk = 1, jpk                      ! control positivity
277         IF( e3w_0  (jk) <= 0._wp .OR. e3t_0  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 '    )
278         IF( gdepw_0(jk) <  0._wp .OR. gdept_0(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' )
279      END DO
280      !
281   END SUBROUTINE zgr_z
282
283
284   SUBROUTINE zgr_bat
285      !!----------------------------------------------------------------------
286      !!                    ***  ROUTINE zgr_bat  ***
287      !!
288      !! ** Purpose :   set bathymetry both in levels and meters
289      !!
290      !! ** Method  :   read or define mbathy and bathy arrays
291      !!       * level bathymetry:
292      !!      The ocean basin geometry is given by a two-dimensional array,
293      !!      mbathy, which is defined as follow :
294      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
295      !!                              at t-point (ji,jj).
296      !!                            = 0  over the continental t-point.
297      !!      The array mbathy is checked to verified its consistency with
298      !!      model option. in particular:
299      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
300      !!                  along closed boundary.
301      !!            mbathy must be cyclic IF jperio=1.
302      !!            mbathy must be lower or equal to jpk-1.
303      !!            isolated ocean grid points are suppressed from mbathy
304      !!                  since they are only connected to remaining
305      !!                  ocean through vertical diffusion.
306      !!      ntopo=-1 :   rectangular channel or bassin with a bump
307      !!      ntopo= 0 :   flat rectangular channel or basin
308      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
309      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
310      !!
311      !! ** Action  : - mbathy: level bathymetry (in level index)
312      !!              - bathy : meter bathymetry (in meters)
313      !!----------------------------------------------------------------------
314      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
315      INTEGER  ::   inum                      ! temporary logical unit
316      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
317      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
318      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
319      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
320      INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data
321      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data
322      !!----------------------------------------------------------------------
323
324      IF(lwp) WRITE(numout,*)
325      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
326      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
327
328      !                                               ! ================== !
329      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
330         !                                            ! ================== !
331         !                                            ! global domain level and meter bathymetry (idta,zdta)
332         !
333         IF( ntopo == 0 ) THEN                        ! flat basin
334            IF(lwp) WRITE(numout,*)
335            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
336            idta(:,:) = jpkm1                            ! before last level
337            zdta(:,:) = gdepw_0(jpk)                     ! last w-point depth
338            h_oce     = gdepw_0(jpk)
339         ELSE                                         ! bump centered in the basin
340            IF(lwp) WRITE(numout,*)
341            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
342            ii_bump = jpidta / 2                           ! i-index of the bump center
343            ij_bump = jpjdta / 2                           ! j-index of the bump center
344            r_bump  = 50000._wp                            ! bump radius (meters)       
345            h_bump  =  2700._wp                            ! bump height (meters)
346            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
347            IF(lwp) WRITE(numout,*) '            bump characteristics: '
348            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
349            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
350            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
351            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
352            !                                       
353            DO jj = 1, jpjdta                              ! zdta :
354               DO ji = 1, jpidta
355                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
356                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
357                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
358               END DO
359            END DO
360            !                                              ! idta :
361            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
362               idta(:,:) = jpkm1
363            ELSE                                                ! z-coordinate (zco or zps): step-like topography
364               idta(:,:) = jpkm1
365               DO jk = 1, jpkm1
366                  WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) )   idta(:,:) = jk
367               END DO
368            ENDIF
369         ENDIF
370         !                                            ! set GLOBAL boundary conditions
371         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
372         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
373            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
374            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp
375         ELSEIF( jperio == 2 ) THEN
376            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
377            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp
378            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
379            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp
380         ELSE
381            ih = 0                                  ;      zh = 0._wp
382            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
383            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
384            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
385            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
386            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
387         ENDIF
388
389         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
390         mbathy(:,:) = 0                                   ! set to zero extra halo points
391         bathy (:,:) = 0._wp                               ! (require for mpp case)
392         DO jj = 1, nlcj                                   ! interior values
393            DO ji = 1, nlci
394               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
395               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
396            END DO
397         END DO
398         !
399         !                                            ! ================ !
400      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
401         !                                            ! ================ !
402         !
403         IF( ln_zco )   THEN                          ! zco : read level bathymetry
404            CALL iom_open ( 'bathy_level.nc', inum ) 
405            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
406            CALL iom_close( inum )
407            mbathy(:,:) = INT( bathy(:,:) )
408            !                                                ! =====================
409            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
410               !                                             ! =====================
411               IF( nn_cla == 0 ) THEN
412                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
413                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
414                  DO ji = mi0(ii0), mi1(ii1)
415                     DO jj = mj0(ij0), mj1(ij1)
416                        mbathy(ji,jj) = 15
417                     END DO
418                  END DO
419                  IF(lwp) WRITE(numout,*)
420                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
421                  !
422                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
423                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
424                  DO ji = mi0(ii0), mi1(ii1)
425                     DO jj = mj0(ij0), mj1(ij1)
426                        mbathy(ji,jj) = 12
427                     END DO
428                  END DO
429                  IF(lwp) WRITE(numout,*)
430                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
431               ENDIF
432               !
433            ENDIF
434            !
435         ENDIF
436         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
437            CALL iom_open ( 'bathy_meter.nc', inum ) 
438            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
439            CALL iom_close( inum )
440            !                                                ! =====================
441            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
442               ii0 = 142   ;   ii1 = 142                     ! =====================
443               ij0 =  51   ;   ij1 =  53                     
444               DO ji = mi0(ii0), mi1(ii1)                    ! Close Halmera Strait
445                  DO jj = mj0(ij0), mj1(ij1)
446                     bathy(ji,jj) = 0._wp
447                  END DO
448               END DO
449               IF(lwp) WRITE(numout,*)
450               IF(lwp) WRITE(numout,*) '      orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
451            ENDIF
452            !                                                ! =====================
453            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
454               !                                             ! =====================
455              IF( nn_cla == 0 ) THEN
456                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
457                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
458                 DO ji = mi0(ii0), mi1(ii1)
459                    DO jj = mj0(ij0), mj1(ij1)
460                       bathy(ji,jj) = 284._wp
461                    END DO
462                 END DO
463                 IF(lwp) WRITE(numout,*)
464                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
465                 !
466                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
467                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
468                 DO ji = mi0(ii0), mi1(ii1)
469                    DO jj = mj0(ij0), mj1(ij1)
470                       bathy(ji,jj) = 137._wp
471                    END DO
472                 END DO
473                 IF(lwp) WRITE(numout,*)
474                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
475              ENDIF
476              !
477           ENDIF
478            !
479        ENDIF
480         !                                            ! =============== !
481      ELSE                                            !      error      !
482         !                                            ! =============== !
483         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
484         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
485      ENDIF
486      !
487      !                                               ! =========================== !
488      IF( nclosea == 0 ) THEN                         !   NO closed seas or lakes   !
489         DO jl = 1, jpncs                             ! =========================== !
490            DO jj = ncsj1(jl), ncsj2(jl)
491               DO ji = ncsi1(jl), ncsi2(jl)
492                  mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry
493                  bathy (ji,jj) = 0._wp               
494               END DO
495            END DO
496         END DO
497      ENDIF
498      !
499      !                                               ! =========================== !
500      !                                               !     set a minimum depth     !
501      !                                               ! =========================== !
502      IF( rn_hmin < 0._wp ) THEN    ;   ik = INT( rn_hmin ) + 1                                    ! from a nb of level
503      ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
504      ENDIF
505      zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels
506      WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
507      ELSE WHERE                     ;   bathy(:,:) = MIN(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
508      END WHERE
509      IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
510      !
511   END SUBROUTINE zgr_bat
512
513
514   SUBROUTINE zgr_bat_zoom
515      !!----------------------------------------------------------------------
516      !!                    ***  ROUTINE zgr_bat_zoom  ***
517      !!
518      !! ** Purpose : - Close zoom domain boundary if necessary
519      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
520      !!
521      !! ** Method  :
522      !!
523      !! ** Action  : - update mbathy: level bathymetry (in level index)
524      !!----------------------------------------------------------------------
525      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
526      !!----------------------------------------------------------------------
527      !
528      IF(lwp) WRITE(numout,*)
529      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
530      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
531      !
532      ! Zoom domain
533      ! ===========
534      !
535      ! Forced closed boundary if required
536      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
537      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
538      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
539      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
540      !
541      ! Configuration specific domain modifications
542      ! (here, ORCA arctic configuration: suppress Med Sea)
543      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
544         SELECT CASE ( jp_cfg )
545         !                                        ! =======================
546         CASE ( 2 )                               !  ORCA_R2 configuration
547            !                                     ! =======================
548            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
549            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
550            ij0 =  98   ;   ij1 = 110
551            !                                     ! =======================
552         CASE ( 05 )                              !  ORCA_R05 configuration
553            !                                     ! =======================
554            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
555            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
556            ij0 = 314   ;   ij1 = 370 
557         END SELECT
558         !
559         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
560         !
561      ENDIF
562      !
563   END SUBROUTINE zgr_bat_zoom
564
565
566   SUBROUTINE zgr_bat_ctl
567      !!----------------------------------------------------------------------
568      !!                    ***  ROUTINE zgr_bat_ctl  ***
569      !!
570      !! ** Purpose :   check the bathymetry in levels
571      !!
572      !! ** Method  :   The array mbathy is checked to verified its consistency
573      !!      with the model options. in particular:
574      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
575      !!                  along closed boundary.
576      !!            mbathy must be cyclic IF jperio=1.
577      !!            mbathy must be lower or equal to jpk-1.
578      !!            isolated ocean grid points are suppressed from mbathy
579      !!                  since they are only connected to remaining
580      !!                  ocean through vertical diffusion.
581      !!      C A U T I O N : mbathy will be modified during the initializa-
582      !!      tion phase to become the number of non-zero w-levels of a water
583      !!      column, with a minimum value of 1.
584      !!
585      !! ** Action  : - update mbathy: level bathymetry (in level index)
586      !!              - update bathy : meter bathymetry (in meters)
587      !!----------------------------------------------------------------------
588      INTEGER ::   ji, jj, jl                    ! dummy loop indices
589      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
590      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace
591      !!----------------------------------------------------------------------
592
593      IF(lwp) WRITE(numout,*)
594      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
595      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
596
597      !                                          ! Suppress isolated ocean grid points
598      IF(lwp) WRITE(numout,*)
599      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
600      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
601      icompt = 0
602      DO jl = 1, 2
603         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
604            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
605            mbathy(jpi,:) = mbathy(  2  ,:)
606         ENDIF
607         DO jj = 2, jpjm1
608            DO ji = 2, jpim1
609               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
610                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
611               IF( ibtest < mbathy(ji,jj) ) THEN
612                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
613                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
614                  mbathy(ji,jj) = ibtest
615                  icompt = icompt + 1
616               ENDIF
617            END DO
618         END DO
619      END DO
620      IF( icompt == 0 ) THEN
621         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
622      ELSE
623         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
624      ENDIF
625      IF( lk_mpp ) THEN
626         zbathy(:,:) = FLOAT( mbathy(:,:) )
627         CALL lbc_lnk( zbathy, 'T', 1._wp )
628         mbathy(:,:) = INT( zbathy(:,:) )
629      ENDIF
630
631      !                                          ! East-west cyclic boundary conditions
632      IF( nperio == 0 ) THEN
633         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
634         IF( lk_mpp ) THEN
635            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
636               IF( jperio /= 1 )   mbathy(1,:) = 0
637            ENDIF
638            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
639               IF( jperio /= 1 )   mbathy(nlci,:) = 0
640            ENDIF
641         ELSE
642            IF( ln_zco .OR. ln_zps ) THEN
643               mbathy( 1 ,:) = 0
644               mbathy(jpi,:) = 0
645            ELSE
646               mbathy( 1 ,:) = jpkm1
647               mbathy(jpi,:) = jpkm1
648            ENDIF
649         ENDIF
650      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
651         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
652         mbathy( 1 ,:) = mbathy(jpim1,:)
653         mbathy(jpi,:) = mbathy(  2  ,:)
654      ELSEIF( nperio == 2 ) THEN
655         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
656      ELSE
657         IF(lwp) WRITE(numout,*) '    e r r o r'
658         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
659         !         STOP 'dom_mba'
660      ENDIF
661
662      !  Boundary condition on mbathy
663      IF( .NOT.lk_mpp ) THEN 
664!!gm     !!bug ???  think about it !
665         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
666         zbathy(:,:) = FLOAT( mbathy(:,:) )
667         CALL lbc_lnk( zbathy, 'T', 1._wp )
668         mbathy(:,:) = INT( zbathy(:,:) )
669      ENDIF
670
671      ! Number of ocean level inferior or equal to jpkm1
672      ikmax = 0
673      DO jj = 1, jpj
674         DO ji = 1, jpi
675            ikmax = MAX( ikmax, mbathy(ji,jj) )
676         END DO
677      END DO
678!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
679      IF( ikmax > jpkm1 ) THEN
680         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
681         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
682      ELSE IF( ikmax < jpkm1 ) THEN
683         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
684         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
685      ENDIF
686
687      IF( lwp .AND. nprint == 1 ) THEN      ! control print
688         WRITE(numout,*)
689         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
690         WRITE(numout,*) ' ------------------'
691         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
692         WRITE(numout,*)
693      ENDIF
694      !
695   END SUBROUTINE zgr_bat_ctl
696
697
698   SUBROUTINE zgr_bot_level
699      !!----------------------------------------------------------------------
700      !!                    ***  ROUTINE zgr_bot_level  ***
701      !!
702      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
703      !!
704      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
705      !!
706      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
707      !!                                     ocean level at t-, u- & v-points
708      !!                                     (min value = 1 over land)
709      !!----------------------------------------------------------------------
710      INTEGER ::   ji, jj   ! dummy loop indices
711      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
712      !!----------------------------------------------------------------------
713      !
714      IF(lwp) WRITE(numout,*)
715      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
716      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
717      !
718      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
719      !                                     ! bottom k-index of W-level = mbkt+1
720      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
721         DO ji = 1, jpim1
722            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
723            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
724         END DO
725      END DO
726      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
727      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
728      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
729      !
730   END SUBROUTINE zgr_bot_level
731
732
733   SUBROUTINE zgr_zco
734      !!----------------------------------------------------------------------
735      !!                  ***  ROUTINE zgr_zco  ***
736      !!
737      !! ** Purpose :   define the z-coordinate system
738      !!
739      !! ** Method  :   set 3D coord. arrays to reference 1D array
740      !!----------------------------------------------------------------------
741      INTEGER  ::   jk
742      !!----------------------------------------------------------------------
743      !
744      DO jk = 1, jpk
745         fsdept(:,:,jk) = gdept_0(jk)
746         fsdepw(:,:,jk) = gdepw_0(jk)
747         fsde3w(:,:,jk) = gdepw_0(jk)
748         fse3t (:,:,jk) = e3t_0(jk)
749         fse3u (:,:,jk) = e3t_0(jk)
750         fse3v (:,:,jk) = e3t_0(jk)
751         fse3f (:,:,jk) = e3t_0(jk)
752         fse3w (:,:,jk) = e3w_0(jk)
753         fse3uw(:,:,jk) = e3w_0(jk)
754         fse3vw(:,:,jk) = e3w_0(jk)
755      END DO
756      !
757   END SUBROUTINE zgr_zco
758
759
760   SUBROUTINE zgr_zps
761      !!----------------------------------------------------------------------
762      !!                  ***  ROUTINE zgr_zps  ***
763      !!                     
764      !! ** Purpose :   the depth and vertical scale factor in partial step
765      !!      z-coordinate case
766      !!
767      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
768      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
769      !!      a partial step representation of bottom topography.
770      !!
771      !!        The reference depth of model levels is defined from an analytical
772      !!      function the derivative of which gives the reference vertical
773      !!      scale factors.
774      !!        From  depth and scale factors reference, we compute there new value
775      !!      with partial steps  on 3d arrays ( i, j, k ).
776      !!
777      !!              w-level: gdepw(i,j,k)  = fsdep(k)
778      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
779      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
780      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
781      !!
782      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
783      !!      we find the mbathy index of the depth at each grid point.
784      !!      This leads us to three cases:
785      !!
786      !!              - bathy = 0 => mbathy = 0
787      !!              - 1 < mbathy < jpkm1   
788      !!              - bathy > gdepw(jpk) => mbathy = jpkm1 
789      !!
790      !!        Then, for each case, we find the new depth at t- and w- levels
791      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
792      !!      and f-points.
793      !!
794      !!        This routine is given as an example, it must be modified
795      !!      following the user s desiderata. nevertheless, the output as
796      !!      well as the way to compute the model levels and scale factors
797      !!      must be respected in order to insure second order accuracy
798      !!      schemes.
799      !!
800      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
801      !!         - - - - - - -   gdept, gdepw and e3. are positives
802      !!     
803      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
804      !!----------------------------------------------------------------------
805      INTEGER  ::   ji, jj, jk       ! dummy loop indices
806      INTEGER  ::   ik, it           ! temporary integers
807      INTEGER  ::   nlbelow          ! temporary integer
808      LOGICAL  ::   ll_print         ! Allow  control print for debugging
809      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
810      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
811      REAL(wp) ::   zmax             ! Maximum depth
812      REAL(wp) ::   zdiff            ! temporary scalar
813      REAL(wp) ::   zrefdep          ! temporary scalar
814      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace
815      !!---------------------------------------------------------------------
816
817      IF(lwp) WRITE(numout,*)
818      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
819      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
820      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
821
822      ll_print = .FALSE.                   ! Local variable for debugging
823!!    ll_print = .TRUE.
824     
825      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
826         WRITE(numout,*)
827         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
828         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
829      ENDIF
830
831
832      ! bathymetry in level (from bathy_meter)
833      ! ===================
834      zmax = gdepw_0(jpk) + e3t_0(jpk)          ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
835      mbathy(:,:) = jpkm1                       ! initialize mbathy to the maximum ocean level available
836      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
837
838      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
839      ! find the number of ocean levels such that the last level thickness
840      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
841      ! e3t_0 is the reference level thickness
842      DO jk = jpkm1, 1, -1
843         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
844         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
845      END DO
846
847      ! Scale factors and depth at T- and W-points
848      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
849         gdept(:,:,jk) = gdept_0(jk)
850         gdepw(:,:,jk) = gdepw_0(jk)
851         e3t  (:,:,jk) = e3t_0  (jk)
852         e3w  (:,:,jk) = e3w_0  (jk)
853      END DO
854      !
855      DO jj = 1, jpj
856         DO ji = 1, jpi
857            ik = mbathy(ji,jj)
858            IF( ik > 0 ) THEN               ! ocean point only
859               ! max ocean level case
860               IF( ik == jpkm1 ) THEN
861                  zdepwp = bathy(ji,jj)
862                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
863                  ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
864                  e3t(ji,jj,ik  ) = ze3tp
865                  e3t(ji,jj,ik+1) = ze3tp
866                  e3w(ji,jj,ik  ) = ze3wp
867                  e3w(ji,jj,ik+1) = ze3tp
868                  gdepw(ji,jj,ik+1) = zdepwp
869                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
870                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
871                  !
872               ELSE                         ! standard case
873                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
874                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
875                  ENDIF
876!gm Bug?  check the gdepw_0
877                  !       ... on ik
878                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
879                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
880                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
881                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
882                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
883                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) )   &
884                     &                     * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
885                  !       ... on ik+1
886                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
887                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
888                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
889               ENDIF
890            ENDIF
891         END DO
892      END DO
893      !
894      it = 0
895      DO jj = 1, jpj
896         DO ji = 1, jpi
897            ik = mbathy(ji,jj)
898            IF( ik > 0 ) THEN               ! ocean point only
899               e3tp (ji,jj) = e3t(ji,jj,ik  )
900               e3wp (ji,jj) = e3w(ji,jj,ik  )
901               ! test
902               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
903               IF( zdiff <= 0._wp .AND. lwp ) THEN
904                  it = it + 1
905                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
906                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
907                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
908                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
909               ENDIF
910            ENDIF
911         END DO
912      END DO
913
914      ! Scale factors and depth at U-, V-, UW and VW-points
915      DO jk = 1, jpk                        ! initialisation to z-scale factors
916         e3u (:,:,jk) = e3t_0(jk)
917         e3v (:,:,jk) = e3t_0(jk)
918         e3uw(:,:,jk) = e3w_0(jk)
919         e3vw(:,:,jk) = e3w_0(jk)
920      END DO
921      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
922         DO jj = 1, jpjm1
923            DO ji = 1, fs_jpim1   ! vector opt.
924               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
925               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
926               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
927               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
928            END DO
929         END DO
930      END DO
931      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
932      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
933      !
934      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
935         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_0(jk)
936         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_0(jk)
937         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_0(jk)
938         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_0(jk)
939      END DO
940     
941      ! Scale factor at F-point
942      DO jk = 1, jpk                        ! initialisation to z-scale factors
943         e3f(:,:,jk) = e3t_0(jk)
944      END DO
945      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
946         DO jj = 1, jpjm1
947            DO ji = 1, fs_jpim1   ! vector opt.
948               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
949            END DO
950         END DO
951      END DO
952      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
953      !
954      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
955         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_0(jk)
956      END DO
957!!gm  bug ? :  must be a do loop with mj0,mj1
958      !
959      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
960      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
961      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
962      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
963      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
964
965      ! Control of the sign
966      IF( MINVAL( e3t  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
967      IF( MINVAL( e3w  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
968      IF( MINVAL( gdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
969      IF( MINVAL( gdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
970     
971      ! Compute gdep3w (vertical sum of e3w)
972      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
973      DO jk = 2, jpk
974         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
975      END DO
976       
977      !                                               ! ================= !
978      IF(lwp .AND. ll_print) THEN                     !   Control print   !
979         !                                            ! ================= !
980         DO jj = 1,jpj
981            DO ji = 1, jpi
982               ik = MAX( mbathy(ji,jj), 1 )
983               zprt(ji,jj,1) = e3t   (ji,jj,ik)
984               zprt(ji,jj,2) = e3w   (ji,jj,ik)
985               zprt(ji,jj,3) = e3u   (ji,jj,ik)
986               zprt(ji,jj,4) = e3v   (ji,jj,ik)
987               zprt(ji,jj,5) = e3f   (ji,jj,ik)
988               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
989            END DO
990         END DO
991         WRITE(numout,*)
992         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
993         WRITE(numout,*)
994         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
995         WRITE(numout,*)
996         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
997         WRITE(numout,*)
998         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
999         WRITE(numout,*)
1000         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1001         WRITE(numout,*)
1002         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1003      ENDIF 
1004      !
1005   END SUBROUTINE zgr_zps
1006
1007
1008   FUNCTION fssig( pk ) RESULT( pf )
1009      !!----------------------------------------------------------------------
1010      !!                 ***  ROUTINE eos_init  ***
1011      !!       
1012      !! ** Purpose :   provide the analytical function in s-coordinate
1013      !!         
1014      !! ** Method  :   the function provide the non-dimensional position of
1015      !!                T and W (i.e. between 0 and 1)
1016      !!                T-points at integer values (between 1 and jpk)
1017      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1018      !!----------------------------------------------------------------------
1019      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1020      REAL(wp)             ::   pf   ! sigma value
1021      !!----------------------------------------------------------------------
1022      !
1023      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1024         &     - TANH( rn_thetb * rn_theta                                )  )   &
1025         & * (   COSH( rn_theta                           )                      &
1026         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1027         & / ( 2._wp * SINH( rn_theta ) )
1028      !
1029   END FUNCTION fssig
1030
1031
1032   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1033      !!----------------------------------------------------------------------
1034      !!                 ***  ROUTINE eos_init  ***
1035      !!
1036      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1037      !!
1038      !! ** Method  :   the function provides the non-dimensional position of
1039      !!                T and W (i.e. between 0 and 1)
1040      !!                T-points at integer values (between 1 and jpk)
1041      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1042      !!----------------------------------------------------------------------
1043      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1044      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1045      REAL(wp)             ::   pf1   ! sigma value
1046      !!----------------------------------------------------------------------
1047      !
1048      IF ( rn_theta == 0 ) then      ! uniform sigma
1049         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1050      ELSE                        ! stretched sigma
1051         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1052            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1053            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1054      ENDIF
1055      !
1056   END FUNCTION fssig1
1057
1058
1059   SUBROUTINE zgr_sco
1060      !!----------------------------------------------------------------------
1061      !!                  ***  ROUTINE zgr_sco  ***
1062      !!                     
1063      !! ** Purpose :   define the s-coordinate system
1064      !!
1065      !! ** Method  :   s-coordinate
1066      !!         The depth of model levels is defined as the product of an
1067      !!      analytical function by the local bathymetry, while the vertical
1068      !!      scale factors are defined as the product of the first derivative
1069      !!      of the analytical function by the bathymetry.
1070      !!      (this solution save memory as depth and scale factors are not
1071      !!      3d fields)
1072      !!          - Read bathymetry (in meters) at t-point and compute the
1073      !!         bathymetry at u-, v-, and f-points.
1074      !!            hbatu = mi( hbatt )
1075      !!            hbatv = mj( hbatt )
1076      !!            hbatf = mi( mj( hbatt ) )
1077      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
1078      !!         function and its derivative given as function.
1079      !!            gsigt(k) = fssig (k    )
1080      !!            gsigw(k) = fssig (k-0.5)
1081      !!            esigt(k) = fsdsig(k    )
1082      !!            esigw(k) = fsdsig(k-0.5)
1083      !!      This routine is given as an example, it must be modified
1084      !!      following the user s desiderata. nevertheless, the output as
1085      !!      well as the way to compute the model levels and scale factors
1086      !!      must be respected in order to insure second order a!!uracy
1087      !!      schemes.
1088      !!
1089      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1090      !!----------------------------------------------------------------------
1091      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1092      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1093      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars
1094      REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace
1095      REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
1096      !!
1097      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigw3
1098      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigt3
1099      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsi3w3
1100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigt3
1101      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigw3
1102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtu3
1103      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtv3
1104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtf3
1105      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwu3
1106      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwv3
1107      !!
1108      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1109      !!----------------------------------------------------------------------
1110
1111      REWIND( numnam )                        ! Read Namelist namzgr_sco : sigma-stretching parameters
1112      READ  ( numnam, namzgr_sco )
1113
1114      IF(lwp) THEN                            ! control print
1115         WRITE(numout,*)
1116         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1117         WRITE(numout,*) '~~~~~~~~~~~'
1118         WRITE(numout,*) '   Namelist namzgr_sco'
1119         WRITE(numout,*) '      sigma-stretching coeffs '
1120         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
1121         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
1122         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
1123         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
1124         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
1125         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
1126         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
1127         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
1128      ENDIF
1129
1130      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
1131      esigt3  = 0._wp   ;   esigw3  = 0._wp 
1132      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
1133      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
1134
1135      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1136      hifu(:,:) = rn_sbot_min
1137      hifv(:,:) = rn_sbot_min
1138      hiff(:,:) = rn_sbot_min
1139
1140      !                                        ! set maximum ocean depth
1141      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1142
1143      DO jj = 1, jpj
1144         DO ji = 1, jpi
1145           IF( bathy(ji,jj) > 0._wp ) THEN
1146              bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1147           ENDIF
1148         END DO
1149      END DO
1150      !                                        ! =============================
1151      !                                        ! Define the envelop bathymetry   (hbatt)
1152      !                                        ! =============================
1153      ! use r-value to create hybrid coordinates
1154      DO jj = 1, jpj
1155         DO ji = 1, jpi
1156            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1157         END DO
1158      END DO
1159      !
1160      ! Smooth the bathymetry (if required)
1161      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1162      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1163      !
1164      jl = 0
1165      zrmax = 1._wp
1166      !                                                     ! ================ !
1167      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
1168         !                                                  ! ================ !
1169         jl = jl + 1
1170         zrmax = 0._wp
1171         zmsk(:,:) = 0._wp
1172         DO jj = 1, nlcj
1173            DO ji = 1, nlci
1174               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1175               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1176               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1177               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1178               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1179               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1180               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
1181               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1182               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
1183            END DO
1184         END DO
1185         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1186         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1187         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
1188         DO jj = 1, nlcj
1189            DO ji = 1, nlci
1190                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1191            END DO
1192         END DO
1193         !
1194         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1195         !
1196         DO jj = 1, nlcj
1197            DO ji = 1, nlci
1198               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
1199               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
1200               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
1201               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
1202               IF( zmsk(ji,jj) == 1._wp ) THEN
1203                  ztmp(ji,jj) =   (                                                                                   &
1204             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
1205             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
1206             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
1207             &                    ) / (                                                                               &
1208             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
1209             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
1210             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
1211             &                        )
1212               ENDIF
1213            END DO
1214         END DO
1215         !
1216         DO jj = 1, nlcj
1217            DO ji = 1, nlci
1218               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1219            END DO
1220         END DO
1221         !
1222         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
1223         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
1224         DO jj = 1, nlcj
1225            DO ji = 1, nlci
1226               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
1227            END DO
1228         END DO
1229         !                                                  ! ================ !
1230      END DO                                                !     End loop     !
1231      !                                                     ! ================ !
1232      !
1233      !                                        ! envelop bathymetry saved in hbatt
1234      hbatt(:,:) = zenv(:,:) 
1235      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1236         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1237         DO jj = 1, jpj
1238            DO ji = 1, jpi
1239               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
1240               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1241            END DO
1242         END DO
1243      ENDIF
1244      !
1245      IF(lwp) THEN                             ! Control print
1246         WRITE(numout,*)
1247         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1248         WRITE(numout,*)
1249         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1250         IF( nprint == 1 )   THEN       
1251            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1252            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1253         ENDIF
1254      ENDIF
1255
1256      !                                        ! ==============================
1257      !                                        !   hbatu, hbatv, hbatf fields
1258      !                                        ! ==============================
1259      IF(lwp) THEN
1260         WRITE(numout,*)
1261         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1262      ENDIF
1263      hbatu(:,:) = rn_sbot_min
1264      hbatv(:,:) = rn_sbot_min
1265      hbatf(:,:) = rn_sbot_min
1266      DO jj = 1, jpjm1
1267        DO ji = 1, jpim1   ! NO vector opt.
1268           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1269           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1270           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1271              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1272        END DO
1273      END DO
1274      !
1275      ! Apply lateral boundary condition
1276!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1277      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1278      DO jj = 1, jpj
1279         DO ji = 1, jpi
1280            IF( hbatu(ji,jj) == 0._wp ) THEN
1281               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1282               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1283            ENDIF
1284         END DO
1285      END DO
1286      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1287      DO jj = 1, jpj
1288         DO ji = 1, jpi
1289            IF( hbatv(ji,jj) == 0._wp ) THEN
1290               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1291               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1292            ENDIF
1293         END DO
1294      END DO
1295      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1296      DO jj = 1, jpj
1297         DO ji = 1, jpi
1298            IF( hbatf(ji,jj) == 0._wp ) THEN
1299               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1300               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1301            ENDIF
1302         END DO
1303      END DO
1304
1305!!bug:  key_helsinki a verifer
1306      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1307      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1308      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1309      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1310
1311      IF( nprint == 1 .AND. lwp )   THEN
1312         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1313            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1314         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1315            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1316         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1317            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1318         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1319            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1320      ENDIF
1321!! helsinki
1322
1323      !                                            ! =======================
1324      !                                            !   s-ccordinate fields     (gdep., e3.)
1325      !                                            ! =======================
1326      !
1327      ! non-dimensional "sigma" for model level depth at w- and t-levels
1328
1329      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
1330         !                         ! below rn_hc, with uniform sigma in shallower waters
1331         DO ji = 1, jpi
1332            DO jj = 1, jpj
1333
1334               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1335                  DO jk = 1, jpk
1336                     gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1337                     gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1338                  END DO
1339               ELSE ! shallow water, uniform sigma
1340                  DO jk = 1, jpk
1341                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1342                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1343                  END DO
1344               ENDIF
1345               IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1346               !
1347               DO jk = 1, jpkm1
1348                  esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1349                  esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1350               END DO
1351               esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
1352               esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
1353               !
1354               ! Coefficients for vertical depth as the sum of e3w scale factors
1355               gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
1356               DO jk = 2, jpk
1357                  gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1358               END DO
1359               !
1360               DO jk = 1, jpk
1361                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1362                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1363                  gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1364                  gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1365                  gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
1366               END DO
1367               !
1368            END DO   ! for all jj's
1369         END DO    ! for all ji's
1370
1371         DO ji = 1, jpi
1372            DO jj = 1, jpj
1373               DO jk = 1, jpk
1374                  esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
1375                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1376                  esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
1377                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1378                  esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
1379                     &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
1380                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1381                  esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
1382                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1383                  esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
1384                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1385                  !
1386                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1387                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1388                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1389                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1390                  !
1391                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1392                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1393                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1394               END DO
1395            END DO
1396         END DO
1397         !
1398      ELSE   ! not ln_s_sigma
1399         !
1400         DO jk = 1, jpk
1401           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1402           gsigt(jk) = -fssig( REAL(jk,wp)        )
1403         END DO
1404         IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1405         !
1406         ! Coefficients for vertical scale factors at w-, t- levels
1407!!gm bug :  define it from analytical function, not like juste bellow....
1408!!gm        or betteroffer the 2 possibilities....
1409         DO jk = 1, jpkm1
1410            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
1411            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1412         END DO
1413         esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
1414         esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
1415
1416!!gm  original form
1417!!org DO jk = 1, jpk
1418!!org    esigt(jk)=fsdsig( FLOAT(jk)     )
1419!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1420!!org END DO
1421!!gm
1422         !
1423         ! Coefficients for vertical depth as the sum of e3w scale factors
1424         gsi3w(1) = 0.5_wp * esigw(1)
1425         DO jk = 2, jpk
1426            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1427         END DO
1428!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1429         DO jk = 1, jpk
1430            zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1431            zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1432            gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
1433            gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
1434            gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
1435         END DO
1436!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1437         DO jj = 1, jpj
1438            DO ji = 1, jpi
1439               DO jk = 1, jpk
1440                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1441                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1442                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1443                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1444                 !
1445                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1446                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1447                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1448               END DO
1449            END DO
1450         END DO
1451         !
1452      ENDIF ! ln_s_sigma
1453
1454
1455      !
1456!!    H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
1457
1458      fsdept(:,:,:) = gdept (:,:,:)
1459      fsdepw(:,:,:) = gdepw (:,:,:)
1460      fsde3w(:,:,:) = gdep3w(:,:,:)
1461      fse3t (:,:,:) = e3t   (:,:,:)
1462      fse3u (:,:,:) = e3u   (:,:,:)
1463      fse3v (:,:,:) = e3v   (:,:,:)
1464      fse3f (:,:,:) = e3f   (:,:,:)
1465      fse3w (:,:,:) = e3w   (:,:,:)
1466      fse3uw(:,:,:) = e3uw  (:,:,:)
1467      fse3vw(:,:,:) = e3vw  (:,:,:)
1468!!
1469      ! HYBRID :
1470      DO jj = 1, jpj
1471         DO ji = 1, jpi
1472            DO jk = 1, jpkm1
1473               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1474               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
1475            END DO
1476         END DO
1477      END DO
1478      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1479         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1480
1481      !                                               ! =============
1482      IF(lwp) THEN                                    ! Control print
1483         !                                            ! =============
1484         WRITE(numout,*) 
1485         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1486         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
1487         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1488      ENDIF
1489      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1490         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
1491         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
1492            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
1493         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
1494            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
1495            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
1496            &                          ' w ', MINVAL( fse3w (:,:,:) )
1497
1498         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
1499            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
1500         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
1501            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
1502            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
1503            &                          ' w ', MAXVAL( fse3w (:,:,:) )
1504      ENDIF
1505      !
1506      IF(lwp) THEN                                  ! selected vertical profiles
1507         WRITE(numout,*)
1508         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1509         WRITE(numout,*) ' ~~~~~~  --------------------'
1510         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1511         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
1512            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1513         DO jj = mj0(20), mj1(20)
1514            DO ji = mi0(20), mi1(20)
1515               WRITE(numout,*)
1516               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1517               WRITE(numout,*) ' ~~~~~~  --------------------'
1518               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1519               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1520                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1521            END DO
1522         END DO
1523         DO jj = mj0(74), mj1(74)
1524            DO ji = mi0(100), mi1(100)
1525               WRITE(numout,*)
1526               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1527               WRITE(numout,*) ' ~~~~~~  --------------------'
1528               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1529               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1530                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1531            END DO
1532         END DO
1533      ENDIF
1534
1535!!gm bug?  no more necessary?  if ! defined key_helsinki
1536      DO jk = 1, jpk
1537         DO jj = 1, jpj
1538            DO ji = 1, jpi
1539               IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1540                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1541                  CALL ctl_stop( ctmp1 )
1542               ENDIF
1543               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1544                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1545                  CALL ctl_stop( ctmp1 )
1546               ENDIF
1547            END DO
1548         END DO
1549      END DO
1550!!gm bug    #endif
1551      !
1552   END SUBROUTINE zgr_sco
1553
1554
1555   !!======================================================================
1556END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.