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

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

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

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

First attempt to put dynamic allocation on the trunk

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