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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2695

Last change on this file since 2695 was 2695, checked in by rblod, 14 years ago

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