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

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

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

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

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