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 on Ticket #936 – Attachment – NEMO

Ticket #936: domzgr.F90

File domzgr.F90, 80.6 KB (added by gm, 12 years ago)

use of clo_bat in domzgr

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