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

Last change on this file since 2712 was 2712, checked in by rfurner, 13 years ago

bug fix for s coordinates and use of rn_hmin
  • Property svn:keywords set to Id
File size: 79.8 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 ( .not. ln_sco ) THEN
503         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
504         ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
505         ENDIF
506         zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels
507         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
508         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
509         END WHERE
510         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
511      ENDIF
512      !
513   END SUBROUTINE zgr_bat
514
515
516   SUBROUTINE zgr_bat_zoom
517      !!----------------------------------------------------------------------
518      !!                    ***  ROUTINE zgr_bat_zoom  ***
519      !!
520      !! ** Purpose : - Close zoom domain boundary if necessary
521      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
522      !!
523      !! ** Method  :
524      !!
525      !! ** Action  : - update mbathy: level bathymetry (in level index)
526      !!----------------------------------------------------------------------
527      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
528      !!----------------------------------------------------------------------
529      !
530      IF(lwp) WRITE(numout,*)
531      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
532      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
533      !
534      ! Zoom domain
535      ! ===========
536      !
537      ! Forced closed boundary if required
538      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
539      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
540      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
541      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
542      !
543      ! Configuration specific domain modifications
544      ! (here, ORCA arctic configuration: suppress Med Sea)
545      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
546         SELECT CASE ( jp_cfg )
547         !                                        ! =======================
548         CASE ( 2 )                               !  ORCA_R2 configuration
549            !                                     ! =======================
550            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
551            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
552            ij0 =  98   ;   ij1 = 110
553            !                                     ! =======================
554         CASE ( 05 )                              !  ORCA_R05 configuration
555            !                                     ! =======================
556            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
557            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
558            ij0 = 314   ;   ij1 = 370 
559         END SELECT
560         !
561         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
562         !
563      ENDIF
564      !
565   END SUBROUTINE zgr_bat_zoom
566
567
568   SUBROUTINE zgr_bat_ctl
569      !!----------------------------------------------------------------------
570      !!                    ***  ROUTINE zgr_bat_ctl  ***
571      !!
572      !! ** Purpose :   check the bathymetry in levels
573      !!
574      !! ** Method  :   The array mbathy is checked to verified its consistency
575      !!      with the model options. in particular:
576      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
577      !!                  along closed boundary.
578      !!            mbathy must be cyclic IF jperio=1.
579      !!            mbathy must be lower or equal to jpk-1.
580      !!            isolated ocean grid points are suppressed from mbathy
581      !!                  since they are only connected to remaining
582      !!                  ocean through vertical diffusion.
583      !!      C A U T I O N : mbathy will be modified during the initializa-
584      !!      tion phase to become the number of non-zero w-levels of a water
585      !!      column, with a minimum value of 1.
586      !!
587      !! ** Action  : - update mbathy: level bathymetry (in level index)
588      !!              - update bathy : meter bathymetry (in meters)
589      !!----------------------------------------------------------------------
590      INTEGER ::   ji, jj, jl                    ! dummy loop indices
591      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
592      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace
593      !!----------------------------------------------------------------------
594
595      IF(lwp) WRITE(numout,*)
596      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
597      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
598
599      !                                          ! Suppress isolated ocean grid points
600      IF(lwp) WRITE(numout,*)
601      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
602      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
603      icompt = 0
604      DO jl = 1, 2
605         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
606            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
607            mbathy(jpi,:) = mbathy(  2  ,:)
608         ENDIF
609         DO jj = 2, jpjm1
610            DO ji = 2, jpim1
611               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
612                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
613               IF( ibtest < mbathy(ji,jj) ) THEN
614                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
615                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
616                  mbathy(ji,jj) = ibtest
617                  icompt = icompt + 1
618               ENDIF
619            END DO
620         END DO
621      END DO
622      IF( icompt == 0 ) THEN
623         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
624      ELSE
625         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
626      ENDIF
627      IF( lk_mpp ) THEN
628         zbathy(:,:) = FLOAT( mbathy(:,:) )
629         CALL lbc_lnk( zbathy, 'T', 1._wp )
630         mbathy(:,:) = INT( zbathy(:,:) )
631      ENDIF
632
633      !                                          ! East-west cyclic boundary conditions
634      IF( nperio == 0 ) THEN
635         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
636         IF( lk_mpp ) THEN
637            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
638               IF( jperio /= 1 )   mbathy(1,:) = 0
639            ENDIF
640            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
641               IF( jperio /= 1 )   mbathy(nlci,:) = 0
642            ENDIF
643         ELSE
644            IF( ln_zco .OR. ln_zps ) THEN
645               mbathy( 1 ,:) = 0
646               mbathy(jpi,:) = 0
647            ELSE
648               mbathy( 1 ,:) = jpkm1
649               mbathy(jpi,:) = jpkm1
650            ENDIF
651         ENDIF
652      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
653         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
654         mbathy( 1 ,:) = mbathy(jpim1,:)
655         mbathy(jpi,:) = mbathy(  2  ,:)
656      ELSEIF( nperio == 2 ) THEN
657         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
658      ELSE
659         IF(lwp) WRITE(numout,*) '    e r r o r'
660         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
661         !         STOP 'dom_mba'
662      ENDIF
663
664      !  Boundary condition on mbathy
665      IF( .NOT.lk_mpp ) THEN 
666!!gm     !!bug ???  think about it !
667         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
668         zbathy(:,:) = FLOAT( mbathy(:,:) )
669         CALL lbc_lnk( zbathy, 'T', 1._wp )
670         mbathy(:,:) = INT( zbathy(:,:) )
671      ENDIF
672
673      ! Number of ocean level inferior or equal to jpkm1
674      ikmax = 0
675      DO jj = 1, jpj
676         DO ji = 1, jpi
677            ikmax = MAX( ikmax, mbathy(ji,jj) )
678         END DO
679      END DO
680!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
681      IF( ikmax > jpkm1 ) THEN
682         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
683         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
684      ELSE IF( ikmax < jpkm1 ) THEN
685         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
686         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
687      ENDIF
688
689      IF( lwp .AND. nprint == 1 ) THEN      ! control print
690         WRITE(numout,*)
691         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
692         WRITE(numout,*) ' ------------------'
693         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
694         WRITE(numout,*)
695      ENDIF
696      !
697   END SUBROUTINE zgr_bat_ctl
698
699
700   SUBROUTINE zgr_bot_level
701      !!----------------------------------------------------------------------
702      !!                    ***  ROUTINE zgr_bot_level  ***
703      !!
704      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
705      !!
706      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
707      !!
708      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
709      !!                                     ocean level at t-, u- & v-points
710      !!                                     (min value = 1 over land)
711      !!----------------------------------------------------------------------
712      INTEGER ::   ji, jj   ! dummy loop indices
713      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
714      !!----------------------------------------------------------------------
715      !
716      IF(lwp) WRITE(numout,*)
717      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
718      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
719      !
720      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
721      !                                     ! bottom k-index of W-level = mbkt+1
722      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
723         DO ji = 1, jpim1
724            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
725            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
726         END DO
727      END DO
728      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
729      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
730      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
731      !
732   END SUBROUTINE zgr_bot_level
733
734
735   SUBROUTINE zgr_zco
736      !!----------------------------------------------------------------------
737      !!                  ***  ROUTINE zgr_zco  ***
738      !!
739      !! ** Purpose :   define the z-coordinate system
740      !!
741      !! ** Method  :   set 3D coord. arrays to reference 1D array
742      !!----------------------------------------------------------------------
743      INTEGER  ::   jk
744      !!----------------------------------------------------------------------
745      !
746      DO jk = 1, jpk
747         fsdept(:,:,jk) = gdept_0(jk)
748         fsdepw(:,:,jk) = gdepw_0(jk)
749         fsde3w(:,:,jk) = gdepw_0(jk)
750         fse3t (:,:,jk) = e3t_0(jk)
751         fse3u (:,:,jk) = e3t_0(jk)
752         fse3v (:,:,jk) = e3t_0(jk)
753         fse3f (:,:,jk) = e3t_0(jk)
754         fse3w (:,:,jk) = e3w_0(jk)
755         fse3uw(:,:,jk) = e3w_0(jk)
756         fse3vw(:,:,jk) = e3w_0(jk)
757      END DO
758      !
759   END SUBROUTINE zgr_zco
760
761
762   SUBROUTINE zgr_zps
763      !!----------------------------------------------------------------------
764      !!                  ***  ROUTINE zgr_zps  ***
765      !!                     
766      !! ** Purpose :   the depth and vertical scale factor in partial step
767      !!      z-coordinate case
768      !!
769      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
770      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
771      !!      a partial step representation of bottom topography.
772      !!
773      !!        The reference depth of model levels is defined from an analytical
774      !!      function the derivative of which gives the reference vertical
775      !!      scale factors.
776      !!        From  depth and scale factors reference, we compute there new value
777      !!      with partial steps  on 3d arrays ( i, j, k ).
778      !!
779      !!              w-level: gdepw(i,j,k)  = fsdep(k)
780      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
781      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
782      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
783      !!
784      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
785      !!      we find the mbathy index of the depth at each grid point.
786      !!      This leads us to three cases:
787      !!
788      !!              - bathy = 0 => mbathy = 0
789      !!              - 1 < mbathy < jpkm1   
790      !!              - bathy > gdepw(jpk) => mbathy = jpkm1 
791      !!
792      !!        Then, for each case, we find the new depth at t- and w- levels
793      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
794      !!      and f-points.
795      !!
796      !!        This routine is given as an example, it must be modified
797      !!      following the user s desiderata. nevertheless, the output as
798      !!      well as the way to compute the model levels and scale factors
799      !!      must be respected in order to insure second order accuracy
800      !!      schemes.
801      !!
802      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
803      !!         - - - - - - -   gdept, gdepw and e3. are positives
804      !!     
805      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
806      !!----------------------------------------------------------------------
807      INTEGER  ::   ji, jj, jk       ! dummy loop indices
808      INTEGER  ::   ik, it           ! temporary integers
809      LOGICAL  ::   ll_print         ! Allow  control print for debugging
810      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
811      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
812      REAL(wp) ::   zmax             ! Maximum depth
813      REAL(wp) ::   zdiff            ! temporary scalar
814      REAL(wp) ::   zrefdep          ! temporary scalar
815      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace
816      !!---------------------------------------------------------------------
817
818      IF(lwp) WRITE(numout,*)
819      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
820      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
821      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
822
823      ll_print = .FALSE.                   ! Local variable for debugging
824!!    ll_print = .TRUE.
825     
826      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
827         WRITE(numout,*)
828         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
829         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
830      ENDIF
831
832
833      ! bathymetry in level (from bathy_meter)
834      ! ===================
835      zmax = gdepw_0(jpk) + e3t_0(jpk)          ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
836      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
837      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
838      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
839      END WHERE
840
841      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
842      ! find the number of ocean levels such that the last level thickness
843      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
844      ! e3t_0 is the reference level thickness
845      DO jk = jpkm1, 1, -1
846         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
847         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
848      END DO
849
850      ! Scale factors and depth at T- and W-points
851      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
852         gdept(:,:,jk) = gdept_0(jk)
853         gdepw(:,:,jk) = gdepw_0(jk)
854         e3t  (:,:,jk) = e3t_0  (jk)
855         e3w  (:,:,jk) = e3w_0  (jk)
856      END DO
857      !
858      DO jj = 1, jpj
859         DO ji = 1, jpi
860            ik = mbathy(ji,jj)
861            IF( ik > 0 ) THEN               ! ocean point only
862               ! max ocean level case
863               IF( ik == jpkm1 ) THEN
864                  zdepwp = bathy(ji,jj)
865                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
866                  ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
867                  e3t(ji,jj,ik  ) = ze3tp
868                  e3t(ji,jj,ik+1) = ze3tp
869                  e3w(ji,jj,ik  ) = ze3wp
870                  e3w(ji,jj,ik+1) = ze3tp
871                  gdepw(ji,jj,ik+1) = zdepwp
872                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
873                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
874                  !
875               ELSE                         ! standard case
876                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
877                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
878                  ENDIF
879!gm Bug?  check the gdepw_0
880                  !       ... on ik
881                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
882                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
883                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
884                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
885                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
886                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) )   &
887                     &                     * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
888                  !       ... on ik+1
889                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
890                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
891                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
892               ENDIF
893            ENDIF
894         END DO
895      END DO
896      !
897      it = 0
898      DO jj = 1, jpj
899         DO ji = 1, jpi
900            ik = mbathy(ji,jj)
901            IF( ik > 0 ) THEN               ! ocean point only
902               e3tp (ji,jj) = e3t(ji,jj,ik  )
903               e3wp (ji,jj) = e3w(ji,jj,ik  )
904               ! test
905               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
906               IF( zdiff <= 0._wp .AND. lwp ) THEN
907                  it = it + 1
908                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
909                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
910                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
911                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
912               ENDIF
913            ENDIF
914         END DO
915      END DO
916
917      ! Scale factors and depth at U-, V-, UW and VW-points
918      DO jk = 1, jpk                        ! initialisation to z-scale factors
919         e3u (:,:,jk) = e3t_0(jk)
920         e3v (:,:,jk) = e3t_0(jk)
921         e3uw(:,:,jk) = e3w_0(jk)
922         e3vw(:,:,jk) = e3w_0(jk)
923      END DO
924      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
925         DO jj = 1, jpjm1
926            DO ji = 1, fs_jpim1   ! vector opt.
927               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
928               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
929               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
930               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
931            END DO
932         END DO
933      END DO
934      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
935      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
936      !
937      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
938         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_0(jk)
939         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_0(jk)
940         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_0(jk)
941         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_0(jk)
942      END DO
943     
944      ! Scale factor at F-point
945      DO jk = 1, jpk                        ! initialisation to z-scale factors
946         e3f(:,:,jk) = e3t_0(jk)
947      END DO
948      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
949         DO jj = 1, jpjm1
950            DO ji = 1, fs_jpim1   ! vector opt.
951               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
952            END DO
953         END DO
954      END DO
955      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
956      !
957      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
958         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_0(jk)
959      END DO
960!!gm  bug ? :  must be a do loop with mj0,mj1
961      !
962      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
963      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
964      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
965      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
966      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
967
968      ! Control of the sign
969      IF( MINVAL( e3t  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
970      IF( MINVAL( e3w  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
971      IF( MINVAL( gdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
972      IF( MINVAL( gdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
973     
974      ! Compute gdep3w (vertical sum of e3w)
975      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
976      DO jk = 2, jpk
977         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
978      END DO
979       
980      !                                               ! ================= !
981      IF(lwp .AND. ll_print) THEN                     !   Control print   !
982         !                                            ! ================= !
983         DO jj = 1,jpj
984            DO ji = 1, jpi
985               ik = MAX( mbathy(ji,jj), 1 )
986               zprt(ji,jj,1) = e3t   (ji,jj,ik)
987               zprt(ji,jj,2) = e3w   (ji,jj,ik)
988               zprt(ji,jj,3) = e3u   (ji,jj,ik)
989               zprt(ji,jj,4) = e3v   (ji,jj,ik)
990               zprt(ji,jj,5) = e3f   (ji,jj,ik)
991               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
992            END DO
993         END DO
994         WRITE(numout,*)
995         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
996         WRITE(numout,*)
997         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
998         WRITE(numout,*)
999         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1000         WRITE(numout,*)
1001         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1002         WRITE(numout,*)
1003         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1004         WRITE(numout,*)
1005         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1006      ENDIF 
1007      !
1008   END SUBROUTINE zgr_zps
1009
1010
1011   FUNCTION fssig( pk ) RESULT( pf )
1012      !!----------------------------------------------------------------------
1013      !!                 ***  ROUTINE eos_init  ***
1014      !!       
1015      !! ** Purpose :   provide the analytical function in s-coordinate
1016      !!         
1017      !! ** Method  :   the function provide the non-dimensional position of
1018      !!                T and W (i.e. between 0 and 1)
1019      !!                T-points at integer values (between 1 and jpk)
1020      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1021      !!----------------------------------------------------------------------
1022      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1023      REAL(wp)             ::   pf   ! sigma value
1024      !!----------------------------------------------------------------------
1025      !
1026      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1027         &     - TANH( rn_thetb * rn_theta                                )  )   &
1028         & * (   COSH( rn_theta                           )                      &
1029         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1030         & / ( 2._wp * SINH( rn_theta ) )
1031      !
1032   END FUNCTION fssig
1033
1034
1035   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1036      !!----------------------------------------------------------------------
1037      !!                 ***  ROUTINE eos_init  ***
1038      !!
1039      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1040      !!
1041      !! ** Method  :   the function provides the non-dimensional position of
1042      !!                T and W (i.e. between 0 and 1)
1043      !!                T-points at integer values (between 1 and jpk)
1044      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1045      !!----------------------------------------------------------------------
1046      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1047      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1048      REAL(wp)             ::   pf1   ! sigma value
1049      !!----------------------------------------------------------------------
1050      !
1051      IF ( rn_theta == 0 ) then      ! uniform sigma
1052         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1053      ELSE                        ! stretched sigma
1054         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1055            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1056            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1057      ENDIF
1058      !
1059   END FUNCTION fssig1
1060
1061
1062   SUBROUTINE zgr_sco
1063      !!----------------------------------------------------------------------
1064      !!                  ***  ROUTINE zgr_sco  ***
1065      !!                     
1066      !! ** Purpose :   define the s-coordinate system
1067      !!
1068      !! ** Method  :   s-coordinate
1069      !!         The depth of model levels is defined as the product of an
1070      !!      analytical function by the local bathymetry, while the vertical
1071      !!      scale factors are defined as the product of the first derivative
1072      !!      of the analytical function by the bathymetry.
1073      !!      (this solution save memory as depth and scale factors are not
1074      !!      3d fields)
1075      !!          - Read bathymetry (in meters) at t-point and compute the
1076      !!         bathymetry at u-, v-, and f-points.
1077      !!            hbatu = mi( hbatt )
1078      !!            hbatv = mj( hbatt )
1079      !!            hbatf = mi( mj( hbatt ) )
1080      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
1081      !!         function and its derivative given as function.
1082      !!            gsigt(k) = fssig (k    )
1083      !!            gsigw(k) = fssig (k-0.5)
1084      !!            esigt(k) = fsdsig(k    )
1085      !!            esigw(k) = fsdsig(k-0.5)
1086      !!      This routine is given as an example, it must be modified
1087      !!      following the user s desiderata. nevertheless, the output as
1088      !!      well as the way to compute the model levels and scale factors
1089      !!      must be respected in order to insure second order a!!uracy
1090      !!      schemes.
1091      !!
1092      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1093      !!----------------------------------------------------------------------
1094      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1095      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1096      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars
1097      REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace
1098      REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
1099      !!
1100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigw3
1101      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigt3
1102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsi3w3
1103      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigt3
1104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigw3
1105      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtu3
1106      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtv3
1107      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtf3
1108      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwu3
1109      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwv3
1110      !!
1111      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1112      !!----------------------------------------------------------------------
1113
1114      REWIND( numnam )                        ! Read Namelist namzgr_sco : sigma-stretching parameters
1115      READ  ( numnam, namzgr_sco )
1116
1117      IF(lwp) THEN                            ! control print
1118         WRITE(numout,*)
1119         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1120         WRITE(numout,*) '~~~~~~~~~~~'
1121         WRITE(numout,*) '   Namelist namzgr_sco'
1122         WRITE(numout,*) '      sigma-stretching coeffs '
1123         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
1124         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
1125         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
1126         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
1127         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
1128         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
1129         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
1130         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
1131      ENDIF
1132
1133      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
1134      esigt3  = 0._wp   ;   esigw3  = 0._wp 
1135      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
1136      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
1137
1138      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1139      hifu(:,:) = rn_sbot_min
1140      hifv(:,:) = rn_sbot_min
1141      hiff(:,:) = rn_sbot_min
1142
1143      !                                        ! set maximum ocean depth
1144      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1145
1146      DO jj = 1, jpj
1147         DO ji = 1, jpi
1148           IF( bathy(ji,jj) > 0._wp ) THEN
1149              bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1150           ENDIF
1151         END DO
1152      END DO
1153      !                                        ! =============================
1154      !                                        ! Define the envelop bathymetry   (hbatt)
1155      !                                        ! =============================
1156      ! use r-value to create hybrid coordinates
1157      DO jj = 1, jpj
1158         DO ji = 1, jpi
1159            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1160         END DO
1161      END DO
1162      !
1163      ! Smooth the bathymetry (if required)
1164      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1165      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1166      !
1167      jl = 0
1168      zrmax = 1._wp
1169      !                                                     ! ================ !
1170      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
1171         !                                                  ! ================ !
1172         jl = jl + 1
1173         zrmax = 0._wp
1174         zmsk(:,:) = 0._wp
1175         DO jj = 1, nlcj
1176            DO ji = 1, nlci
1177               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1178               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1179               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1180               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1181               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1182               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1183               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
1184               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1185               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
1186            END DO
1187         END DO
1188         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1189         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1190         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
1191         DO jj = 1, nlcj
1192            DO ji = 1, nlci
1193                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1194            END DO
1195         END DO
1196         !
1197         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1198         !
1199         DO jj = 1, nlcj
1200            DO ji = 1, nlci
1201               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
1202               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
1203               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
1204               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
1205               IF( zmsk(ji,jj) == 1._wp ) THEN
1206                  ztmp(ji,jj) =   (                                                                                   &
1207             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
1208             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
1209             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
1210             &                    ) / (                                                                               &
1211             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
1212             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
1213             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
1214             &                        )
1215               ENDIF
1216            END DO
1217         END DO
1218         !
1219         DO jj = 1, nlcj
1220            DO ji = 1, nlci
1221               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1222            END DO
1223         END DO
1224         !
1225         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
1226         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
1227         DO jj = 1, nlcj
1228            DO ji = 1, nlci
1229               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
1230            END DO
1231         END DO
1232         !                                                  ! ================ !
1233      END DO                                                !     End loop     !
1234      !                                                     ! ================ !
1235      !
1236      !                                        ! envelop bathymetry saved in hbatt
1237      hbatt(:,:) = zenv(:,:) 
1238      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1239         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1240         DO jj = 1, jpj
1241            DO ji = 1, jpi
1242               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
1243               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1244            END DO
1245         END DO
1246      ENDIF
1247      !
1248      IF(lwp) THEN                             ! Control print
1249         WRITE(numout,*)
1250         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1251         WRITE(numout,*)
1252         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1253         IF( nprint == 1 )   THEN       
1254            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1255            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1256         ENDIF
1257      ENDIF
1258
1259      !                                        ! ==============================
1260      !                                        !   hbatu, hbatv, hbatf fields
1261      !                                        ! ==============================
1262      IF(lwp) THEN
1263         WRITE(numout,*)
1264         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1265      ENDIF
1266      hbatu(:,:) = rn_sbot_min
1267      hbatv(:,:) = rn_sbot_min
1268      hbatf(:,:) = rn_sbot_min
1269      DO jj = 1, jpjm1
1270        DO ji = 1, jpim1   ! NO vector opt.
1271           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1272           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1273           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1274              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1275        END DO
1276      END DO
1277      !
1278      ! Apply lateral boundary condition
1279!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1280      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1281      DO jj = 1, jpj
1282         DO ji = 1, jpi
1283            IF( hbatu(ji,jj) == 0._wp ) THEN
1284               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1285               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1286            ENDIF
1287         END DO
1288      END DO
1289      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1290      DO jj = 1, jpj
1291         DO ji = 1, jpi
1292            IF( hbatv(ji,jj) == 0._wp ) THEN
1293               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1294               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1295            ENDIF
1296         END DO
1297      END DO
1298      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1299      DO jj = 1, jpj
1300         DO ji = 1, jpi
1301            IF( hbatf(ji,jj) == 0._wp ) THEN
1302               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1303               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1304            ENDIF
1305         END DO
1306      END DO
1307
1308!!bug:  key_helsinki a verifer
1309      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1310      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1311      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1312      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1313
1314      IF( nprint == 1 .AND. lwp )   THEN
1315         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1316            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1317         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1318            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1319         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1320            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1321         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1322            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1323      ENDIF
1324!! helsinki
1325
1326      !                                            ! =======================
1327      !                                            !   s-ccordinate fields     (gdep., e3.)
1328      !                                            ! =======================
1329      !
1330      ! non-dimensional "sigma" for model level depth at w- and t-levels
1331
1332      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
1333         !                         ! below rn_hc, with uniform sigma in shallower waters
1334         DO ji = 1, jpi
1335            DO jj = 1, jpj
1336
1337               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1338                  DO jk = 1, jpk
1339                     gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1340                     gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1341                  END DO
1342               ELSE ! shallow water, uniform sigma
1343                  DO jk = 1, jpk
1344                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1345                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1346                  END DO
1347               ENDIF
1348               IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1349               !
1350               DO jk = 1, jpkm1
1351                  esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1352                  esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1353               END DO
1354               esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
1355               esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
1356               !
1357               ! Coefficients for vertical depth as the sum of e3w scale factors
1358               gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
1359               DO jk = 2, jpk
1360                  gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1361               END DO
1362               !
1363               DO jk = 1, jpk
1364                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1365                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1366                  gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1367                  gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1368                  gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
1369               END DO
1370               !
1371            END DO   ! for all jj's
1372         END DO    ! for all ji's
1373
1374         DO ji = 1, jpi
1375            DO jj = 1, jpj
1376               DO jk = 1, jpk
1377                  esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
1378                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1379                  esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
1380                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1381                  esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
1382                     &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
1383                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1384                  esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
1385                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1386                  esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
1387                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1388                  !
1389                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1390                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1391                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1392                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1393                  !
1394                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1395                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1396                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
1397               END DO
1398            END DO
1399         END DO
1400         !
1401      ELSE   ! not ln_s_sigma
1402         !
1403         DO jk = 1, jpk
1404           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1405           gsigt(jk) = -fssig( REAL(jk,wp)        )
1406         END DO
1407         IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1408         !
1409         ! Coefficients for vertical scale factors at w-, t- levels
1410!!gm bug :  define it from analytical function, not like juste bellow....
1411!!gm        or betteroffer the 2 possibilities....
1412         DO jk = 1, jpkm1
1413            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
1414            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1415         END DO
1416         esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
1417         esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
1418
1419!!gm  original form
1420!!org DO jk = 1, jpk
1421!!org    esigt(jk)=fsdsig( FLOAT(jk)     )
1422!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1423!!org END DO
1424!!gm
1425         !
1426         ! Coefficients for vertical depth as the sum of e3w scale factors
1427         gsi3w(1) = 0.5_wp * esigw(1)
1428         DO jk = 2, jpk
1429            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1430         END DO
1431!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1432         DO jk = 1, jpk
1433            zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1434            zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1435            gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
1436            gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
1437            gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
1438         END DO
1439!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1440         DO jj = 1, jpj
1441            DO ji = 1, jpi
1442               DO jk = 1, jpk
1443                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1444                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1445                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1446                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1447                 !
1448                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1449                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1450                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1451               END DO
1452            END DO
1453         END DO
1454         !
1455      ENDIF ! ln_s_sigma
1456
1457
1458      !
1459!!    H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
1460
1461      fsdept(:,:,:) = gdept (:,:,:)
1462      fsdepw(:,:,:) = gdepw (:,:,:)
1463      fsde3w(:,:,:) = gdep3w(:,:,:)
1464      fse3t (:,:,:) = e3t   (:,:,:)
1465      fse3u (:,:,:) = e3u   (:,:,:)
1466      fse3v (:,:,:) = e3v   (:,:,:)
1467      fse3f (:,:,:) = e3f   (:,:,:)
1468      fse3w (:,:,:) = e3w   (:,:,:)
1469      fse3uw(:,:,:) = e3uw  (:,:,:)
1470      fse3vw(:,:,:) = e3vw  (:,:,:)
1471!!
1472      ! HYBRID :
1473      DO jj = 1, jpj
1474         DO ji = 1, jpi
1475            DO jk = 1, jpkm1
1476               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1477               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
1478            END DO
1479         END DO
1480      END DO
1481      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1482         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1483
1484      !                                               ! =============
1485      IF(lwp) THEN                                    ! Control print
1486         !                                            ! =============
1487         WRITE(numout,*) 
1488         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1489         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
1490         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1491      ENDIF
1492      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1493         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
1494         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
1495            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
1496         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
1497            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
1498            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
1499            &                          ' w ', MINVAL( fse3w (:,:,:) )
1500
1501         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
1502            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
1503         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
1504            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
1505            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
1506            &                          ' w ', MAXVAL( fse3w (:,:,:) )
1507      ENDIF
1508      !
1509      IF(lwp) THEN                                  ! selected vertical profiles
1510         WRITE(numout,*)
1511         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1512         WRITE(numout,*) ' ~~~~~~  --------------------'
1513         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1514         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
1515            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1516         DO jj = mj0(20), mj1(20)
1517            DO ji = mi0(20), mi1(20)
1518               WRITE(numout,*)
1519               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1520               WRITE(numout,*) ' ~~~~~~  --------------------'
1521               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1522               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1523                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1524            END DO
1525         END DO
1526         DO jj = mj0(74), mj1(74)
1527            DO ji = mi0(100), mi1(100)
1528               WRITE(numout,*)
1529               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1530               WRITE(numout,*) ' ~~~~~~  --------------------'
1531               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1532               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1533                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1534            END DO
1535         END DO
1536      ENDIF
1537
1538!!gm bug?  no more necessary?  if ! defined key_helsinki
1539      DO jk = 1, jpk
1540         DO jj = 1, jpj
1541            DO ji = 1, jpi
1542               IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1543                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1544                  CALL ctl_stop( ctmp1 )
1545               ENDIF
1546               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1547                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1548                  CALL ctl_stop( ctmp1 )
1549               ENDIF
1550            END DO
1551         END DO
1552      END DO
1553!!gm bug    #endif
1554      !
1555   END SUBROUTINE zgr_sco
1556
1557
1558   !!======================================================================
1559END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.