Ticket #1127: domzgr.F90

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