New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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