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

Last change on this file since 2460 was 2460, checked in by gm, 13 years ago

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